如何用Python GDAL读取MODIS HDF文件的位置坐标？

原因分析 #

1. 库的定位差异：pyhdf是专门针对HDF4格式的科学数据集读取库，会直接枚举文件内所有数据集（包括Longitude/Latitude这类辅助坐标数据集）；而GDAL核心是处理栅格地理数据，默认优先识别带有标准地理参考的栅格子数据集，对于MODIS HDF中与主变量同维度的经纬度数组，不会自动在GetSubDatasets()结果中展示。
2. 代码变量错误：你的代码中file.GetSubDatasets()是错误的，应该使用gdal_ds.GetSubDatasets()，不过即使修正后，GDAL默认返回的子数据集也可能不包含经纬度，因为它们不属于GDAL默认识别的"栅格波段"范畴。

用GDAL提取经纬度坐标的方法 #

方法1：直接构造子数据集路径打开 #

MODIS HDF4的子数据集路径遵循固定格式：HDF4_EOS:EOS_SWATH:"文件名.hdf":数据集组名:数据集名。对于MYD04_3K文件，经纬度的数据集组是mod04，所以可以直接构造路径读取：

from osgeo import gdal

filename = 'MYD04_3K.A2016001.0310.061.2018059005952.hdf'

# 读取经度数组
lon_ds = gdal.Open(f'HDF4_EOS:EOS_SWATH:"{filename}":mod04:Longitude')
lon_array = lon_ds.ReadAsArray()

# 读取纬度数组
lat_ds = gdal.Open(f'HDF4_EOS:EOS_SWATH:"{filename}":mod04:Latitude')
lat_array = lat_ds.ReadAsArray()

# 验证维度是否匹配
print(f"经纬度维度：{lon_array.shape}")  # 应输出(676, 451)

方法2：查看完整子数据集列表获取路径 #

如果不确定数据集路径，可以用gdal.Info()输出文件的完整元数据，从中找到经纬度对应的子数据集路径：

# 以JSON格式输出文件完整信息
info = gdal.Info(filename, format='json')

# 遍历所有子数据集，找到经纬度项
for subds in info['subDatasets']:
    print(f"路径：{subds['name']}，描述：{subds['description']}")

从输出中找到Longitude和Latitude对应的路径后，再用gdal.Open()读取即可。

后续地图绘制思路 #

拿到目标变量（如Solar_Zenith）和对应的经纬度数组后，可结合matplotlib、cartopy完成地图绘制：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# 读取目标变量（以Solar_Zenith为例）
solar_ds = gdal.Open(f'HDF4_EOS:EOS_SWATH:"{filename}":mod04:Solar_Zenith')
solar_data = solar_ds.ReadAsArray()

# 创建带地理投影的绘图对象
fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={'projection': ccrs.PlateCarree()})

# 绘制数据，指定坐标投影
im = ax.pcolormesh(lon_array, lat_array, solar_data, transform=ccrs.PlateCarree(), cmap='viridis')

# 添加海岸线和网格
ax.coastlines(resolution='10m')
ax.gridlines(draw_labels=True)

# 添加颜色条和标题
plt.colorbar(im, ax=ax, label='Solar Zenith Angle (degrees)')
plt.title('MYD04_3K Solar Zenith Angle')
plt.show()

内容的提问来源于stack exchange，提问作者Puneet Verma