CODEWITHSUNDEEP
matlabcoordinatesinterpolationnetcdflarge-data
Terranees
Terranees

Reputation: 143

Retrieving data on coordinates which or not on the data grid through interpolation

I'm using Matlab to read a large (NetCDF) data set with information about a magnetic field. The data set is a three-dimensional array of 318x562x554 and I can retrieve have three one-dimensional array (318x1, 562x1 and 554x1) with each axis values of the coordinates. I would like to know the magnetic field values on points that do not fit on the data set grid. These points are in this case trajectory coordinates of a spacecraft placed in a two-dimensional array (3xn,n depends on how many coordinates you have).

x = ncread(file,'X_axis');
y = ncread(file,'Y_axis');
z = ncread(file,'Z_axis');
Bx = ncread(file,'Bx');

[x2,y2,z2] = meshgrid(y,x,z);              

length = numel(interval_ET_5000);
Bx_intp = zeros(1,length);
for i = 1:length
    [xi,yi,zi] = meshgrid(position_MEX_Mars_5000(1,i),...
                          position_MEX_Mars_5000(2,i),...
                          position_MEX_Mars_5000(3,i));
    F = interp3(x2,y2,z2,Bx,xi,yi,zi);
    Bx_intp(i) = F;
end

I have tried many things that didn't even work. This 'works' but not correct because the values in Bx_intp are way to high. Also because of the doing coordinates one at the time in a for loop makes it very slow, a normal run is about 3500 coordinates.

So basicly what I am looking for is a reverse scatteredInterpolant. This function accepts random data points and you interpolate the values on a meshgrid. But now I have a regular grid and I want interpolation on random points.

Upvotes: 0

Views: 247

Answers (1)

Terranees
Terranees

Reputation: 143

Thanks for the tip Ashish Uthama! I got it working with the code below. For other people with the same problem. You need ndgrid instead of meshgrid for griddedInterpolant and the coordinates need to be monotonic increasing. 

x = ncread(file,'X_axis');
y = ncread(file,'Y_axis');
z = ncread(file,'Z_axis');
Bx = ncread(file,'Bx');

[x2,y2,z2] = ndgrid(x,y,z);

F = griddedInterpolant(x2,y2,z2,Bx,'linear','none');
Bx_intp = F(position_MEX_Mars_5000(1,i),...
            position_MEX_Mars_5000(2,i),...
            position_MEX_Mars_5000(3,i));

Upvotes: 0

Related Questions