NetCDF file grid box extraction

I'm new to Matlab programming (usually use IDL) but I've been tasked with extracting meteorological data from netcdf files at specific gridsquares. I've looked at the structure of the file with ncread but my problem is that the netcdf is on a 193x by 130y grid that does not match any global grid. The gridpoint I'm looking for is at: lat 28.046, lon 86.902 in decimal degrees. How do I find where this is in the netcdf file?

An example of one of the netcdf files (test.nc) can be found here: https://drive.google.com/folderview?id=0B-amO6HuYM_yfjhFT3djOXF5dmhwaWc5b3Q4cXZoVkhMVFdJeU8xZ20tZXVpcXZTdzlVTEU&usp=sharing

Thanks

Answer 1

[Chad Greene]

I think the issue is that the lat/lon grid is not regular--they're in rotated coordinates. If you have the Mapping Toolbox, here's the first "slice" of tas:

lati = 28.046;
loni = 86.902;
tas = ncread('test.nc','tas'); 
rlat = ncread('test.nc','rlat'); 
rlon = ncread('test.nc','rlon'); 

origin = newpole(80,-123);

[rlati,rloni] = rotatem(lati,loni,origin,'forward','degrees'); 

tas1 = squeeze(tas(:,:,1)); 

interp2(rlat,rlon,tas1,rlati,rloni)
ans =
        261.41

[Chad Greene]

If you do not have the Mapping Toolbox and this is the only location you're interested in, the rotate coordinates are (19.2798, 27.7826).

rlati = 19.2798;
rloni = 27.7826;
tas = ncread('test.nc','tas'); 
rlat = ncread('test.nc','rlat'); 
rlon = ncread('test.nc','rlon'); 
t = datenum(1949,12,ncread('test.nc','time')); 

tasi = NaN(size(t)); 
for k = 1:length(t)
    tasi(k) = interp2(rlat,rlon,tas(:,:,k),rlati,rloni); 
end

plot(t,tasi,'linewidth',2)
box off
datetick('x','mmm ''yy')
axis tight

[James Douglas]

Thanks for your help, this seems to work well. The final thing I need to is to print to a file with TAS in one column and time in 2 other columns. However, time is currently 'days since 1949-12-00' whereas I need it to have 'year' in one column and 'day of year' in the next. Are there any inbuilt matlab function that can do this?

ideally, I'd like to implement this so that it can work for several netcdf files. Foe example, the 'test.nc' file I provided is only the year 1980 and I have several other files from 1981-1985, 1986-1990, etc. to 2006-2010. Eventually, I need one file containing TAS for all these years, year, and day of year.

Any ideas or help is very much appreciated, thanks!

[James Douglas]

i should add that my code at the moment is as follows:

lati = 28.046;
loni = 86.902;

%filename
fn = 'test.nc'

%COORDINATE ROTATION
tas = ncread(fn,'tas'); 
rlat = ncread(fn,'rlat'); 
rlon = ncread(fn,'rlon'); 

origin = newpole(80,-123);
[rlati,rloni] = rotatem(lati,loni,origin,'forward','degrees'); 

[rlati,rloni] %print

t = datenum(1949,12,ncread(fn,'time')); 
tasi = NaN(size(t)); 

for k = 1:length(t)
    tasi(k) = interp2(rlat,rlon,tas(:,:,k),rlati,rloni); 
end

tasic=tasi-273.15; %convert to C

plot(t,tasic,'linewidth',1.5)
box off
datetick('x','mmm'' yyyy')
axis tight| |

[Mohammad Abouali]

if you want to use rlon/rlat as Chad said here you need to change your destination point form 86.902,28.046 also to the rotated version, which seems you are doing.

But if you don't want to do that the not-rotated version of the coordinates are stored in lat/lon variable, and you can use the code I provided (interp2 doesn't work for non-rotated version).

It turns out to be faster to perform time-series interpolation using the ConstructPolyInterpolant2D, and non-rotated version. I designed that function particularly for this sort of application, that you have the same source and destination grid, but you want to interpolate a time series or you want to interpolate lot's of different fields. If you have just few interpolation that function might not be suitable but if you have 366 interpolation as it shows it works faster.

Answer 2

[Mohamed Elshamy]

Thanks for this information. Is there a way to reproject the whole field. I already found the lat and lon but they are matrices not vectors so interp2 did not accept them.

[Walter Roberson]

? interp2 does fine on 2D arrays for any of its arguments other than the interpolation method specification.

[Chad Greene]

A common problem with interp2 and geo data is the order. It's the issue of whether the grid was created by

[latgrid,longrid] = meshgrid(lat,lon);

or

[longrid,latgrid] = meshgrid(lon,lat);

Different datasets build their grids in different ways, so you might need to do

zi = interp2(latgrid,longrid,z,lati,loni);

or

zi = interp2(longrid,latgrid,z,loni,lati);

Answer 3

[Ashish Uthama]

I got this far.. then realized no exact match exists for that lat/lon. So, conceptually, what do you want to do? Pick a point closest to your required lat/lon? (closest in terms of euclidean distance? or is there another metric to measure the "distance" between two lats/lons?)

I used MAD below:

%%Inspect
ncFile = 'test.nc';
ncdisp(ncFile);

%%Read non-uniform (?) lat and lons
alllat = ncread(ncFile, 'lat');
alllon = ncread(ncFile, 'lon');

%%Find the index of the required lat/lon
lat = 28.046;
lon = 86.902;

%%Look for exact match
latInd = find(alllat==lat)
lonInd = find(alllon==lon)

%%Look for closest match (use Mean Absolute Difference)
madlatlon = abs(alllat-lat) + abs(alllon-lon);
[~, closestInd] = min(madlatlon(:)); % note - this is a linear index

closestLat = alllat(closestInd)
closestLon = alllon(closestInd)

%%Read that value:
[rlonInd, rlatInd] = ind2sub(size(alllat), closestInd); % convert to x,y index

% read all time for this lat/lon
tas = ncread(ncFile,'tas', [rlonInd, rlatInd, 1], [1 1 inf]);
tas = squeeze(tas); % make data single dimension
plot(tas); title(['tas at' num2str(closestLat) ' lat, ' num2str(closestLon) ' lon']);

Answer 4

[Chad Greene]

Use

t = datenum(1949,12,ncread('test.nc','time'));

to get time in Matlab's datenum format, which is days since year 0. A date vector contains year, month, and day, and is given by

dv = datevec(t);

The year would be

year = dv(:,1);

and

doy = datenum([dv(:,1:3), zeros(size(dv,1),3)]) - datenum([dv(:,1),...
    ones(size(dv,1),1), zeros(size(dv,1), 4)]);