Question with interp3

Hi all,

I got an error using interp3 as:

%  Error using griddedInterpolant
%  The grid vectors do not define a grid of points that match the given values.

My code is as:

addpath('nctoolbox-1.1.3');
load('matlab.mat')
setup_nctoolbox;
grib = ncgeodataset('2018_Jan.grib');
grib.variables;
grib.size('10_metre_U_wind_component_surface');
U_wind1    = grib.geovariable('10_metre_U_wind_component_surface');
grib.dimensions('10_metre_U_wind_component_surface');
size(U_wind1);
class(U_wind1);
U_wind      = U_wind1.data(:, :, :);
U_wind      = permute(U_wind,[2 3 1]);
U_wind      = squeeze(double(U_wind));
U_wind           = flipud(U_wind);
skip             = length(U_wind)/2;
U_wind           = circshift(U_wind,[0,skip]);

coords_latitude  = -90:0.75:90;
coords_longitude = -180:0.75:179.25;
coords_lat       = coords_latitude';
coords_lon       = coords_longitude';

path_time_total_grib = path_time_total/6;

for i=1:size(path_lat_value,1)
    for j=1:size(path_lat_value,2)
            U_wind_route(i,j) = interp3(coords_lon,coords_lat,...
                coords_time,U_wind,path_lon_value(i,j),path_lat_value(i,j),path_time_total_grib(i,j));  
    end
end

I also attach the path_lon_value, path_lat_value, path_time_total as matlab.mat and Grib file as this google drive link https://drive.google.com/open?id=1ifuEnEnSThdhsq6FHbUXDfy-DyhhB_Ba

Thank you for your help

Answer 1

[meghannmarie]

First error I saw was you were using lat/lon vectors that were 0.75 degree resolution where the wind data was 1 degree resolution. I instead read the lat/lon from the file to make sure the data matched the coordinates.

Next, the data's longitude range was 0:360 where the path longitude rage was -180:180. So I did the modulus of the longitude path data by 360 to fix that.

Lastly, you do not want to loop through data. Just do interp3 on complete dataset, much faster.

Try this:

addpath('nctoolbox-1.1.3');
load('matlab.mat')
setup_nctoolbox;
grib = ncgeodataset('2018_Jan.grib');
U_wind1 = grib.geovariable('10_metre_U_wind_component_surface');
U_wind = U_wind1.data(:, :, :);
U_wind = permute(U_wind,[2 3 1]);

time = grib.geovariable('time');
grib_time = double(time.data(:));

lat = grib.geovariable('lat');
grib_lat = double(lat.data(:));

lon = grib.geovariable('lon');
grib_lon = double(lon.data(:));

path_time_total_grib = path_time_total/6;
path_lon_value_grib = mod(path_lon_value,360);
U_wind_route = interp3( grib_lon,grib_lat,grib_time, U_wind,...
                path_lon_value_grib,path_lat_value,path_time_total_grib);

[Trung Ngo]

Hi meghann,

Thank you so much for your answer, I forgot to put the rest of the converting code back. Since my long and lat dataset is as -180 to 180 and -90 to 90, respectively.

Do you know how to convert path_lon_value_grib from [0 360] to [-180 180]?

Should it be:

%Convert dataset from [0 360] back to [-180 180]
path_lon_value_180 = mod((path_lon_value_grib +180),360)-180

Also what is your recommendation on the interpolation method for 1 degree resolution at sea? I have several options and not sure which one to choose.

• 'linear' (default)
• 'nearest'
• 'cubic'
• 'spline'
• 'makima'

