Display data of a particular variable present in a netcdf file

Hello people!

I am a newbie in MATLAB who has just started to use for the past month or so am facing issue while learning to work with oceanic netcdf files for my project purpose.

I have a netcdf file as an output from my oceanographic simulation and I want to get a single ouput for a particular variable from my output file for a specific latitude longitude, time index and depth.

eg: I want to get the value of u momentum velocity from the variable u for latitude "-8.22426367031282" and longitude "115.309657931722" , level index = 1 (for my s_rho value) , and a specific oceantime say " 01-Apr-2008-15:00:00" or any such timestep, when I input all these values say u_vel[115.309657931722 ,-8.22426367031282, 1 , '01-Apr-2008-15:00:00'] , I hope to get an ouput like 0.25477523 (any units).

Here are my dimension of the output file and the variable:

Output file: Dimensions:

xi_rho = 100

xi_u = 99

xi_v = 100

xi_psi = 99

eta_rho = 50

eta_u = 50

eta_v = 49

eta_psi = 49

N = 15

s_rho = 15

s_w = 16

tracer = 1

boundary = 4

ocean_time = 7 (UNLIMITED)

u

Size: 99x50x15x273

Dimensions: xi_u,eta_u,s_rho,ocean_time

Datatype: single

Attributes:

long_name = 'u-momentum component'

units = 'meter second-1'

time = 'ocean_time'

grid = 'grid'

location = 'edge1'

coordinates = 'lon_u lat_u s_rho ocean_time'

v

Size: 100x49x15x7

Dimensions: xi_v,eta_v,s_rho,ocean_time

Datatype: single

Attributes:

long_name = 'v-momentum component'

units = 'meter second-1'

time = 'ocean_time'

grid = 'grid'

location = 'edge2'

coordinates = 'lon_v lat_v s_rho ocean_time'

field = 'v-velocity, scalar, series'

I wrote a code for an attempt but every time I execute my MATLAB crashes and also I am not sure if what I am doing will give me any result as well: If anybody please suggest me the right script for my intended job, that would be of great great help.

file = 'Lombok_roms_his.nc';
latitude = ncread(file,'lat_u');
longitude = ncread(file,'lon_u');
mask = ncread(file,'mask_rho');
time = ncread(file,'ocean_time');
u_vel = ncread(file,'u');
v_vel = ncread(file,'v');
%% this is the part which I try to exceute in command window%%
u_vel[115.309657931722 ,-8.22426367031282, 1 , '01-Apr-2008-15:00:00']
v_vel[115.309657931722 ,-8.22426367031282, 1 , '01-Apr-2008-15:00:00']

Answer 1

[Walter Roberson]

latitude is effectively "plaid" -- there are 50 columns and 50 different latitudes (to within round-off). But you cannot just use unique(): there are over 2800 "unique" latitudes, all of which are the same as other latitudes to within 4e-15 .

The fact that it is effectively "plaid" allows us to do interpolation over it to figure out where any particular latitude is represented.

Likewise, longitude is effectively plaid but in the other direction, with there being 99 rows and 99 distinct longitude values to within round-off error.

The latitude range in the file is -8.96036167119211 to -8.2089282952945.

The longitude range in the file is 115.175227369057 to 116.320809555245.

You are hoping for information about latitude 1.2417, longitude 103.8, depth 4... but that data simply does not exist in the file.

Is it correct that depth = 4 corresponds to the 4th entry along the s_rho dimension?

(My internal reference for this matter is T0098919)

[Shankhaneel Basak]

@Walter Roberson Sir, I am sorry there was a interpretation error in my objective statement. In reference to the matter ticket number T0098919, I am updating the actual question with correct cordinates. the value depth = 4, is for my 4th value in s_rho matrix of dimension 15x1. As you say there are 2800 unique latitudes and similarly many unique longitudes. Then how should I be moving forward with respect to getting the specific value of 'u' variable for the say , latitude "-8.22426367031282" and longitude "115.309657931722" , level index = 1 (for my s_rho value) , and a specific oceantime say " 01-Apr-2008-15:00:00" or any such timestep?

I tried to execute the previous attempt with these cordinates and the result is same, no output fom MATLAB 2019 keeps running infinitely. What should I be writing in the script in order to get my desired output?

I am expecting an output something like u = 0.25477523 when I pass the command as per my previous script attempt:

u_vel[115.309657931722 ,-8.22426367031282, 1 , '01-Apr-2008-15:00:00']

Kind regards

[Walter Roberson]

