Which depth does any temperature have in the each point of 3D-field?

Hi,

I have a 3D matrix with a water temperature (X x Y x Z dimensions where X and Y are latitude and longitude, Z is a deep). Also, I have a vector with depths (1 x Z dimensions). I need to know which depth does any temperature (f.g., 3 °C) have in the each point of field.

Depths vector has a different discreteness:

depth = [1 5 10 25 50 100];

Temperature array does not always contain 3 °C value, f.g.:

 temp(1,1,:) = [5.3 5.0 4.3 3.3 2.7 2.5]);

As a result, I need a X x Y depth field at which the 3 °C temperature was observed.

I think I should use linear interpolation between relevant levels but I have no idea how exactly I can do it. It will be enough if you help me determine the depth at one point of field.

I hope anyone can help me,

Thanks,

Elena

[KSSV]

REad about slice.

[Elena Novoselova]

Thank you for your answer.

I read about it but this isn't exactly what I want. Anyway, I'll remember this function because it's usfull for my study.

Answer 1

[Guillaume]

You just need to create an interpolant, Z = F(X, Y, T) and then query it for all X and Y with T = 3. I'm assuming you have vectors for X and Y as well:

%X, Y, Z vectors of longitude and latitude and depth
%temperature: 3D matrix of temperatures
[Xfull, Yfull, Zfull] = ndgrid(X, Y, Z);  %creates 3D matrices of X, Y and Z coordinates
TempInterpolant = griddedInterpolant(Xfull, Yfull, temperatures, Zfull);  %creates Z = F(X, Y, T)

%now you can query the interpolant wherever you want:
[Xquery, Yquery] = ndgrid(X, Y);  %creates 2D arrays of X and Y
Tquery = 3;
ZatTquery = TempInterpolant(Xquery, Yquery, repmat(Tquery, size(Xquery)))

[Elena Novoselova]

Thank you for quick and clear answer. It's really easier than I thougth.

But I have one problem:

"Error using griddedInterpolant. The coordinates of the input points must be finite values; Inf and NaN are not permitted."

I think this is due to the fact that the temperature matrix has NaNs at several bottom levels because bottom flor isn't constant. I can't cut all levels which connect NaNs because sometimes I have 2000 m depth and sometimes only 50 m. Also, I have a land points. This's the second reason of NaNs.

Do you know what can I do with that?

[Guillaume]

You'd have to remove the NaN temperatures and the corresponding elements from the Xfull, Yfull, and Zfull arrays. As a result, you'd have to switch to a scatteredInterpolant since you would no longer have gridded data:

%X, Y, Z vectors of longitude and latitude and depth
%temperatures: 3D matrix of temperatures with maybe some NaNs
[Xfull, Yfull, Zfull] = ndgrid(X, Y, Z);  %creates 3D matrices of X, Y and Z coordinates
validtemp = ~isnan(temperatures);
TempInterpolant = scatteredInterpolant(Xfull(validtemp), Yfull(validtemp), temperatures(validtemp), Zfull(validtemp));

%Query same as before:
[Xquery, Yquery] = ndgrid(X, Y);  %creates 2D arrays of X and Y
Tquery = 3;
ZatTquery = TempInterpolant(Xquery, Yquery, repmat(Tquery, size(Xquery)))

[Elena Novoselova]

I have the last question. How can I return NaNs in the 'ZatTquery'?

I use this but maybe there's easier way to do it. I always use cycles but I know that it slows down the program significantly. So I would like to learn how to program it without cycles.

ZatTquery(ZatTquery<0) = NaN; % Delete if 'ZatTquery' is negative depth
ZatTquery(~validtemp(:,:,1)) = NaN; % Delete land 
for i = 1:length(lat) % Delete if 'ZatTquery' is more than real depth in point
    for j = 1:length(lon)
        depth_point = depth(squeeze(validtemp(i,j,:)));
        if ~isempty(depth_point) && ZatTquery(i,j) > depth_point(end) 
            ZatTquery(i,j) = NaN;
        end
    end
end

If you know I would be grateful for the help.

[Guillaume]

If I understood correctly, and assumed that depth(i, j, :) is ordered in increasing depth, your loop is equivalent to:

ZatTquery(Zatquery > max(depth, [], 3)) = NaN;

However, I think (untested) that you may be able to get the NaN directly if you specify 'none' as the extrapolation method when you create the interpolant:

TempInterpolant = scatteredInterpolant(Xfull(validtemp), Yfull(validtemp), temperatures(validtemp), Zfull(validtemp), 'linear', 'none');

[Elena Novoselova]

Thank you!