Python操作卫星数据时格式解析与地理投影变换的算法实现

Python操作卫星数据：从格式解析到地理投影的实战指南

大家好，作为一名长期与遥感数据打交道的开发者，我深知处理卫星数据时，格式解析和投影变换是两道绕不开的“坎”。新手往往会被HDF、NetCDF等格式搞得晕头转向，或者在将数据可视化时发现图像“歪了”或位置不对。今天，我就结合自己的踩坑经验，带大家用Python一步步攻克这些难题，实现从原始数据到正确空间信息的完整流程。

第一步：理解卫星数据的常见格式与解析库

卫星数据不像普通的图片或文本，它通常包含多维数据（如不同波段）、丰富的元数据（如拍摄时间、轨道参数）和地理空间信息。最常见的格式是HDF（特别是HDF4/HDF5）和NetCDF。刚开始时，我曾试图用文本编辑器打开它们，结果自然是乱码一片。

实战选择： 在Python生态中，我们主要依赖以下几个库：

• h5py / pyhdf：用于读取HDF5和HDF4格式。HDF5现在更主流，h5py的接口也非常友好。
• netCDF4：用于读取NetCDF格式，这是气象和海洋卫星数据的常客。
• rasterio / GDAL：更通用的地理空间数据处理库，它们不仅能读格式，还能直接处理投影。对于初学者，rasterio的API更简洁。
• xarray：一个高级工具，特别适合处理带标签的多维数组（比如带经纬度坐标的栅格数据），它能优雅地封装NetCDF数据。

让我们先从一个具体的HDF5文件（比如MODIS卫星的L2级数据）开始。假设我们已经下载了文件 MOD021KM.A2024120.0230.061.hdf。

# 首先安装必要的库
pip install h5py numpy matplotlib rasterio

import h5py
import numpy as np

# 打开HDF5文件
file_path = 'MOD021KM.A2024120.0230.061.hdf'
with h5py.File(file_path, 'r') as f:
    # 首先，查看文件结构，这是关键的一步！
    def print_structure(name, obj):
        print(name, type(obj))
    f.visititems(print_structure) # 这会打印出所有数据集和组的路径

    # 假设我们找到了反射率数据集的路径，例如 'EV_250_Aggr1km_RefSB'
    dataset_path = 'EV_250_Aggr1km_RefSB'
    if dataset_path in f:
        data = f[dataset_path][:]  # 将数据读入NumPy数组
        # 读取缩放因子和偏移量（卫星数据常用存储方式：实际值 = 缩放因子 * 存储值 + 偏移量）
        scale_factor = f[dataset_path].attrs.get('scale_factor', 1.0)
        add_offset = f[dataset_path].attrs.get('add_offset', 0.0)
        calibrated_data = data.astype(np.float32) * scale_factor + add_offset
        print(f"数据形状: {calibrated_data.shape}, 数据类型: {calibrated_data.dtype}")

踩坑提示： 直接读取的数据数组往往是整型（为了压缩存储），必须根据属性（`scale_factor`, `add_offset`）进行定标换算才能得到真实的物理值（如反射率）。忽略这一步，后续所有分析都是错的。

第二步：提取地理定位信息与理解原始投影

拿到数据数组只是第一步，它只是一堆数字。要让这些数字对应到地球上的位置，我们需要地理定位信息。通常，同一个HDF文件里会有独立的经纬度数据集，名字可能是`Latitude`, `Longitude`或`geolocation/lat`等。

# 继续在同一个with语句块中操作
    # 查找并读取经纬度数据
    lat_path = 'Latitude'
    lon_path = 'Longitude'
    if lat_path in f and lon_path in f:
        latitude = f[lat_path][:]
        longitude = f[lon_path][:]
        print(f"经纬度网格形状: 经度 {longitude.shape}, 纬度 {latitude.shape}")

现在，我们有了`calibrated_data`（比如某个波段的反射率）和对应的`latitude`、`longitude`网格。但请注意，很多卫星数据（尤其是中低分辨率全球数据）的原始投影不是常见的经纬度（WGS84）。MODIS数据的Level-1B产品通常使用正弦曲线投影（Sinusoidal）。这意味着`latitude`和`longitude`数组本身是地理坐标，但数据像素的排列是在正弦投影的平面坐标系下。直接使用经纬度网格进行拼接或与GIS软件中的WGS84图层叠加，可能会产生微小错位。

第三步：核心算法——地理投影变换

投影变换的本质，是将数据从一个坐标系统（源投影）转换到另一个坐标系统（目标投影）。我们最常用的目标投影是地理坐标系WGS84（EPSG:4326）或网络墨卡托Web Mercator（EPSG:3857，用于在线地图）。

这里我推荐使用rasterio和pyproj来完成。我们需要知道源投影的信息。对于MODIS正弦投影，其PROJ字符串是固定的。

from rasterio.transform import from_origin
from rasterio.crs import CRS
import pyproj

# 假设我们从元数据中获知了图像左上角在正弦投影下的坐标和像素分辨率
# 这些信息有时存储在HDF文件的全局属性中，需要仔细查阅产品文档
# 这里以假设值演示
sinu_ulx = -20015109.354  # 左上角X坐标（米）
sinu_uly = 10007554.677   # 左上角Y坐标（米）
pixel_size = 1000.0       # 像素大小1000米

