我有一个 modis hdf 文件，mcd19a2 产品。这些文件是使用正弦投影存储的。然后将图块划分为 1200 x 1200 网格（或者在 1200 x 1200 矩阵中，你可以说）。经度和纬度数据也被存储，但不是在 WG84 纬度/经度值中。它们以米为单位。我想从特定文件中提取 lat 和 on 信息，然后将其转换为 wg84 lat/lon。我正在使用 MATLAB 版本 2021。
这是一个函数，我使用的是“findindex”。问题是瓷砖是在正弦投影中进行的，即它具有不同的左上角和左下角经度值，而从代码中，我得到了相似的值。
我正在使用的文件 MCD19A2.A2021121.h24v06.006.2021124142622.hdf https://drive.google.com/file/d/1g1hnNFIeSef8C17L_6Mq8g0RaJxe1mW1/view?usp=sharing
关于上下纬度/经度的纬度/经度信息如下：
iv ih 0: ll lon ll lat 1: ul lon ul lat 2: ur lon ur lat 3: lr lon lr lat
6 24 63.6095 19.9322 69.0253 30.0331 80.8440 29.9993 74.4929 19.9002
我用来转换的代码；
%Delhi coordinates
% LatLon=[28.8855033,76.837768;28.416234,77.329407]
filename='./MCD19A2.A2021121.h24v06.006.2021124142622.hdf';
proj=1;
% if nargin < 3
% disp(['Error! three mandatory input arguments are required for this '...
% 'function, please see Function Description for more information']);
% return
% end
projType={'Sinusoidal'; 'Global Cylindrical Equal-Area'};
if proj>2
coordX_idx=0; coordY_idx=0; projType='Wrong projType';
msgbox(sprintf([' "proj = %d" does not match any acceptable projection '...
'types, only 1: Sinusoidal, and 2: Cylindrical Equal-Area Scalable '...
'Earth Grid (EASE) are accepted by this function; '...
' please try again'],proj),...
'Invalid projection type requested!','Error','modal');
return %abort code if mismatch detected
end
% Specify Earth's Equatorial Radius based on the Projection Type:
if proj==1
R=6371007.181; %Radius of Earth in meters (Sinusoidal)
%False Easting: 0.0
%False Northing: 0.0
%ISin NZone: 86400 (see
%http://geospatialmethods.org/documents/ppgc/ppgc.html)
elseif proj ==2
R = 6371228; %m R of Earth (EASE & all spherical projns, no eccent.)
%All spherical projns, eccentricity should not be specified, and
%the default equatorial radius is 6371228 m.
else %more projections to work out:
R=6378137.0; %R Earth (m), WGS84 ellipsoid (eccentricity: 0.081819190843)
%All elliptical projtns: eccent. 0.082271673 & R 6378206.4 (defaults).
end
% Circumference of a circle: %C=pi*d so C=pi*R*2
C=pi*(R*2); %circumference of Earth in Equator (meters)
dist1degree=C/360;%metric distance equivalent to 1 degree in Equator
lambda0=0;
% read HDF EOS file information from the given file:
% http://www.mathworks.com/company/newsletters/digest/nov02/earth_pt3.html
% left-right-top-botton bounds and rows-columns of the HDF file are needed
geoInfo = hdfinfo(filename,'eos'); % HDF file info in EOS mode
xDim=geoInfo.Grid(1).Columns; % columns of data
yDim=geoInfo.Grid(1).Rows; % rows of data
ULx0=geoInfo.Grid(1).UpperLeft(1); % UpperLeft bound X
LRx0=geoInfo.Grid(1).LowerRight(1); % LowerRight bound X
ULy0=geoInfo.Grid(1).UpperLeft(2); % UpperLeft bound Y
LRy0=geoInfo.Grid(1).LowerRight(2); % LowerRight bound Y
% add offset:
xOffset =0;%offset (false Easting)
yOffset =0;%offset (false Northing)
% reformat UL/LR corners
LRx=abs(LRx0 - ULx0); %let LRx as the largest x in dataset
LRy=abs(LRy0 - ULy0); %let LRy as the largest y in dataset
ULx=ULx0 - ULx0; %turn ULx as origin and let it be 0
ULy=ULy0 - ULy0; %turn ULy as origin ane let it be 0
%LLy=LRy; %LLy is the same as LRy (not needed)
%LLx=ULx; %LLx is the same as ULx (not needed)
%ULx=ULx0;ULy=ULy0;LRx=LRx0;LRy=LRy0; %no change (diagnostic)
% calculate resolution across X and Y:
xRes=(LRx-ULx)/xDim; %resolution along x (longitude) in meters(positive)
%yRes=(ULy-LLy)/yDim; %resolution along y (latitude) in meters(negative)
yRes=(LRy-ULy)/yDim; %resolution along y (latitude) in meters(positive)
% determine where to start searching points x-y; since the coordinates
% start from UL corner of pixel, but we want center of pixel, xRes/2 or
% yRes/2 will be added to the start2look point:
xStartLooking=ULx+xRes/2; %start looking from UL corner
yStartLooking=ULy+yRes/2; %start looking from UL corner
n=size(LatLon,1); %number of points to look for
coordX_idx=zeros(n,1); %point's x index in HDF file: initialize to 0
coordY_idx=zeros(n,1); %point's y index in HDF file: initialize to 0
mdist2point=ones(n,1)*(abs(xRes)+abs(yRes))*10;
% start looking for all points based on Lat-Lon coordinates:
for i2=1:n
phi=LatLon(i2,1); %latitude coordinate
lambda=LatLon(i2,2); %longitude coordinate
if proj == 1 % Sinusoidal Projection:
x=(lambda-lambda0)*cos(phi*pi/180);%lon in Sin proj in local lat
y=phi; %latitude in Sin projection
x=abs(x*dist1degree-ULx0); %valid for north america as in discussion (conclded above equator)
y=abs(y*dist1degree-ULy0); %
% x=abs(x)*dist1degree-abs(ULx0); %x in m in Global Sin proj (for
% newzealand, below equator)
% y=abs(y)*dist1degree-abs(ULy0); %y in m in Global Sin proj
xMax=180*cos(phi*pi/180); %max longitude in Sin proj in local lat
%cLocal=C*(xMax*2)/360; %circumference of Earth in local latitude (m)
%dist1degLocal = cLocal/360;%metric dist = 1 deg in local lat
elseif proj ==2 % EASE projection:
% Formula (http://nsidc.org/data/ease/ease_grid.html):
% �r = r0 + R/C * lambda * cos(30)
% �s = s0 - R/C * sin(phi) / cos(30)
% �h = cos(phi) / cos(30)
% �k = cos(30) / cos(phi)
% r Column coordinate
% s Row coordinate
% h Particular scale along meridians
% k Particular scale along parallels
% lambda Longitude in radians (lambda*pi/180)
% phi Latitude in radians (phi*pi/180)
% R Radius of the Earth = 6371.228 km in Equator
% C Nominal cell size
% r0 Map origin column
% s0 Map origin row
%Note: true latitude (30) must also be converted to radians
r0=floor(xDim/2);
s0=floor(yDim/2);
r_x = r0 + R/xRes * (lambda*pi/180) * cos(30*pi/180); %x coord in EASE proj
s_y = s0 - R/yRes * sin(phi*pi/180) / cos(30*pi/180); %y coord in EASE proj
%h=cos(phi*pi/180)/cos(30*pi/180);%scale along meridians(not needed)
%k=cos(30*pi/180)/cos(phi*pi/180);%scale along parallels(not needed)
x=r_x*xRes; % x in Metric scale
y=s_y*yRes; % y in Metric scale
end
% Now search for the nearest pixel to the given coordinates:
xMetric=x+xOffset; %x coordinate in meters adjusted for add-offset
yMetric=y+yOffset; %y coordinate in meters adjusted for add-offset
distInitial=(abs(xRes)+abs(yRes))*10; %approx initial dist based on grid res
pixelX=xStartLooking; %initial pixel x to check the distance
for i=1:xDim
pixelY=yStartLooking; %initial pixel y to check the distance
for j=1:yDim
sq1=(pixelY-yMetric)^2; % lat part of dist func
sq2=(pixelX-xMetric)^2; % lon part of dist func
dFound=sqrt(sq1+sq2); % distance function
sq1fo(i,j)=sq1;
sq1foww(i,j)=pixelY;
sq2fo(i,j)=sq2;
sq2foww(i,j)=pixelX;
dfo(i,j)=dFound; %store distances (diagnostic)
if dFound<distInitial
coordX_idx(i2)=i; % point's x index in HDF file (returned)
coordY_idx(i2)=j; % point's y index in HDF file (returned)
distInitial=dFound; % replace initial dist with found dist
end
pixelY=pixelY+yRes; % move 1 pixel forward along Lat
end
pixelX=pixelX+xRes; % move 1 pixel forward along Lon
end
mdist2point(i2)=distInitial; %distance from given point to its closest
%pixel (indexes of which are given in coordX_idx & coordY_idx)
end
% lat_deg= sq1foww+ULx0 ;
% lat_deg1=(sq1foww./dist1degree)+ULx0;
% check_lat=(1200*xRes)+ULx0;
converting back to lat and lon
for i=coordY_idx
ws_lat=abs((sq2foww-ULy0)./dist1degree);%Lattitude conversion of metres sin proj
% to wsg84 datum coordinate
% system
end
%
for i=1:xDim
ws_lon=abs((sq1foww1+ULx0)./dist1degree); %longitude conversion from sin proj
factor_lat=cos(ws_lat*pi/180); % to wsg84 datum cooridnate system
factor_lat1=factor_lat.';
ws_lon1=ws_lon./factor_lat1; % final conversion of lon
end
请建议。如果有任何其他简单的方法可以做到这一点。