Interpolating by latitude and depth

I am very new to Matlab - coming from an R background and i could use some help getting started.

We sampled at numerous 'stations' along a latitudinal gradient. At each station, we collected samples from 10-20 depths. So i have a dataset that has latitude, longitude, station number, depth, and a whole range of variables (nutrients etc.).

I need to do a 2D interpolation of the data. 1) To interpolate the nutrient data from the surface to 4000 m at 1m intervals and also 2) across the latitude -60 to -58 degrees. What i have done is:

data = readtable('mock_data.csv')

[val, q] = unique(data.nut1);

z_profile = [1:4000]; 
Xq = -60:0.01:-58;

if length(q) > 0
   Vq = interp2(data.latitude(q), data.depth(q), data.nut1(q), Xq, z_profile); 
 
end

But i get an error

Error using griddedInterpolant

Interpolation requires at least two sample points for each grid dimension.

Error in interp2>makegriddedinterp (line 226)

F = griddedInterpolant(varargin{:});

Error in interp2 (line 126)

F = makegriddedinterp({X, Y}, V, method,extrap);

I used the example from herehttps://au.mathworks.com/help/matlab/ref/interp2.html Vq = interp2(X,Y,V,Xq,Yq) but im wondering if this error is because i have a number of "stations' that much each be treated independently?

Ive added a mock dataset, if anyone could point me in the right direction that would be greatly appreciated!

Answer 1

[Star Strider]

The scatteredInterpolant function would be best for this problem.

Try this —

T1 = readtable('mock_data.csv')

T1 = 15×5 table
    Station    depth    Lat    Nut1     Nut2 
    _______    _____    ___    ____    ______

       1          5     -60      2      50.41
       1         10     -60      4     123.21
       1        100     -60      5     171.61
       3          5     -59      3      82.81
       3          9     -59      6     228.01
       3         20     -59      7     292.41
       3         50     -59      8     364.81
       3        100     -59      9     445.21
       3        150     -59     10     533.61
       4          2     -58      3      82.81
       4          5     -58      4     123.21
       4          8     -58      7     292.41
       4         20     -58      8     364.81
       4         50     -58     10     533.61
       4        100     -58     12     734.41

FNut1 = scatteredInterpolant(T1.Lat, T1.depth, T1.Nut1);

Xq = -60:0.01:-58;                                          % Latitude Vector
Yq = linspace(0, 150, numel(Xq));                           % Depth Vector

[Lat_mtx,Dpt_mtx] = ndgrid(Xq, Yq);
Nut1_mtx = FNut1(Lat_mtx, Dpt_mtx);

figure
surfc(Lat_mtx, Dpt_mtx, Nut1_mtx)
grid on
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 1 (mol/L)')

FNut1 = scatteredInterpolant(T1.Lat, T1.depth, T1.Nut2);

Xq = -60:0.01:-58;                                          % Latitude Vector
Yq = linspace(0, 150, numel(Xq));                           % Depth Vector

[Lat_mtx,Dpt_mtx] = ndgrid(Xq, Yq);
Nut1_mtx = FNut1(Lat_mtx, Dpt_mtx);

figure
surfc(Lat_mtx, Dpt_mtx, Nut1_mtx)
grid on
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 2 (mol/L)')

You will have to repeat this for each nutirent (or other variable), however that is not difficult.

.

[Lavenia]

Thank you!

I managed to use your example on my dataset. May i please know why you suggest scatteredInterpolant instead of interp2? and what does numel(Xq) mean?

Thank you!

[Star Strider]

My pleasure!

I use scatteredInterpolant because your data are scattered, not gridded. If they were gridded, other options (such as interp2) would be available. The griddata function might also be an option, however in my experience, scatteredInterpolant is more robust. It is also possible to use scatteredInterpolant if you want to create cross-sectional data for any sort of section, for example across different depths and different latitudes.

In creating ‘Yq’ I use numel(Xq), the number of elements in ‘Xq’ to be certain that the lengths are equal. That is not absolutely necessary for the ndgrid function, however it then produces square matrices.

I created my grid over all the stations, because there was not enough data in any one station (in the example file) to do the interpolation over all the depths. I note that the stations are on different latitudes, so separating them by station may not be necessary, since the latitudes would separate them automatically.

If my Answer helped you solve your problem, please Accept it!

.

[Lavenia]

Sorry, im still trying to wrap my head around this. Re scatter vs gridded -> Why is the data considered scattered?

