Data I am importing via NCREAD is not similar to values from third party websites
조회 수: 3 (최근 30일)
이전 댓글 표시
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
댓글 수: 4
dpb
2014년 12월 29일
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?
답변 (3개)
Chad Greene
2014년 12월 29일
편집: Chad Greene
2014년 12월 29일
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/172692/image.png)
Above I used two of my own functions, mm2in and rgbmap, which are both available on the File Exchange.
댓글 수: 2
Chad Greene
2014년 12월 30일
Ah, oops. Above, the for loop should contain this:
DailyPrecip(:,:,k) = nansum(data(:,:,floor(time)==uniqueDays(k)),3);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/180917/image.png)
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.
Chad Greene
2014년 12월 30일
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])
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/172717/image.png)
댓글 수: 9
참고 항목
카테고리
Help Center 및 File Exchange에서 Modify Image Colors에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!