我在 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;