Surface bounded by a polygon

I have a polygon and a surface that covers only part of this polygon (this is in 2-D x and y axis). Is there any way to delimit the surface by the polygon? I want only the part that is inside the polygon.

Thanks.

Answer 1

[Chad Greene]

Without any sample data or illustrations it's mighty difficult to know exactly what you mean, but it sounds like inpolygon is the function you're looking for. Just feed it the x,y locations of your surface (generated by meshgrid) and it will tell you which corresponding surface locations are inside the polygon.

[Sarah Andrade]

Hi Chad! This is an example of what I'm trying to say (my code is too big):

lat = [-22.738889 
-23.462778
-23.433889
-21.595833
-22.905833]

long = [-45.590833
-46.532778
-45.070833
-46.888889
-47.060833]

variavel = [5793.4
7857.8
11163.7
4893.9
6879.6]

rangeZ = floor(min(long)):.01:ceil(max(long));
rangeX = floor(min(lat)):0.01:ceil(max(lat));
[X,Z] = meshgrid(rangeX,rangeZ);
Y = griddata(lat,long,variavel,X,Z,'cubic');

lat2 = [-23.4
-23
-22
-23.4]

long2 = [-47
-45.2
-46.2
-47]

figure(11)
surface(X,Z,Y)
axis tight
shading interp
colorbar
colormap(jet)
hold on
plot(lat2,long2)

I want to cut the surface inside the triangle and stay with it only. Using 'inpolygon' is there any way to erase the surface outside the polygon? I didn't find it before.

Thank you!

[Chad Greene]

Hi Sarah,

Thanks for providing a minimal working example. That helps tremendously. Before we go any further, I want to make a fundamental change to your code. For this kind of stuff we often think of longitudes as being like x values and latitudes as similar to y values. However, the use of variable names X and Y can imply that the spatial dimensions have been projected into meters, so it gets confusing to mix lat, long, X, and Y all on the same map with the same units. If I were faced with your problem I'd start by plotting the raw scattered data in unprojected coordinates, but letting long and lat be the horizontal and vertical axes of a plot, respectively (note, this is a change from your method):

figure
scatter(long,lat,60,variavel,'filled'); 
hold on
plot(long2,lat2)
xlabel 'longitude'
ylabel 'latitude'

I've made a habit of typically using lowercase variable names (lat, long, x, y) for scattered data, and uppwercase (LAT, LONG, X, Y) for gridded data. That's a bit more intuitive for me than letting gridded longitudes become Z and latitudes X. So I'd rewrite your gridding section like this:

% No need to change variable names, just change their cases: 
[LONG,LAT] = meshgrid(floor(min(long)):.01:ceil(max(long)),...
   floor(min(lat)):0.01:ceil(max(lat)));

VARIAVEL = griddata(long,lat,variavel,LONG,LAT,'cubic');

s = surface(LONG,LAT,VARIAVEL); 
shading interp

Now use inpolygon to find which grid cells are inside the triangle:

% These points are inside the polygon: 
in = inpolygon(LONG,LAT,long2,lat2); 

% Set points not inside (~in) the polygon to NaN: 
VARIAVEL(~in) = NaN; 

% delete the surface we just plotted: 
delete(s) 

% Plot the new surface: 
s = surface(LONG,LAT,VARIAVEL); 
shading interp

[Sarah Andrade]

Thank you so much, Chad!!

[Daniel Faller]

Thank very much, Chad! I really had some use for your answer!