There are a lot of NaN in u_vel, presumably corresponding to places not selected by the mask variable, such as places on land for a sea dataset. Your target location row is 12 1/2, column 49, which is a portion of the file that is nan in all depths and times. What would you like to do in that case?

[Walter Roberson]

The ocean times stored in that nc file are based upon 'seconds since 2011-01-01 00:00:00' and since none of the entries are negative, none of them can be in 2008. The stored times are

   01-Jan-2020 00:00:00
   01-Jan-2020 00:30:00
   01-Jan-2020 01:00:00
   01-Jan-2020 01:30:00
   01-Jan-2020 02:00:00
   01-Jan-2020 02:30:00
   01-Jan-2020 03:00:00

All UTC.

Here is demonstration code. You probably want this made into a function.

Note that your desired date syntax '01-Apr-2008-15:00:00' does not use any of the standard formats. Before making a function out of this, a programmer would want to know what the permitted date formats for input are to be: adjustments can be made for formats known in advance, but having to be able to handle any arbitrary date format is difficult.

   file = 'Lombok_roms_his.nc';

%The latitude range in the file is -8.96036167119211 to -8.2089282952945.
%The longitude range in the file is 115.175227369057 to 116.320809555245.

latitude = ncread(file,'lat_u');
longitude = ncread(file,'lon_u');
mask = ncread(file,'mask_rho');
ocean_time = ncread(file,'ocean_time');
u_vel = ncread(file,'u');

unique_lat  = uniquetol(latitude);    %corresponds to columns
unique_long = uniquetol(longitude);   %corresponds to rows
caltime = seconds(ocean_time) + datetime('2011-01-01 00:00:00', 'TimeZone', 'UTC');

target_lat = -8.22426367031282;
target_long = 115.309657931722;
depth_idx = 1;
target_time = datetime('01-Jan-2020 01:00:00', 'TimeZone', 'UTC');

time_idx = round(interp1(caltime, 1:length(caltime), target_time));
if isnan(time_idx)
    error('target time is not in the range stored in the file');
end

u_slice = u_vel(:,:,depth_idx,time_idx);

target_u_value = interp2(unique_lat, unique_long, u_slice, target_lat, target_long);

if isnan(target_u_value)
    error('location is over land or otherwise not available');
else
   fprintf('predicted u_val is %g\n', target_u_value);
end

[Walter Roberson]

In some cases, the query location might be adjacent to an area over land. For example if the query latitude were -8.2242601 but 8.2242605 was nan and 8.224260 was defined, then interp2() in that form would use the nan to interpolate.

For systems like this there is always the question of where to define the boundaries as occuring. If the exact lat/long combination in the tables is on a rock but 1cm over it plunges deep, then what should be stored in the table? If data exists in one square but not the next square over, should you assume that data is valid to half-way toward the next square?

[Shankhaneel Basak]

@Walter Roberson sir, sorry for bugging you again. In reference to the matter ticket number T0098919, as an extended query, I would like to ask 3 different questions. As per your latest comment which says that if I were to pass -8.2242601 and it's nearest -8.224605 is corresponding a NaN value, then interpolation using interp2() would be done using the NaN value itself, could we define the method of interpolation such as cubic or spline in the syntax of interp2() as additional parameter inside interp2() function? Or will it not accept any method inside interp2() involving NaN's ?

My second query or rather the most important one is that, in cases where we know cordinates precisely upto 2 digits after decimal itself and not upto 10 or 14 digits as in this case, how would we get to know whether the value corresponding to it is a NaN or a value? I see an if statement is passed in the code but could we not know it beforehand? Like you identified in my query that my target location is full of NaNs at all depths and all time stamps. How to know at which precision of latitude and longitude my value is present and where it's a Nan? For example if I am supposed to find out the value at somewhere roughly [-8.23,115.3] and -8.22426367031282 and 115.309657931722 has a NaN in it, how shall we determine which combination of lat,lon in this precision would have a value?

My third query would be that it seems in my attempt I've used lat_u and lon_u in my ncread command. In your sample demo script, I see you also made use of lat_u and lon_u to get the latitude and longitude variable and then unique latitude and longitude variable. Would this not create an error in the portion of your demo script where these are passed on into the target_u_value variable as my 'u' variable has the dimensions of xi_u,eta_u,s_rho,ocean_time? Will it not cause any dimension anomaly?

u_slice = u_vel(:,:,depth_idx,time_idx);
target_u_value = interp2(unique_lat, unique_long, u_slice, target_lat, target_long);

