Extraction NetCDF time series to point

18 次查看（过去 30 天）
Phat Pumchawsaun
Phat Pumchawsaun 2020-1-7
回答： Mir Talas Mahammad Diganta 2021-7-29
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 个评论
Meg Noah
Meg Noah 2020-1-8
Please attach the 'precip.mon.mean.nc' file or give a link to where it may be downloaded on the internet.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Meg Noah
Meg Noah 2020-1-8
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 个评论
显示 2更早的评论
Meg Noah
Meg Noah 2020-1-8
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 2020-1-8
编辑：Phat Pumchawsaun 2020-1-9
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

请先登录，再进行评论。

更多回答（2 个）

Meg Noah
Meg Noah 2020-1-9
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 个评论
Phat Pumchawsaun
Phat Pumchawsaun 2020-1-9
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 2020-1-9
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));

请先登录，再进行评论。

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?

类别

在 Help Center 和 File Exchange 中查找有关 NetCDF 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by