0

我正在处理 MERRA 数据Matlab。我有兴趣绘制MERRA300.prod.assim.tavg1_2d_rad_Nx 产品的指定坐标的数据。

关于如何提取特定位置数据（基于纬度/经度），我需要您的帮助吗？我附上了我处理过的代码。我也在这里和这里看到了例子。

file_name = 'MERRA300.prod.assim.tavg1_2d_rad_Nx.20150101.hdf';
GRID_NAME = 'EOSGRID'
file_id = hdfgd('open', file_name, 'rdonly');
grid_name='EOSGRID';
grid_id = hdfgd('attach', file_id, grid_name);
% datafield_name='U500';
datafield_name='ALBEDO';

[data, fail] = hdfgd('readfield', grid_id, datafield_name, [], [], []);
data=squeeze(double(data(:,:,1)));
data=data';

% size(data)

datafield_name='XDim';
[lon, status] = hdfgd('readfield', grid_id, datafield_name, [], [], []);
lon=double(lon);
datafield_name='YDim';
[lat, status] = hdfgd('readfield', grid_id, datafield_name, [], [], []);
lat=double(lat);
hdfgd('detach', grid_id);
hdfgd('close', file_id);
SD_id = hdfsd('start',file_name, 'rdonly');

% datafield_name='U500';
datafield_name='ALBEDO';

% size(lat)
% size(lon)
sds_index = hdfsd('nametoindex', SD_id, datafield_name);

sds_id = hdfsd('select',SD_id, sds_index);
fillvalue_index = hdfsd('findattr', sds_id, '_FillValue');
[fillvalue, status] = hdfsd('readattr',sds_id, fillvalue_index);

missingvalue_index = hdfsd('findattr', sds_id, 'missing_value');
[missingvalue, status] = hdfsd('readattr',sds_id, missingvalue_index);

long_name_index = hdfsd('findattr', sds_id, 'long_name');
[long_name, status] = hdfsd('readattr',sds_id, long_name_index);

units_index = hdfsd('findattr', sds_id, 'units');
[units, status] = hdfsd('readattr',sds_id, units_index);

scale_index = hdfsd('findattr', sds_id, 'scale_factor');
[scale, status] = hdfsd('readattr',sds_id, scale_index);

scale = double(scale);
offset_index = hdfsd('findattr', sds_id, 'add_offset');
[offset, status] = hdfsd('readattr',sds_id, offset_index);

offset = double(offset);
hdfsd('endaccess', sds_id);
data(data==fillvalue) = NaN;
data(data==missingvalue) = NaN;
data = data*scale + offset ;
latlim = [floor(min(min(lat))),ceil(max(max(lat)))];
lonlim = [floor(min(min(lon))),ceil(max(max(lon)))];
min_data = floor(min(min(data)));
max_data = ceil(max(max(data)));

  % Find indexes for region along longitude.
 lon_region = (lon > 18.0 & lon < 45.0);
 i = find(lon_region,1,'first');
 j = find(lon_region,1,'last');

 % Subset data using the above indices.
 % lat = lat(i:j);
 % lon = lon(i:j);
 % time = time(i:j);
 % data = data(:,i:j);

 % Find indexes region along latitude.
 lat_region = (lat > 18.0 & lat < 45.0);
 i = find(lat_region,1,'first');
 j = find(lat_region,1,'last');

 % Subset data using the above indices.
 lat = lat(i:j);
 lon = lon(i:j);
 % time = time(i:j);
 % data = data(:,i:j);


f = figure('Name', file_name,'visible','on');

axesm('MapProjection','eqdcylin','Frame','on','Grid','on', ...
      'MapLatLimit',latlim,'MapLonLimit',lonlim, ...
      'MeridianLabel','on','ParallelLabel','on')
coast = load('coast.mat');

% surfacem(lat,lon,data);
surfm(lat,lon,data);
colormap('Jet');
caxis([min_data max_data]); 

% Change the value if you want to have more than 10 tick marks.
ntickmarks = 10;
granule = floor((max_data - min_data) / ntickmarks);

h=colorbar('YTick', min_data:granule:max_data);

plotm(coast.lat,coast.long,'k')

title({file_name; ...
       [long_name  [' at TIME=0']]}, ...
       'Interpreter', 'None', 'FontSize', 16,'FontWeight','bold');

set (get(h, 'title'), 'string', units, 'FontSize', 16,'FontWeight','bold');

% Use the following if your screen isn't too big (> 1024 x 768). 
scrsz = get(0,'ScreenSize');

% The following fixed-size screen size will look better in PNG if
% your screen is too large.
scrsz = [1 1 800 600];
set(f,'position',scrsz,'PaperPositionMode','auto');
4

0 回答 0