Extracting part of grid from meshgrid

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???

[Sara]

Why do you use griddata? is the original DEM missing some values? Have you tried:

   profiler on
  ...your code
   profiler off
  profiler viewer

to check what is the slowest step?

How big is Data?

For the clipping, do you want a rectangular region or something else?

[sumant]

I used griddata because, that is what I found while googling, the DEM is in xyz format, and so needs to be converted to a grid before being able to plot it. I am very open to other suggestions. I did try your suggestion and it is indeed the griddata, which is the slowest step. The data is ~250 MB in xyz format with several null values denoted by -9999. I replaced them by Nan.

For the clipping, I am mostly interested in a rectangular region (or may be trapezoidal).

[Sara]

Have you tried repmat? if you have nx data repeated ny times (or the other way around), that function will create the 2d array for you without the interpolation step done by griddata.

[Chad Greene]

May I ask which data set you're using? If xyz data is scattered it likely needs to be gridded. But it's possible that your data set is already gridded.

Answer 1

[Chad Greene]

1. Instead of griddata you may want to consider gridfit because it allows much larger inputs and does a very nice job of it. An improved version of gridfit found here extends gridfit to bicubic interpolation, which may be more appropriate for a DEM.

2. You can use inpolygon to clip extents of data bound by x1,x2,...,xn. An example of using inpolygon to clip an Antarctic DEM to an ice shelf is shown in the pdf documentation of the Bedmap2 Toolbox for Matlab. You don't need the Mapping Toolbox for many of the tools in that toolbox.

3. You don't have the Mapping Toolbox, but I've heard nice things about a free Matlab toolbox called m_map. I have not used m_map; I'd be interested in reading a comparison of m_map and the official toolbox.

[sumant]

Chad: I tried installing Bedmap2 toolbox twice, but it does not allows me to - unless, I have a mapping toolbox. So, is the "inpolygon" accessible without installing the bedmap2 toolbox?

[Sara]

Why using inpolygon? Can't you create the full space of coordinates with meshgrid(xr,yr) and then check which points satisfy your conditions with e.g., find(xr > x1 & xr < x2 & ... and so on?

[Chad Greene]

Hi Sumant,

Sorry you had trouble with the Bedmap2 Toolbox. Someday I may make a version that does not require the Mapping Toolbox.

The inpolygon function does not require any special toolboxes. Follow the example on page 47 of the Bedmap2 toolbox documentation, but swap lat and lon for x and y.