MATLAB Answers

Extraction NetCDF time series to point

64 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
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
Phat Pumchawsaun
Phat Pumchawsaun on 8 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
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!