0

我在 HDF 文件中有一些 MODIS 数据(MOD17A2)。数据都显示一个相同的网格(1200*1200 像素)一年。我只需要为每个 HDF 文件检索一个像素数据(我的目标点),最后为我的目标点创建一年的时间序列数据。现在我可以在 matlab 中查看 1200*1200 矩阵数据,并知道 4 个网格角和网格中心的纬度/经度。如何在 1200*1200 矩阵中定位我的目标点(已知纬度/经度)?

clear all; clc;
%Opening the HDF-EOS2 Grid File
FILE_NAME='MOD17A2.A2012273.h09v04.005.2013230033352.hdf';
file_id = hdfgd('open', FILE_NAME, 'rdonly');

%Reading Data from a Data Field
GRID_NAME='MOD_Grid_MOD17A2';
grid_id = hdfgd('attach', file_id, GRID_NAME);

DATAFIELD_NAME='Gpp_1km';

[data1, fail] = hdfgd('readfield', grid_id, DATAFIELD_NAME, [], [], []);

%Convert M-D data to 2-D data
data=data1;

%Convert the data to double type for plot
data=double(data);


% This file contains coordinate variables that will not properly plot. 
% To properly display the data, the latitude/longitude must be remapped.

[xdimsize, ydimsize, upleft, lowright, status] = hdfgd('gridinfo', grid_id);
%Reading attributes from the data field
SD_id = hdfsd('start',FILE_NAME, 'rdonly');
DATAFIELD_NAME='Gpp_1km';

sds_index = hdfsd('nametoindex', SD_id, DATAFIELD_NAME);

sds_id = hdfsd('select',SD_id, sds_index);

%Reading filledValue from the data field
fillvalue_index = hdfsd('findattr', sds_id, '_FillValue');
[fillvalue, status] = hdfsd('readattr',sds_id, fillvalue_index);
%The _FillValue in the file contains value 32767 but the actual value
%observed from the data is 32766.
fillvalue = 32761:32767;

%Reading valid_range from the data field
valid_range_index = hdfsd('findattr', sds_id, 'valid_range');
[valid_range, status] = hdfsd('readattr',sds_id, valid_range_index);

%Reading units from the data field
units_index = hdfsd('findattr', sds_id, 'units');
[units, status] = hdfsd('readattr',sds_id, units_index);

%Reading scale_factor from the data field
scale_index = hdfsd('findattr', sds_id, 'scale_factor');
[scale, status] = hdfsd('readattr',sds_id, scale_index);

%Reading add_offset from the data field
offset_index = hdfsd('findattr', sds_id, 'add_offset');
[offset, status] = hdfsd('readattr',sds_id, offset_index);

% Reading long_name from the data field
long_name_index = hdfsd('findattr', sds_id, 'long_name');
[long_name, status] = hdfsd('readattr',sds_id, long_name_index);

%Terminate access to the corresponding data set
hdfsd('endaccess', sds_id);
%Closing the File
hdfsd('end', SD_id);

%Replacing the filled value with NaN
for i=1:length(fillvalue)
data(data==fillvalue(i)) = NaN;
end
%Replacing the out of range values with NaN
data(data < valid_range(1)) = NaN;
data(data > valid_range(2)) = NaN;

%Multiplying scale and adding offset, the equation is scale *(data-offset).
data = (data-offset) * scale*1000;
4

1 回答 1

1

您可以使用 HDF-EOS2 转储器工具先在 ASCII 文本文件中生成纬度/经度,然后在 MATLAB 中读回它们。利用此技术的完整代码可在

http://hdfeos.org/zoo/LPDAAC/MOD17A2_PsnNet_1km.m

上述代码使用以下 HDF-EOS2 转储工具选项以 ASCII 格式从 MOD17A2 转储纬度和经度。

$eos2dump -c1 MOD17A2.A2007113.h11v09.005.2007136163924.hdf > lat_MOD17A2.A2007113.h11v09.005.2007136163924.output
$eos2dump -c2 MOD17A2.A2007113.h11v09.005.2007136163924.hdf > lon_MOD17A2.A2007113.h11v09.005.2007136163924.output

然后,您可以从下图中检查您感兴趣的位置,或从通过读取 ASCII 文件创建的 MATLAB 数组中搜索纬度/经度值。

上述 MATLAB 代码的绘图

于 2013-12-19T04:54:58.780 回答