Extraction NetCDF time series to point

30 views (last 30 days)
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
  2 Comments
Meg Noah
Meg Noah on 8 Jan 2020
Please attach the 'precip.mon.mean.nc' file or give a link to where it may be downloaded on the internet.
Phat Pumchawsaun
Phat Pumchawsaun on 8 Jan 2020
Hi,
The link for data is as above.
Thanks

Sign in to comment.

Accepted Answer

Meg Noah
Meg Noah on 8 Jan 2020
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']);
RainRate.png
  4 Comments
Meg Noah
Meg Noah on 8 Jan 2020
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
Phat Pumchawsaun on 8 Jan 2020
Edited: Phat Pumchawsaun on 9 Jan 2020
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
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

Sign in to comment.

More Answers (2)

Meg Noah
Meg Noah on 9 Jan 2020
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:
GlobalMean.png
  2 Comments
Phat Pumchawsaun
Phat Pumchawsaun on 9 Jan 2020
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
Meg Noah on 9 Jan 2020
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));

Sign in to comment.


Mir Talas Mahammad Diganta
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?

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!