# 1. 定义源投影（MODIS正弦投影）
src_crs = CRS.from_proj4("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m")
# 2. 定义目标投影（例如WGS84）
dst_crs = CRS.from_epsg(4326)  # WGS84

# 3. 创建转换器
transformer = pyproj.Transformer.from_crs(src_crs, dst_crs, always_xy=True)

# 4. 转换图像四个角点的坐标（从源投影到地理坐标）
# 计算图像右下角坐标
height, width = calibrated_data.shape  # 假设是单波段
sinu_lrx = sinu_ulx + width * pixel_size
sinu_lry = sinu_uly - height * pixel_size

# 转换左上角和右下角
ul_lon, ul_lat = transformer.transform(sinu_ulx, sinu_uly)
lr_lon, lr_lat = transformer.transform(sinu_lrx, sinu_lry)
print(f"转换后地理范围: 左上({ul_lat:.2f}, {ul_lon:.2f}), 右下({lr_lat:.2f}, {lr_lon:.2f})")

更实用的方法：使用rasterio进行重投影与重采样
上面的方法转换了角点，但我们需要的是整个图像网格在新的投影下重新排列。这需要重采样。我们可以用rasterio.warp.reproject函数。

from rasterio.warp import calculate_default_transform, reproject, Resampling

# 计算从源投影到目标投影的最优变换参数和新图像尺寸
transform, width_dst, height_dst = calculate_default_transform(
    src_crs, dst_crs, width, height,
    left=sinu_ulx, right=sinu_lrx, top=sinu_uly, bottom=sinu_lry
)

# 准备一个空数组来存放重投影后的数据
data_reprojected = np.zeros((height_dst, width_dst), dtype=np.float32)

# 执行重投影
reproject(
    source=calibrated_data,          # 源数据数组
    destination=data_reprojected,    # 目标数组
    src_transform=from_origin(sinu_ulx, sinu_uly, pixel_size, pixel_size), # 源图像的空间变换
    src_crs=src_crs,                 # 源投影
    dst_transform=transform,         # 计算得到的目标图像变换
    dst_crs=dst_crs,                 # 目标投影
    resampling=Resampling.nearest    # 重采样方法，最近邻法保真度，双线性插值更平滑
)
print(f"重投影后数据形状: {data_reprojected.shape}")

实战经验： 选择Resampling.nearest（最近邻）可以保持原始值不变，适合分类产品；选择Resampling.bilinear（双线性）会使图像更平滑，适合连续变量（如反射率、温度）。重投影是计算密集型操作，大数据时比较耗时。

第四步：验证与输出

完成变换后，如何验证？最简单的办法是用matplotlib叠加一个海岸线或与已知的WGS84矢量数据对比。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs  # 一个强大的地理绘图库

fig = plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())  # 创建WGS84坐标系的地图轴
ax.coastlines(resolution='10m', color='red', linewidth=0.5)  # 添加海岸线

# 显示重投影后的数据，并指定其坐标系（现在已经是WGS84）
extent = [ul_lon, lr_lon, lr_lat, ul_lat]  # matplotlib的extent顺序是 [左， 右， 下， 上]
img = ax.imshow(data_reprojected, extent=extent, origin='upper',
                cmap='viridis', transform=ccrs.PlateCarree())
plt.colorbar(img, ax=ax, label='反射率')
ax.gridlines(draw_labels=True)  # 添加经纬度网格线
plt.title('MODIS数据重投影至WGS84')
plt.show()

如果海岸线与你的图像轮廓完美贴合，那么恭喜你，投影变换成功了！最后，我们可以将处理好的数据保存为更通用的格式，如GeoTIFF。

import rasterio

# 使用计算得到的目标变换参数和CRS保存为GeoTIFF
output_profile = {
    'driver': 'GTiff',
    'height': height_dst,
    'width': width_dst,
    'count': 1,  # 单波段
    'dtype': data_reprojected.dtype,
    'crs': dst_crs,
    'transform': transform,
    'compress': 'LZW'  # 压缩以减小文件体积
}

with rasterio.open('reprojected_data.tif', 'w', **output_profile) as dst:
    dst.write(data_reprojected, 1)
print("GeoTIFF文件已保存，可在QGIS/ArcGIS中直接打开。")

总结与核心要点

回顾整个流程：格式解析 -> 数据定标 -> 定位信息提取 -> 投影定义 -> 坐标变换/重投影 -> 验证输出。每一步都至关重要。

最后的忠告：
1. 永远先读元数据：投影参数、定标参数、填充值等信息都藏在文件属性里，不读明白寸步难行。
2. 理解你的数据：清楚你处理的是哪颗卫星、哪个产品、哪个级别（L1B, L2, L3），它们的默认投影是什么。
3. 善用专业库：不要试图自己写底层投影算法，pyproj、rasterio、GDAL是经过千锤百炼的工具。
4. 可视化验证：人眼是最后的裁判，一张叠加了基础地图的图片能立刻告诉你对错。

处理卫星数据就像解谜，虽然初期会遇到各种格式和坐标系的“拦路虎”，但一旦掌握了这套方法，你就能自由地让卫星数据在正确的位置“说话”了。希望这篇指南能帮你少走弯路，祝大家探索愉快！