How to grid data based on distance, rather than longitude and latitude?

If I want to grid my data into longitude and latitude grids over the whole planet, I can simply do this:

[X, Y] = meshgrid([20:1:380],[-90:1:90]);
Z = griddata(lon, lat, parameter, X, Y);

The issue is that the earth is a globe. Distance is squeezed towards the poles. How do I do the gridding based on distance?

My advisor has asked me to do the gridding based on 1550km * 740km.

Please help.

[Image Analyst]

You can't map rectangles over the spherical planet without gaps or overlapping. I don't see anyway around it. Ask him to show you a picture of where he wants these rectangles placed.

Answer 1

[Kelly Kearney]

You could calculate the geographic distance from your data points to each point in the grid (using distance.m in the Mapping Toolbox or any one of many 3rd-party toolboxes, like m_map or GSW). But then you need to use a different interpolation algorithm than griddata, since griddata (and scatteredInterpolant) assumes Euclidean distance between points in its calculations. Something like inverse distance weighting might be appropriate.

The other option would be to project your data into x/y space, using a projection that preserves distance pretty well. This is usually my solution when working with regional datasets, but I'm not sure if this could be easily translated to a global scale, since all global projections include some manner of distortion. Perhaps you could break the interpolation into smaller regions?

[Leon]

Thank you very much! I can use divagrid.

Would you please elaborate the details of doing so. Right now, I only see how to calculate the distance.

[Kelly Kearney]

I'm not familiar with divagrid... which of my suggestions are you trying to implement: distance-based interpolation or projection?

[Leon]

I'm not particular with either of them, as long as I can get it done.

Thanks.

[Kelly Kearney]

Using my inverse distance weighting function, you can get an interpolation based on real geographic distance:

lat = rand(1000,1)*180 - 90;
lon = rand(1000,1)*360 + 20;
parameter = rand(1000,1);

[X, Y] = meshgrid([20:1:380],[-90:1:90]);
Z1 = griddata(lon, lat, parameter, X, Y);
Z2 = invdistweight(parameter, 5, [X(:) Y(:)], [lon lat], 'geo'); 
Z2 = reshape(Z2, size(X));

subplot(2,1,1);
pcolor(X,Y,Z1);
shading flat;
hold on;
scatter(lon, lat, [], parameter, 'filled');
title('griddata, linear interpolation');

subplot(2,1,2);
pcolor(X,Y,Z2);
shading flat;
hold on;
scatter(lon, lat, [], parameter, 'filled');
title('inverse distance weighting');

set(findobj('type', 'scatter'), 'markeredgecolor', 'k');

 -

You can see that near the poles, the interpolated data looks more squished in the IDW case than with griddata, as it should, since points are actually closer in the longitudinal direction than the latitudinal one.

To get exactly what your advisor requested is a little more complicated, since you'll have to create an uneven output grid, with more points at the equator than polewards rows. Depending on the application, that may or may not be worth the trouble... the important aspect in my opinion is to do the gridding based on geographic distance, regardless of the resolution of your output grid.

Edit: Just noticed that the Git repo for the invdistweight function seems to have gone a bit awry (no code there at the moment). I'll try to get that fixed by the end of the day).

[Leon]

Thank you so much, Kelly! Unfortunately, I got the below error.

Also is there a way I can do similar distance based gridding by using other methods? I've been asked to use DIVA gridding.

Thank you.

??? No appropriate method, property, or field addParameter for class inputParser.

Error in ==> invdistweight at 65
p2.addParameter('dist', [],      @(x) validateattributes(x, {'numeric'}, {'2d'}));

Error in ==> Untitled2 at 11
Z2 = invdistweight(parameter, 5, [X(:) Y(:)], [lon lat], 'geo');

[Kelly Kearney]

What version of Matlab are you running? For some reason, Matlab decided to phase out the old addParamValue in favor of addParameter sometime around 2012b, though as far as I can tell the two methods are interchangeable. So if you're running an older version of Matlab, try replacing addParameter with addParamValue on that line.

[Leon]

Thank you! It is working now, but when I switch it to my data, my computer just become frozen. I have global data with up 10,000 rows of data, and I want to grid it on a 1*1 grid.

Looks like "invdistweight" is a whole new gridding method. Is there a way I can do distance gridding with the gridding method I prefer, in this case, DIVA?

Thanks.

[Kelly Kearney]

With that many data points, your interpoint distance matrix is going to be huge, so yeah, IDW is probably not your best bet there. (I've thought about adapting IPDM to allow for geopgraphic distances, and therefore allow my function to deal with larger datasets, but I haven't gotten around to it yet).

Do you actually have a Matlab function to do the DIVA gridding? That's a pretty complex algorithm, and the appropriate format for your input coordinates would depend on exactly how that function calculates distance internally. Honestly, you'd probably be much better off using Ocean Data View, or one of the other interfaces specifically designed for that type of gridding; I'm pretty sure those deal with all the geographic distance issues already.