RE: Extraction of NetCDF to certain point and location

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

Answer 1

[Meg Noah]

Here's a solution:

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']);

It produces this plot: