Meshgrid / Surface clipping

Question

I am working with some antarctic DEM data in Matlab. So, far I have been able to generate a nice looking mesh, with the following basic code:

load (Data.xyz)
X = Data(:,1);
Y = Data(:,2);
Z = Data(:,3);
xr = unique(X);
yr = unique(Y);
gz = zeros(length(yr),length(xr));
gz = griddata(X,Y,Z,xr,yr');

figure
mesh(xr,yr,gz);
hold on
contour3(xr,yr,gz,'k-');
hold off

Now I have a few questions, which I have not been able to answer despite being at it since past couple of days and looking at all forums and googling day and night. I hope you all experts might be able to suggest me something. My questions are:

1. The above code takes a lot of time. Agreed that the DEM for antarctica is large sized and slow response time for a code does not necessarily mean that its incorrect. However, I am totally unable to run this code on my laptop (2.5 GHz/4GB) - its so slow. I am wondering if there are other ways to generate mesh which are faster and more efficient.

2. The second issue is that the above "Data.xyz" contains DEM data from all antarctica. After generating a mesh, I want to clip it based on locations. Say, for e.g., I want to extract mesh data for area bound by x1,y1, x2,y2, x3, y3, and x4,y4. How do I go about doing that? I could not find a suitable function or tool or any user script anywhere which will allow me to do so. Is it possible to cut a mesh in matlab?

I am running Matlab 2012a, and I do not have access to mapping toolbox. Any suggestions???

Answer 1

1.I'd just like to clarify what it is you want your code to do. For the data X,Y,Z, I would assume these are points (X,Y) (not sampled on a grid) each associated with some elevation Z.

When you call

gz = griddata(X,Y,Z,xr,yr');

You are saying that every possible pair (xr(i),yr(j)) are the locations on your grid where you want to sample the surface to create your mesh. I could be wrong but I don't think this is what you want? I think you would rather sample at equally spaced points, instead using something like...

xr = min(X):spacing:max(X);
yr = min(Y):spacing:max(Y);   
gz = griddata(X,Y,Z,xr,yr');  % ' is for the transpose
mesh(xr,yr,gz);

Where spacing a reasonable number for the scale of your data. The way your code is now, you could be taking far more samples than you want, which could be why your code takes so long.

For 2. I think you could get a for loop to just add the points insde the region of interest to three new lists of X,Y,Z values. If your region is rectangular bounded by x_left,x_right and y_left,y_right ...

Xnew = []; Ynew = []; Znew = [];
for i = 1:length(X)
  if ( x_left<X(i) )&&( X(i)<x_right )&&( y_left<Y(i) )&&( Y(i)<y_right )
    Xew = [Xnew, X(i)];
    Ynew = [Ynew, Y(i)];
    Znew = [Znew, Z(i)];
  end
end