Data I am importing via NCREAD is not similar to values from third party websites

Below is the code I am using and attached is the image it should look like. Spatially, the data is similar but value-wise it is significantly different (unit conversions aside and interpolation methods). I cannot understand why the values are so different given they are pulling from the same source...

Could the data be changing within NCREAD? Or am I using NANSUM wrong? Any insight would be appreciated. Ryan

% ************************ GFS RAINFALL DATA ******************************
clear all

variable = [1];   %1 = rain, 2 = maxT, 3 = minT

variables = {'apcpsfc','tmax2m','tmin2m'}; 

latregions = [30 50];
lonregions = [250 290];  %originally 250 290    110W - 70W

run = '00';
date = datestr(now,'yyyy-mm-dd'); 

fullpath = 'http://nomads.ncep.noaa.gov:80/dods/gfs_0p50/gfs20141229/gfs_0p50_00z';

%---- Pull in variables : lat, lon... then rain within those lat/lons -----

'Pulling Lat/Lon/Time Info...'
lat  = ncread(fullpath,'lat',[1],[Inf]);
lon  = ncread(fullpath,'lon',[1],[Inf]);
time = ncread(fullpath,'time',[1],[Inf]);
'Pulling Lat/Lon/Time Info... DONE'

timesegs = numel(time);
increment = 8;              % Since model is every 3hrs  (24/3=8)
days = (timesegs-1)/increment;

lonST = find(lon == lonregions(1));
lonEND= find(lon == lonregions(2));
latST = find(lat == latregions(1));
latEND= find(lat == latregions(2));

   data = ncread(fullpath,'apcpsfc',[1 1 1],[Inf Inf Inf]);%,[1 1 1]);   %Pulls in rain data

regData = data(lonST:lonEND, latST:latEND, :);
regLon = lon(lonST:lonEND);
regLat = lat(latST:latEND);

% For precip take out first time segment b/c no acc. pcp in it 
if variable(1) == 1
 regData = regData(:,:,2:end);
 data    = data(:,:,2:end); 
end

%To put rain values in daily values
Dailyfulldata = zeros(numel(lon),numel(lat),days);
Dailyregdata  = zeros(numel(regLon),numel(regLat),days);

st=1;
endd = increment;

for i=1:days

     Dailyfulldata(:,:,i) = nansum(data(:,:,st:endd),3);
     Dailyregdata(:,:,i) = nansum(regData(:,:,st:endd),3); 

     st = endd + 1;
     endd = endd + increment;

end

% Adjust matrix and lat/lon coordinates for map and plotting

%adjDailyfulldata = zeros(numel(lat),numel(lon),days);
adjDailyregdata  = zeros(numel(regLat),numel(regLon),days);

adjLat = flipud(regLat);
adjLon = regLon.';

for i=1:days
    adjDailyregdata(:,:,i) = Dailyregdata(:,:,i).'; 
end