I know my questions are absolutely silly and rather monotonic. You may chose not to answer them. I was going through NetCDF documentation in MATLAB but didn't find any suitable documentation to handle the kind of task like the one in my question. I thought netcdf.getVar could do the job but it showed bunch of errors. Is there any good resource where I can learn analysing the type of netcdf files as in my question and maybe you know work on getting a time series plot etc. ?

Thank you so much for your generous help.

Kind regards

-Neel

[Shankhaneel Basak]

Mr @Walter Roberson , I am sad to let you know(solely due to my lack of knowledge) that when I try to execute the script, I find that I am always shown the error case: location is over land or otherwise not available . I get these for any lat, lon combination that I use.I do realise that the value of latitude and longitude I am enetering must be having a NaN at that depth index and timestamp and everything is fine with your approach of script. But I am confused as to how to know the exact cordinate(lat,lon) that won't be returning to me a NaN value for u_slice variable.

It would be of great help if you could explain me how to find out which cordinate would have the value of u_slice. For example: If I needed to find out the u_slice value at somewhere roughly [-8.23,115.3] and my chosen -8.22426367031282 and 115.309657931722 has a NaN in it, how shall we determine which combination of lat,lon in this exact precision of decimal for my query would have value?

Kind regards

-Neel

[Walter Roberson]

file = 'Lombok_roms_his.nc';
%The latitude range in the file is -8.96036167119211 to -8.2089282952945.
%The longitude range in the file is 115.175227369057 to 116.320809555245.
latitude = ncread(file,'lat_u');
longitude = ncread(file,'lon_u');
depth_idx = 1;
time_idx = 2;
u_slice = u_vel(:,:,depth_idx,time_idx);
surf(longitude, latitude, u_slice, 'edgecolor','none');
xlabel('long'); ylabel('lat')

u_mask = ncread(file,'mask_u');
idx = find(u_mask);
scatter(longitude(idx),latitude(idx),'*')
xlabel('long'); ylabel('lat')

u_mask is 1 for each location where u_vel has valid data; everywhere that u_mask is 0, then u_vel is nan. So the above has * for each location where there is valid data.

[Walter Roberson]

If I needed to find out the u_slice value at somewhere roughly [-8.23,115.3]

Near 115.3, you have to go to about -8.58 to get valid data.

Near -8.23 you have to go near 115.5 to get valid data.

Look at a map such as https://www.google.com/maps/place/8%C2%B013'48.0%22S+115%C2%B018'00.0%22E/@-8.229542,115.1889255,9.65z/data=!4m5!3m4!1s0x0:0x0!8m2!3d-8.23!4d115.3 to see that your target is solidly on land in Bali.

The nearest coast to the north is about https://www.google.com/maps/place/Pantai+Pura+Dalem+-+Bondalem/@-8.1103747,115.3165848,18.08z/data=!4m5!3m4!1s0x2dd1ece57dc831c1:0x272b5dd988c1ba77!8m2!3d-8.1095006!4d115.3168137 {antai Pura Dalem, at -8.1104, 115.3166 . However, that is outside the data for this nc file.

[Walter Roberson]

u_mask = ncread(file,'mask_u');
idx = find(u_mask);
selected_long = longitude(idx);
selected_lat = latitude(idx);

target_lat = -8.22426367031282;
target_long = 115.309657931722;
dists = pdist2([target_long, target_lat], [selected_long, selected_lat]);

[~, minidx] = min(dists,[],2);
closest_lat = selected_lat(minidx);
closest_long = selected_long(minidx);

fprintf('closest defined point is at lat = %.15g, long = %.15g\n', closest_lat, closest_long)

closest defined point is at lat = -8.22426367031282, long = 115.537605407545

u_slice = u_vel(:,:,depth_idx,time_idx);
target_u_value = interp2(unique_lat, unique_long, u_slice, closest_lat, closest_long);
target_u_value

target_u_value =

-3.45410558114883e-06

[Shankhaneel Basak]

Thanks a lot Mr Walter Roberson. It was such a great help.

Kind regards

-Neel

[Walter Roberson]

Note that that logic about finding the nearest point, interpolates at that nearest point. You should refrain from doing that if you can interpolate where you are instead.

... And if your origin point is 45 km from the nearest point you have data for, then it isn't clear that knowing about the data over there is meaningful. Moving the point could be useful for cases where you are shifting from just near the coast, over to coast that is effectively "right there" if only you had more precise positions (and data).