When i think scattered, i think that its random sampling across depth, latitude etc. but when i think gridded, its samples (data) collected from specific station along a latitudinal gradient. So they are more fixed, rather than random. i think im getting this completely wrong though!

[Torsten]

"Gridded" in MATLAB means that you have data for your nutrients in all these points:

x = [2 5 8 9 20 50 100 150];
y = [-58 -59 -60];
[X,Y] = meshgrid(x,y);
plot(X,Y,'o b')

But your data on such a grid are incomplete.

M = [5	-60	2	50.41
10	-60	4	123.21
100	-60	5	171.61
5	-59	3	82.81
9	-59	6	228.01
20	-59	7	292.41
50	-59	8	364.81
100	-59	9	445.21
150	-59	10	533.61
2	-58	3	82.81
5	-58	4	123.21
8	-58	7	292.41
20	-58	8	364.81
50	-58	10	533.61
100	-58	12	734.41];
;
plot(M(:,1),M(:,2),'o b')

One way to cope with this is to first interpolate 1d for all latitudes in depth direction and the to use "interp2" when your data are completed.

This gives e.g. for latitude = -59 and nutrient 1:

xq = [2 5 8 9 20 50 100 150];
x = M(4:9,1);
nut1 = M(4:9,3);
nut1q = interp1(x,nut1,xq,'linear','extrap')

nut1q = 1×8

0.7500 3.0000 5.2500 6.0000 7.0000 8.0000 9.0000 10.0000

[Lavenia]

oh that makes sense! thank you Thorsten! In R, i was doing a two-step interpolation as well: datav2 <- data %>% group_by(LONGITUDE, LATITUDE, STATION) %>% do(do_interp(.))

I managed to adapt star striders example for my work but i realised that by defining Yq = linspace(0, 150, numel(Xq)); then all the data is inteprolated down to 150m for every location.

Is there a way this can be changed depending on station? For example, some stations the seafloor is at 100m and some stations its at 150m. In the data, the last depth should be taken as the final depth so we cannot interpolate any further.

i guess my question is, how can i adapt xq in your example to not have a fixed end point at 150m but instead be variable depending on each station?

And may i please know what you mean by (4:9,1) and (4:9,3). My R brain is telling me that numbers are going from 4 to 9 with an interval of 1 or 3...

x = M(4:9,1);

nut1 = M(4:9,3);

[Star Strider]

It is considered scattered because the data are not sampled regularly across a specific grid. Here, the sampling is irregular — not all samples are at the same grid points. I added two illustrations demonstrating the difference between grided and non-gridded data.

Example —

T1 = readtable('mock_data.csv');

figure
stem3(T1.Lat, T1.depth, T1.Nut1)
hold on
scatter3(T1.Lat, T1.depth, T1.Nut1, 50, 'filled')
hold off
grid on
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 1 (mol/L)')
title('mock\_data')
axis('padded')

N = 5;
x = linspace(min(T1.Lat), max(T1.Lat), N);
y = linspace(min(T1.depth), max(T1.depth), N);
[X,Y] = ndgrid(x, y);
Z = rand(size(X(:)));

figure
stem3(X(:), Y(:), Z(:))
hold on
scatter3(X(:), Y(:), Z(:), 50, 'filled')
plot3(X, Y, zeros(size(X)), '--k')
plot3(X', Y', zeros(size(X)), '--k')
hold off
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 1 (mol/L)')
title('Gridded Data Illustration 1')
axis('padded')

N = 5;
x = linspace(min(T1.Lat), max(T1.Lat), N);
y = logspace(log10(min(T1.depth)), log10(max(T1.depth)), N);
[X,Y] = ndgrid(x, y);
% Z = rand(size(X(:)));

figure
stem3(X(:), Y(:), Z(:))
hold on
scatter3(X(:), Y(:), Z(:), 50, 'filled')
plot3(X, Y, zeros(size(X)), '--k')
plot3(X', Y', zeros(size(X)), '--k')
hold off
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 1 (mol/L)')
title('Gridded Data Illustration 2')
axis('padded')

I added the dashed grid lines to both the gridded illustrations to provide clarity.

EDIT — (4 Feb 2024 at 14:36)

With respect to extrapolating, if you do not want it to extrapolate beyond the supplied data (for example for depths that do not exist for a particular station), you can set that to 'none' although in that instance, specifying an interpolation method is required.

This revises my original code to avoid extrapolation —

T1 = readtable('mock_data.csv')