contour(regData(:,:,1).',30) %Compare this plot with attachment

[dpb]

Probably more useful to those here would be the actual data as read from the other source.

[Ryan]

Thanks dpb. I'll try saving the data in a .mat format.

[Ryan]

Hi dpb, because of the data limit, I am only able to attach for the variables:

regData (rawData attachment)

regData(:,:,1) (testData attachment)

[dpb]

Which of these are supposed to be commensurate? One is 3D the other 2D. If I try the first plane of regData as compared to testData I find

>> sum(testData(:)==0)
ans =
      5903
>> sum(sum(regData(:,:,1)==0))
ans =
     10295
>>

So, those are clearly incompatible datasets; do you think any of them really should be?

Answer 1

[Chad Greene]

Looks like the apcpsfc data are in kg/m^2, not inches.

[Ryan]

Hi Chad, thanks for your response.

From my understanding, kg/m^2 is millimeters, and the rain values I'm working with in my code are still higher than the image, metric conversions aside.

Answer 2

[Chad Greene]

I get the correct order of magnitude. Note that your image shows a log color scale and most to the data don't exceed 0.5 inches. Here I'm plotting the second day's daily precipitation on a linear color scale. I wonder if your problem is in the way you computed daily totals. Below, I re-summed using unique days.

fullpath = 'http://nomads.ncep.noaa.gov:80/dods/gfs_0p50/gfs20141229/gfs_0p50_00z';
data = ncread(fullpath,'apcpsfc');
lat  = ncread(fullpath,'lat');
lon  = ncread(fullpath,'lon');
time = ncread(fullpath,'time');

% Change dimensions of data to lat x lon x time: 
data = permute(data,[2 1 3]); 

% Get daily sums: 
uniqueDays = unique(floor(time)); 
DailyPrecip = NaN(length(lat),length(lon),length(uniqueDays)); 
for k = 1:length(uniqueDays)
    DailyPrecip(:,:,k) = nansum(data(:,:,time==uniqueDays(k)),3); 
end

figure
worldmap([21 52],[-135 -65]);
geoshow('usastatelo.shp', 'FaceColor', [.5 .5 1]);

pcolorm(lat,lon,mm2in(DailyPrecip(:,:,2)))
cb = colorbar('location','southoutside'); 
xlabel(cb,'mean daily precipitation in inches')
caxis([0 .5])
rgbmap('white','green','dark blue')

Above I used two of my own functions, mm2in and rgbmap, which are both available on the File Exchange.

[Ryan]

Thanks for this! However it seems odd as you take the 8th, 16th, 24th element. From my understanding each timestep is 3hr total rainfall, and I need 24hr totals for Day 1,2,... I've also noticed that the location of rainfall is slightly off compared to third party images. I'm realizing this problem is a lot deeper than I thought... Thanks your help above.

[Chad Greene]

Ah, oops. Above, the for loop should contain this:

DailyPrecip(:,:,k) = nansum(data(:,:,floor(time)==uniqueDays(k)),3);

The difference is the necessary floor command. In that I'm saying, for each unique day, sum all the data whose date given by floor(time) corresponds to that unique day. Without the floor command it was only grabbing the midnight data.

Answer 3

[Chad Greene]

Simplest solution: downsample_ts now supports sums and nansums. Here's the full code from data download to map creation:

% Read data: 
fullpath = 'http://nomads.ncep.noaa.gov:80/dods/gfs_0p50/gfs20141229/gfs_0p50_00z';
data = ncread(fullpath,'apcpsfc');
lat  = ncread(fullpath,'lat');
lon  = ncread(fullpath,'lon');
time = ncread(fullpath,'time');

% Change dimensions of data to lat x lon x time: 
data = permute(data,[2 1 3]); 

% Calculate daily totals: 
DailyPrecip = downsample_ts(data,time,'daily','nansum');

% Plot day 2: 
figure
worldmap([21 52],[-130 -65]);
pcolorm(lat,lon,mm2in(DailyPrecip(:,:,2)))
shading interp
rgbmap('white','lime green','dark green','dark blue','teal')
geoshow('usastatelo.shp', 'edgecolor',rgb('brown'),'facecolor','none');
cb = colorbar('location','southoutside'); 
xlabel(cb,'Total rainfall on Dec. 30, 2014 (inches)')
caxis([0 2])

[Ryan]

Thanks for this Chad it makes a lot more sense. Basically this is now where I'm at in my coding... Spatially I realized this is spot on, but notice the values themselves compared to third party vendors.. the values produced in both your program and mine is significantly higher than the image attached. This is my current endeavor, trying to figure out why this is..

[Chad Greene]

I disagree--the values in this plot seem to match your Rainfall24h.png plot exactly. Note that the third-party vendor contour image you provided uses a strange scale that starts as a log scale, then goes linear. Using clickz, I click in the middle of Colorado and get ~0.26 inches of rainfall. Near Rochester, NY has about 0.085 inches. Dallas has 0.005 inches and the very bottom tip of Texas has 0.04 inches. Hilton Head, SC has 0.395 inches. All of these values match the Rainfall24h.png image exactly.

[Ryan]

I'll double check this Chad and get back to you. Thanks

[Ryan]

Hey Chad, not trying to be argumentative whatsoever but how would you explain the apparent differences over the Atlantic? According to your most recent plot (and mine too) large areas of the Atlantic is well above 1"... verses the third-party vendor plot that has only small sections above 1". The area I'm looking at is between 70W and 65W.

[Chad Greene]

Ryan, using

contourfm(lat,lon,mm2in(DailyPrecip(:,:,2)),[1 1])
caxis([0 2])
colormap([1 1 1;0 0 0])

gives

Indeed, the patch over 1" is larger than the light blue patches in Rainfall24h.png. Keep in mind that both solutions are interpolating smoothish contours to sub-degree resolution from one-degree resolution source data. Are your other vendors using sub-one-degree source data, or simply a different interpolation scheme? Other possible differences might be timing. Is the midnight data applied to the previous day or the next day? That is, does the 3 hours of data timestampped midnight get applied to 9:00 pm to midnight? Or midnight to 3:00 am? Or 10:30 pm to 1:30 pm? Is it a difference of binning times based on GMT versus Eastern Daylight time?

[Ryan]

Hey Chad,

I think I found the solution... The data we are pulling in is 3hr time segments.. so we have data at hour 3,6,9,12,... As I was searching around in the data I discovered that hr 3 is accumulated precipitation for the past 3 hours and that hr 6 is accumulated precip for the past SIX hours. Hr 9 is past 3hrs, hr 12 is past 6 hrs... Why does it seem to be like this?? I have no idea.. but basically rainfall values are now very similar to third party solutions.. Very very strange! This is something I just uncover about 20 minutes ago.. will keep you updated! Thanks again, Chad

[Ryan]

Hey Chad, just an update on this, what I mentioned above is correct as the rain volumes are very similar to third party images... however spatially something seems off. I will post a picture of the differences soon..

[Ryan]

 Oddly, my program matches the third party images when taking in data from a different source called "GFS HD" verses "GFS 0.5 degree". Both are the same resolution and seem to be the same... yet one works and the other does not. All very very strange..

URL we've been using: http://nomads.ncep.noaa.gov:80/dods/gfs_0p50/gfs20150113/gfs_0p50_00z

New URL I tried: http://nomads.ncep.noaa.gov:80/dods/gfs_hd/gfs_hd20150113/gfs_hd_00z

[Ryan]

Oh and another thought, this as least gives me confidence that my code isn't completely out to lunch... but still leaves me scratching my head as far as why the 0.5 degree or 0.25 degree datasets are different