如何使用Python或GDAL将NASA GPM HDF5数据转换为GeoTIFF？

看起来你遇到的是GPM DPR这类条带式（swath）HDF5数据的投影变形问题——这类数据的每个像素都对应独立的经纬度，不是规则的地理格网，直接用GCP的经纬度极值做裁剪范围会把不规则的条带强行拉伸成矩形，自然会出现你看到的怪异图案。我来给你分享两种靠谱的解决方法，亲测有效：

方案一：让GDAL自动处理GCP直接投影 #

其实GDAL已经能识别GPM HDF5里的地理控制点信息，不需要手动提取GCP再指定范围。直接用gdalwarp处理子数据集就行，它会自动根据每个像素的GCP信息完成投影：

import subprocess

# 替换成你的文件路径
input_hdf = "2A.GPM.DPR.V9-20240130.20240505-S231235-E004550.057857.V07C.HDF5"
subdataset_path = f"HDF5:{input_hdf}:/FS/CSF/flagShallowRain"
output_tif = "shallow_rain_flag.tif"

# 构建gdalwarp命令
gdal_warp_cmd = [
    "gdalwarp",
    "-t_srs", "EPSG:4326",  # 目标坐标系，这里用WGS84
    "-r", "near",  # 因为是标志位数据，用最近邻重采样避免改变值
    subdataset_path,
    output_tif
]

try:
    subprocess.run(gdal_warp_cmd, check=True, capture_output=True, text=True)
    print(f"转换完成！输出文件：{output_tif}")
except subprocess.CalledProcessError as e:
    print(f"转换出错了，错误信息：{e.stderr}")

方案二：用经纬度数据集插值生成规则格网 #

如果方案一的结果还是不符合预期，你可以利用HDF5里自带的/Longitude和/Latitude数据集，把每个像素的位置和值组合成点数据，再用gdal_grid插值成规则格网的GeoTIFF，这种方法更适合处理不规则分布的条带数据：

import subprocess

input_hdf = "2A.GPM.DPR.V9-20240130.20240505-S231235-E004550.057857.V07C.HDF5"
output_tif = "shallow_rain_flag_grid.tif"
temp_vrt = "temp_point.vrt"

# 先创建一个VRT文件，把数据和经纬度关联成点要素
vrt_content = f"""<OGRVRTDataSource>
    <OGRVRTLayer name="flagShallowRain">
        <SrcDataSource>{input_hdf}</SrcDataSource>
        <SrcLayer>/FS/CSF/flagShallowRain</SrcLayer>
        <GeometryType>wkbPoint</GeometryType>
        <GeometryField encoding="PointFromColumns" x="/Longitude" y="/Latitude"/>
    </OGRVRTLayer>
</OGRVRTDataSource>"""

# 保存临时VRT文件
with open(temp_vrt, "w") as f:
    f.write(vrt_content)

# 用gdal_grid生成规则格网
gdal_grid_cmd = [
    "gdal_grid",
    "-t_srs", "EPSG:4326",
    "-outsize", "1000", "1000",  # 输出图像的宽高，可根据需求调整
    "-a", "nearest:radius1=1:radius2=1:angle=0:min_points=1",  # 最近邻插值
    temp_vrt,
    output_tif
]

try:
    subprocess.run(gdal_grid_cmd, check=True, capture_output=True, text=True)
    print(f"格网转换完成！输出文件：{output_tif}")
except subprocess.CalledProcessError as e:
    print(f"转换出错了，错误信息：{e.stderr}")

为什么之前的方法会出错？ #

你之前手动提取GCP的经纬度极值作为-te参数，相当于把原本不规则的卫星条带强行“拉平”成矩形框，每个像素的位置被错误拉伸，自然会显示出混乱的图案。GPM这类条带数据的核心是每个像素都有独立的地理坐标，必须让GDAL基于这些坐标来处理，而不是强制套一个矩形范围。

额外注意事项 #

• 确保你的GDAL版本是3.x以上，对HDF5的支持更完善；
• 可以用gdalinfo {input_hdf}命令查看所有子数据集，确认你要处理的路径是否正确；
• 如果处理的是连续数据（比如降水率），可以把重采样方法换成bilinear或cubic，但分类/标志位数据一定要用near最近邻，避免改变原始值。

备注：内容来源于stack exchange，提问作者JWB