How to do interpolation on the map?

Hello all,

In the attached image, the blue circles are weather stations within the Alaska, and red cross points are the gerenrated mesh. I have estimated a parameter called "frost depth" in all weather stations. I want to know how I can do interpolation process to get frost depth values in the mesh grid points. Can anyone help me with this regard?

[Walter Roberson]

The grid points appear to mesh along lines of equal latitude or equal longitude, with there being a wide enough span of latitudes that the curvature of the Earth can make a noticeable difference. I estimate the bottom-most lines of longitude are roughly twice as far apart as would be the case at the top.

The question then becomes which coordinate system you wish to interpolate in. The normal bilinear interpolation performed by interp2() assumes that there is an implied quadralateral that is in linear coordinates, at least to a close-enough approximation. That, for example, 1/10 degree coordinate diffence has the same physical meaning near the top as near the bottom. That is probably not a good enough approximation for your purposes.

[Behrooz Daneshian]

I want to do interpolation in lat and lon coordinate system since I have lat and lon of all weather stations as well as their corresponding values of frost depth.

Answer 1

[Nathan Hardenberg]

You can interpolate the datapoints with the griddata() function. As @Walter Roberson mentioned this also assumes a uniform grid, which is not given when working with world-coordinates. I personally think this is not such a big problem here, since your datapoints are already an estimated parameter and the rest is a simple interpolation.

What you get is a "better" resolution when going further north and your interpolation goes along the latitude lines instead of going the shortest way from point to point. But in the end you have to decide/know if this is good enough for you.

lon = -1*[165, 155, 155, 140, 132, 140, 145, 180, 165]; % data
lat = [55, 60, 70, 67, 62, 60, 70, 52, 65];
depth = [3, 1, 5, 2, 0.5, 3, 4, 0, 1];

[x,y] = meshgrid(-1*(130:180), 50:75);  % your grid

v = griddata(lon, lat, depth, x, y, "natural");  % interpolate
% possible interpolations: "linear", "nearest", "natural", "cubic" and "v4"
% v is now already your interpolated depth

% Visualization
figure(1);
surf(x,y,v)
title('Surf-Plot Depth (as uniform grid)');

x_ = reshape(x, 1, []); y_ = reshape(y, 1, []); v_ = reshape(v, 1, []);  % flatten

figure(2);  % plot depth on map (probably not the best way to visualize)
geodensityplot(y_, x_, v_,"FaceColor","interp")
geobasemap topographic
geolimits([50 71],-[179 130])
title("depth on worldmap")