Extraction NetCDF time series to point

Hi,

I have NetCDF file (file name is precip.mon.mean.nc at this link ftp://ftp.cdc.noaa.gov/Datasets/cmap/std/precip.mon.mean.nc) with four variables as column as;Lat; Lon;Precip;Time. I want to extract monthly grid time series data available from 1979 to 2019 to my station point. My code is as below:

file = ('precip.mon.mean.nc'); %openfile
time = ncread(file,'time');
prep = ncread(file,'precip'); 
lon = ncread(file,'lon');
lat = ncread(file,'lat');
latlim = [19.787500 18.029167]; %N to S
lonlim = [101.904167 103.391667]; %W to E
from_date = '01-01-1979'; to_date = '01-11-2019';
time_datenum = datenum('01-01-1979','dd-mm-yyyy');
date_match = time_datenum >= datenum(from_date) && time_datenum <= datenum(to_date);
selected_prep = prep(:,:,date_match);

Yet, I couldn't get data in select_prep which is the new one I want to store data. Could you guys help me on this.

Thanks

Phat

[Meg Noah]

Please attach the 'precip.mon.mean.nc' file or give a link to where it may be downloaded on the internet.

[Phat Pumchawsaun]

Hi,

The link for data is as above.

Thanks

Answer 1

[Meg Noah]

Here's a solution that is workable:

clc
close all
clear all

% ftp://ftp.cdc.noaa.gov/Datasets/cmap/std/precip.mon.mean.nc
filename = ('precip.mon.mean.nc');
info = ncinfo(filename);
unitsTime = info.Variables(3).Attributes(1);
unitsPrep = info.Variables(4).Attributes(3).Value;
validRangePrep = info.Variables(4).Attributes(2).Value;
labelPrep = info.Variables(4).Attributes(1).Value;

time = ncread(filename,'time');
tDatenum = datenum(1800,1,1,time,0,0);

prep = ncread(filename,'precip'); 
lon = double(ncread(filename,'lon'));
lat = double(ncread(filename,'lat'));
latlim = [18.029167 19.787500 ];  % N to S
lonlim = [101.904167 103.391667]; % W to E

% note there are no longitudes in your limits
% choices are idxLon = 41 lon = 101.2500
%          or idxLon = 42 lon = 103.7500
idxLon = 42;
idxLat = find(latlim(1) <= lat & lat <= latlim(2));

from_datenum = datenum(1979,1,1,0,0,0);
to_datenum = datenum(2019,1,11,0,0,0);
idxDatenum = find(from_datenum <= tDatenum & tDatenum <= to_datenum);

myPrepData = squeeze(prep(idxLon,idxLat,idxDatenum));
myPrepDatenum = tDatenum(idxDatenum);

figure('color','white','Position',[70 500 1200 450]);
hold on
box on
xlim([datenum(1978,1,1,0,0,0) datenum(2021,1,1,0,0,0)]);
set(gca,'xtick',datenum(1979:2:2020,1,1,0,0,0));

plot(myPrepDatenum,myPrepData);
datetick(gca,'x','yyyy','keepticks');
ylabel([labelPrep ' [' unitsPrep ']']);
ylim(validRangePrep);
title(['Latitude = ' num2str(lat(idxLat)) ' \circN Longitude = ' ...
    num2str(lon(idxLon)) ' \circE']);

[Meg Noah]

Note:

unitsTime = info.Variables(3).Attributes(1);

Explains how to interpret and use the time variable. It's shown above.

[Phat Pumchawsaun]

Hi Meg,

Thanks so much for your suggestion. One more thing, how should I type index for lat lon if my position is single point. I've tried as below but it's doesn't work.

idxLon = 102.5;
idxLat = 18.5;

[Meg Noah]

Hi Phat Pumchawsaun,

The idx values are integers. They are not the latitude or longitude value, but the index in to the array of latitude, longitude values of the image. I'm hoping that I've correctly indexed into the array 'prep'. What I should do is create a 2D world map at one timestep to verify that it's not getting the values from 90-lat, since sometimes these get flipped north-south.

so...

LatitudeOfPoint = lat(idxLat);

LongitudeOfPoint = lon(idxLon);

There is no point that is within your range of longitudes. There are two points fairly close to that range. The descretization of the longitudes in the prep array is bigger than your range of site location values.

-Meg

[Phat Pumchawsaun]

Hi Meg,

Thanks so much again. I've modified a bit for your original code to extarct daily data but it doesn't work. The reasons are; (1) it seems like time matrix is minus; (2) it seems like time matrix exceeds array bounds (must not exceed 365); and (3) empty matrix after extraction

Link for data is ftp://ftp.cdc.noaa.gov/Datasets/cpc_global_precip/precip.2019.nc

Code is as folliwings;

clc
close all
clear all
filename = ('precip.2019.nc');
info = ncinfo(filename);
unitsTime = info.Variables(3).Attributes(2);
unitsPrep = info.Variables(4).Attributes(2).Value;
validRangePrep = info.Variables(4).Attributes(9).Value;
labelPrep = info.Variables(4).Attributes(7).Value;
time = ncread(filename,'time');
prep = ncread(filename,'precip'); 
lon = double(ncread(filename,'lon'));
lat = double(ncread(filename,'lat'));
% 
location_lat = knnsearch(lat,18.5);
location_lon = knnsearch(lon,102.5);
tDatenum = datenum('01-01-2019 00:00:00','dd-mm-yyyy HH:MM:SS');
from_date = '01-01-2019 00:00:00'; to_date = '31-12-2019 00:00:00';
idxDatenum = from_date <= tDatenum & tDatenum <= to_datenum;
myPrepData = squeeze(prep(location_lon,location_lat,idxDatenum));
myPrepDatenum = tDatenum(time);
figure('color','white','Position',[70 500 1200 450]);
hold on
box on
xlim([datenum(2019,1,1,0,0,0) datenum(2019,12,31,0,0,0)]);
set(gca,'xtick',datenum(1979:2:2020,1,1,0,0,0));
plot(myPrepDatenum,myPrepData);
datetick(gca,'x','yyyy','keepticks');
ylabel([labelPrep ' [' unitsPrep ']']);
ylim(validRangePrep);
title(['Latitude = ' num2str(lat(location_lat)) ' \circN Longitude = ' ...
    num2str(lon(location_lon)) ' \circE']);

Thanks in advance.

Phat

Answer 2

[Meg Noah]

Hi.

This is a completely different file format. The data are on a finer spatial grid, but there are only 365 'times'. The times are days of the year where the data are the mean values for that day of year, averaged over that day of year over many years. Here's how to read that data. I don't have the stats package. There are 4 grid cells equally distant from your site.

clc
close all
clear all
filename = ('precip.day.1981-2010.ltm.nc');
info = ncinfo(filename);

unitsPrep = info.Variables(5).Attributes(12).Value;
% validRangePrep = info.Variables(5).Attributes(8).Value;
labelPrep = info.Variables(5).Attributes(11).Value;
validRangePrep = info.Variables(5).Attributes(10).Value;

% data are a time series of days averaged from '1981/01/01 - 2010/12/31'
% for each day of the year - we don't care about the year
% time = ncread(filename,'time');
% unitsTime = info.Variables(3).Attributes(2).Value;
% fprintf(1,'Time Units: %s\n',unitsTime);

% time12 = time*(-1);
prep = ncread(filename,'precip'); 
lon = double(ncread(filename,'lon'));
lat = double(ncread(filename,'lat'));
% 
% location_lat = knnsearch(lat,18.5);
% location_lon = knnsearch(lon,102.5);
dist = abs(lat - 18.5);
idxLat = find(dist == min(dist));
dist = abs(lon - 102.5);
idxLon = find(dist == min(dist));
% print out locations found
nsites = length(idxLat)*length(idxLon);
fprintf(1,'Found %d locations near site\n',nsites);
for kLat = 1:length(idxLat)
    for kLon = 1:length(idxLon)
         fprintf(1,'Location: Lat = %f Lon = %f\n', lat(idxLat(kLat)), lon(idxLon(kLon)));
    end
end

% hardwiring getting all 4 sites
DOY = 1:365; DOY = DOY';
GlobalLTM = table(DOY);
GlobalLTM.Site1 = squeeze(prep(idxLon(1),idxLat(1),:));
GlobalLTM.Site2 = squeeze(prep(idxLon(1),idxLat(2),:));
GlobalLTM.Site3 = squeeze(prep(idxLon(2),idxLat(1),:));
GlobalLTM.Site4 = squeeze(prep(idxLon(2),idxLat(2),:));
GlobalLTM.Properties.VariableUnits = {'Day Of Year','mm','mm','mm','mm'};
GlobalLTM.Properties.UserData.Site1 = [lat(idxLat(1)), lon(idxLon(1))];
GlobalLTM.Properties.UserData.Site2 = [lat(idxLat(2)), lon(idxLon(1))];
GlobalLTM.Properties.UserData.Site3 = [lat(idxLat(1)), lon(idxLon(2))];
GlobalLTM.Properties.UserData.Site4 = [lat(idxLat(2)), lon(idxLon(2))];
GlobalLTM.Properties.UserData.Description = {info.Variables(5).Attributes.Value}';

strSite1 = ['Site1 = ' ...
    num2str(GlobalLTM.Properties.UserData.Site1(1)) ...
    ' \circN, ' ...
     num2str(GlobalLTM.Properties.UserData.Site1(2)) ' \circE'];
strSite2 = ['Site2 = ' ...
    num2str(GlobalLTM.Properties.UserData.Site2(1)) ...
    ' \circN, ' ...
     num2str(GlobalLTM.Properties.UserData.Site2(2)) ' \circE'];
strSite3 = ['Site3 = ' ...
    num2str(GlobalLTM.Properties.UserData.Site3(1)) ...
    ' \circN, ' ...
     num2str(GlobalLTM.Properties.UserData.Site3(2)) ' \circE'];
strSite4 = ['Site4 = ' ...
    num2str(GlobalLTM.Properties.UserData.Site4(1)) ...
    ' \circN, ' ...
     num2str(GlobalLTM.Properties.UserData.Site4(2)) ' \circE'];
 
figure('color','white','Position',[70 500 1200 450]);
set(gca,'fontname','Ariel','fontsize',14,'fontweight','bold');
hold on
box on
plot(GlobalLTM.DOY,GlobalLTM.Site1,'DisplayName',strSite1);
plot(GlobalLTM.DOY,GlobalLTM.Site2,'DisplayName',strSite2);
plot(GlobalLTM.DOY,GlobalLTM.Site3,'DisplayName',strSite3);
plot(GlobalLTM.DOY,GlobalLTM.Site4,'DisplayName',strSite4);
legend('location','best');
ylabel([labelPrep ' [' unitsPrep ']']);
title({'CPC Global Precipitation: Daily Long Term Mean total of precipitation'; ...
    'time: mean within days time: mean over days time: mean over years'});
xlim([0 365]);
% tickmarks at first day of each month
ttt = datenum(2020,1:12,1,0,0,0);
ttt = ttt - ttt(1) + 1;
set(gca,'xtick',ttt);
set(gca,'xticklabel',{'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug', ...
    'Sep','Oct','Nov','Dec'});

% save it - make a spreadsheet
save('GlobalLTM.mat','GlobalLTM');
writetable(GlobalLTM,'GlobalLTM.xlsx');

The script produces this plot:

[Phat Pumchawsaun]

Hi Meg,

Bug thanks to suggestion. One more thing, regarding to previouse code, why did you suggest to use index for latitude as 42 and longtitude as the code as followings;

latlim = [18.029167 19.787500 ];  % N to S
lonlim = [101.904167 103.391667]; % W to E
% note there are no longitudes in your limits
% choices are idxLon = 41 lon = 101.2500
%          or idxLon = 42 lon = 103.7500
idxLon = 42;
idxLat = find(latlim(1) <= lat & lat <= latlim(2));

Thanks

Phat

[Meg Noah]

It is a type-o. There are no longitudes in your limits for the annual daily averages. It's a smaller grid. Set the index for the longiutde array to 42. But you can search the latitudes for a nearby one.

This strategy also would work:

dist = abs(lon - 102.5);
idxLon = find(dist == min(dist));

Answer 3

[Mir Talas Mahammad Diganta]

@Meg Noah

hi Meg. I have a similar problem and expecting your help.

I want to plot two time series (please see the attched image) for the entire area with netcdf file. I have six netcdf files corresponding each year (file_link: https://drive.google.com/drive/folders/1iGMwaYh-rqUmZ8uC6mgfKk6gInnsubqN?usp=sharing).

I have also merged these files into one file (file_link: https://drive.google.com/file/d/1UDS4dar0n_K39bk9e0xNlMskOWnG0f88/view?usp=sharing).

Can you help me with the matlab script?