Because if a storm is happening in that area (let's say 1x1 area), I wonder what is the best interpolation method that representing the vessel's position at sea is best fitted in that scenario. Again, I am really appreciated for your help.

SIncerely,

[meghannmarie]

It looks like MATLAB has a unwrapto180 function, but Ive never used it. I usually just keep my original longitudes and make a separate variable with the converted ones. Otherwise it would just be:

% lon is longitude vector
lon_idx = lon > 180;
lon(lon_ idx) = lon(lon_idx) - 360;

I dont have the nctools installed here, but maybe the resolution is an attribute of the geovariable. Otherwise, just read in lat and lon and take a difference:

lon_resolution = unique(diff(lon));
lat_resolution = unique(diff(lat));

I like to use cubic interpolation.

[Trung Ngo]

Thank you for your help, I will use cubic interpolation then

[Trung Ngo]

Hi meghann, I also wonder should I convert the path_lon_value from [-90 90] to [90 -90] like the one from grib file start from [90 -90] insteaed? If it is the case, should I just change the positive value to negative value and vice versa ?

%Convert path time to 4-hour format
path_time_total_grib = path_time_total/6;
%Convert path long and lat
path_lon_value_grib = mod(path_lon_value,360);
path_lat_value_grib = -path_lat_value;

[meghannmarie]

The path values are just query points. It is just pulling the value from your interpolated grid made from the grib file at that specific query point. The order of the query points shouldn't matter.

For instance you could pull values from lats completely out of order (e.g. -90,-88,72,-18,35,80,-40) and it will give you values from the interpolated grid at those points.

Understand?

[Trung Ngo]

Thank you so much for your details answer, I have one last question trying to pull the value from longtitude between 359 and 360 degree. While reading longtitude from the grib file, they are from 0 to 359. I tried to extract weather criteria from the point between 359 and 360 and it returns a NaN value. Should I use this code instead

coords_longitude = 0:1:360;

Sorry for asking too many questions

[meghannmarie]

So your 0° value is your 360°. If you were to to make your longitude coordsto be 0:1:360, you would have to copy the data from 0° and concatonate it onto the end of the data so you longitudes would match your data. Then you would have interpolated values between 359 and 360.

grib_lon_0_360 = 0:1:360;
U_wind_0_360 = cat(1,U_wind,U_wind(1,:,:));

Or if you do not need to extract data between 180° and 181°, you can transform data to be between -180 and 180. If you do convert to this you would have to sort the data to make sure the longitude is increasing:

grib_lon = unwrapto180(grib_lon);
[grib_lon_sorted I] = sort(grib_lon);
U_wind_Sorted = U_wind(I,:,:); 

If you do this do not mod your path longitudes by 180 to keep them in the -180 to 180 range.

I am at work and do not have MATLAB on this computer to test the code above.

[Trung Ngo]

I tried

grib_lon_0_360 = 0:1:360;
U_wind_0_360 = cat(1,U_wind,U_wind(1,:,:));

However, it gave interp error, perhaps concatenate arrays should start with

U_wind(:,1,:) 

I am not sure if this is correct and why cat has to be 2 instead of 181

clearvars;
load('matlab.mat')
setup_nctoolbox;
grib = ncgeodataset('2018_Jan.grib');
U_wind1 = grib.geovariable('10_metre_U_wind_component_surface');
U_wind = U_wind1.data(:, :, :);
U_wind = permute(U_wind,[2 3 1]);
time = grib.geovariable('time');
grib_time = double(time.data(:));
lat = grib.geovariable('lat');
grib_lat = double(lat.data(:));
lon = grib.geovariable('lon');
grib_lon = double(lon.data(:));
path_time_total_grib = path_time_total/6;
path_lon_value_grib = mod(path_lon_value,360);
U_wind_route = interp3(grib_lon,grib_lat,grib_time, U_wind,...
                path_lon_value_grib,path_lat_value,path_time_total_grib);

grib_lon_0_360 = 0:1:360;
grib_lon_0_360 = grib_lon_0_360';
U_wind_0_360 = cat(2,U_wind,U_wind(:,1,:));
U_wind_route_0_360 = interp3(grib_lon_0_360,grib_lat,grib_time, U_wind_0_360,...
                path_lon_value_grib,path_lat_value,path_time_total_grib);

[meghannmarie]

That should be correct.

I didnt have the data to see that longitude is the second dimension. So use 2 as your first input into cat because you are concatonating along the second dimension, use U_wind as second because you are adding data to end of longitudes, and U_wind(:,1,:) as the third input because that is the 1st values for longitude (0°) which is same as 360°.

U_wind_0_360 = cat(2,U_wind,U_wind(:,1,:));

[Trung Ngo]

Thank you so much for your help !!!

Answer 2

[meghannmarie]

Your welcome!