T1 = 15×5 table
    Station    depth    Lat    Nut1     Nut2 
    _______    _____    ___    ____    ______

       1          5     -60      2      50.41
       1         10     -60      4     123.21
       1        100     -60      5     171.61
       3          5     -59      3      82.81
       3          9     -59      6     228.01
       3         20     -59      7     292.41
       3         50     -59      8     364.81
       3        100     -59      9     445.21
       3        150     -59     10     533.61
       4          2     -58      3      82.81
       4          5     -58      4     123.21
       4          8     -58      7     292.41
       4         20     -58      8     364.81
       4         50     -58     10     533.61
       4        100     -58     12     734.41

FNut1 = scatteredInterpolant(T1.Lat, T1.depth, T1.Nut1, 'linear', 'none');

Xq = -60:0.01:-58;                                          % Latitude Vector
Yq = linspace(0, 150, numel(Xq));                           % Depth Vector

[Lat_mtx,Dpt_mtx] = ndgrid(Xq, Yq);
Nut1_mtx = FNut1(Lat_mtx, Dpt_mtx);

figure
surfc(Lat_mtx, Dpt_mtx, Nut1_mtx)
grid on
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 1 (mol/L)')

FNut1 = scatteredInterpolant(T1.Lat, T1.depth, T1.Nut2, 'linear', 'none');

Xq = -60:0.01:-58;                                          % Latitude Vector
Yq = linspace(0, 150, numel(Xq));                           % Depth Vector

[Lat_mtx,Dpt_mtx] = ndgrid(Xq, Yq);
Nut1_mtx = FNut1(Lat_mtx, Dpt_mtx);

figure
surfc(Lat_mtx, Dpt_mtx, Nut1_mtx)
grid on
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 2 (mol/L)')

figure
surfc(Lat_mtx, Dpt_mtx, Nut1_mtx)
grid on
colormap(turbo)
xlabel('Latitude (°)')
ylabel('Depth (m)')
zlabel('Nutrient 2 (mol/L)')
view(140,20)

It would likely be necessary to rotate the plot to see the effects of preventing extrapolation.

.

[Torsten]

i guess my question is, how can i adapt xq in your example to not have a fixed end point at 150m but instead be variable depending on each station?

You can use xq as the array of values that end at the minimum of the maximum depths over all the stations (i.e.100 in this case). Thus you could remove the measurements at d = 150 in this case.

[Lavenia]

ok - i cant remove a measurement unfortunately because im mapping the nutrient concentration to the seafloor so i'd have to have the data interpolated to different maximum depths (i.e., the last depth reading is always the deepest, so it has to be inteprolated up to that point).

[Torsten]

ok - i cant remove a measurement unfortunately because im mapping the nutrient concentration to the seafloor

Thus back to "ScatteredInterpolant", I guess.

And may i please know what you mean by (4:9,1) and (4:9,3). My R brain is telling me that numbers are going from 4 to 9 with an interval of 1 or 3...

x = M(4:9,1);

nut1 = M(4:9,3);

It means that I extract the depths and nutrient1 data for -59 latitude from the data matrix M, namely M(4,1),M(5,1),...,M(9,1) for the depths and M(4,3),M(5,3),...,M(9,3) for nutrient1.

[Torsten]

@Star Strider

So it would be possible to use the original data, specify "none" for the grid points missing for complete "gridded data" and use "GriddedInterpolant" or "interp2" ?

[Lavenia]

Thanks star strider!

i tried your method but i get an error:

Warning: Contour not rendered for non-finite ZData

is that because Yq says 159? Yq = linspace(0, 150, numel(Xq));

[Star Strider]

@Torsten — I believe scatteredInterpolant will interpolate over missing grid points within the limits, and specifying 'none' for the extrapolation method will prevernt it from extrapolating beyond the given data.

@Lavenia — If any of the ‘Z’ data are non-finite (I believe that in this instance they would be NaN values, although they could be Inf), the function does not specifically deal with them and throws an error. (I cannot find anything in the documentation that specifically states what it does with NaN or Inf values.) If there are Inf values, I would first set them to NaN (see footnote) and then use fillmissing to to interpolate the NaN values. Using rmmissing is another option, however that removes the entire row of data containing a NaN value in any column.

Also, the More About documentation section explains the difference between gridded and scattered data in the context of the function.

Footnote: Example code for that would be —

Z = [1 2 Inf 4; 5 Inf 7 NaN; 9 10 11 12]

Z = 3×4

1 2 Inf 4 5 Inf 7 NaN 9 10 11 12

Z(~isfinite(Z)) = NaN

Z = 3×4

1 2 NaN 4 5 NaN 7 NaN 9 10 11 12

Z = fillmissing(Z, 'nearest')

Z = 3×4

1 2 7 4 5 10 7 12 9 10 11 12

Consider different interpolation methods to give the result you want.

.

[Lavenia]

Thank you so much Star Strider! i managed to get it to work with the code below. With this, i do not get the error: Warning: Contour not rendered for non-finite ZData

but, its interpolating well beyond its maximum depth. Is there a way i can interpolate only to the last depth for a given station? Instead of down to 4000m for everything?

And would you have any tips or keywords i can search for matlab specific, to smoothen the interpolation? It makes these really sharp contours. i tried changing from linear to nearest but it didnt make a difference.

And again, thank you so much!

  [val, q] = unique(data.Nut1)
   Yq = [1:1:4000];
   Xq = -65.6:0.1:-62.7;
   % Change the interpolation method from 'linear' to 'nearest'
    F.Method = 'nearest'
    Nut1 = scatteredInterpolant(data.LATITUDE(q), data.PRESSURE(q),    data.Nut1(q));
    [Lat_mtx,Dpt_mtx] = ndgrid(Xq, Yq);
    Nut1_mtx = Nut1(Lat_mtx, Dpt_mtx);
    figure; contourf(Lat_mtx, Dpt_mtx, Nut1_mtx)
     set(gca,'YDir','reverse')
     colorbar; %colormap(cmap);
     hold on
     scatter(data.LATITUDE, data.PRESSURE, 20, data.Nut1, "filled","black")
     cb = colorbar;
     ylabel('Depth (m)')
     ylabel(cb,'Nut1 (umol/kg)');
     colormap(turbo)
     xlabel('Latitude (°)')

[Star Strider]

My pleasure!

If it is interpolating beyond the maximum depth, it is extrapolating. If you do not want it to extrapolate, you have to tell it, and you do that (as I mentioned earlier in the EDIT at the end) by giving it an interpolation method of your choice (I chose 'linear' here) and then 'none' for an extrapolation method:

FNut1 = scatteredInterpolant(T1.Lat, T1.depth, T1.Nut1, 'linear', 'none');

and in your code example:

Nut1 = scatteredInterpolant(data.LATITUDE(q), data.PRESSURE(q), data.Nut1(q), 'linear', 'none');

That should stop it from extrapolating.

‘And would you have any tips or keywords i can search for matlab specific, to smoothen the interpolation?’

That may in part depend on the initial data distribution (more data are better), however you can choose a different interpolation method, specifically 'natural':

Nut1 = scatteredInterpolant(data.LATITUDE(q), data.PRESSURE(q), data.Nut1(q), 'natural', 'none');

That may improve the result.

.

[Lavenia]

Yes, i tried that, and also linear and nearest extrapolation method from here. It does cut the bottom depths where i dont have data, but it also weirdly cuts the sides where i do have data. :( i tried adding more latitude but that didnt work either.

[Star Strider]

I am just guessing since I don’t have your data, however starting ‘Yq’ from 0 instead of 1 may solve that —

Yq = 0:1:4000;

instead of:

Yq = [1:1:4000];

Also, the square brackets are not necessary and can slow your code.

Answer 2

[Walter Roberson]

interp2() expects the first two inputs to be either vectors or 2d arrays. If they are vectors then interp2() automatically expands them into grids.

interp2() expects the third input to be a 2D array.

It is not clear what q is in your code. If it is a scalar, then you would be asking to interp2() a scalar, a scalar, and a scalar -- not something that is suitable for interpolating over.

If your q is a vector of indices or is a logical mask, then latitude(q) and depth(q) are potentially vectors, so that part potentially works. But nut(q) would have to be a vector in that circumstance, not a 2d array.

The first two parameters are effectively there to describe the "marginal" indices of the grid of data that is the third parameter.

[Lavenia]

Hi Walter,

Thanks for the help! I edited the question to give a little more clarification.

The first two inputs (latitude and depth) are vectors. The third input is the nutrient concentration which varies with depth and by latitude.

q should be a vector as it just identifying unique nutrient values and preventing duplicates - but i may be wrong in how im doing this, so any help would be great!

[Walter Roberson]

The first two inputs (latitude and depth) are vectors. The third input is the nutrient concentration which varies with depth and by latitude.

For interp2() purposes, the third input must be length(depth) by length(latitude)

%example
latitude = linspace(0,1,10);
depth = linspace(0,1,15);
nc = rand(length(depth), length(latitude));
interp2(latitude, depth, nc, [0.15], [0.73])

ans = 0.5447