I think you want inpolygon. To get the indices of all data inside the polygon 30 to 40 degrees north and 20 to 10 degrees west, you'd do
ind = inpolygon(lon,lat,[-20 -10],[30 40]);
Then your data measurements in that polygon are: vals(ind)
There's a Mapping Toolbox function called ingeoquad that does the same thing, but inpolygon does it just as well and doesn't require the Mapping Toolbox.
