INSAT静止卫星LST的H5文件重投影至EPSG:4326问题求助

问题分析 #

你遇到的重投影偏移问题，核心原因是手动指定的GEOS投影参数不全，且未正确将原始的GeoX/GeoY转换为地球同步投影的实际坐标（原始的GeoX/GeoY是像素索引，并非投影平面坐标）。

修正步骤与代码 #

1. 完善GEOS投影参数 #

INSAT-3D的地球同步投影需要完整的椭球参数，替换你之前的CRS定义：

geos_crs = "+proj=geos +lon_0=82 +h=35785831 +x_0=0 +y_0=0 +a=6378137 +b=6356752.3142 +units=m +no_defs"

2. 正确映射像素坐标到GEOS投影 #

原始的GeoX/GeoY是像素位置，需要转换为GEOS投影下的实际坐标。根据文件属性的覆盖范围，计算每个像素对应的投影坐标：

import os
import warnings
import numpy as np
import xarray as xr
import rioxarray as rxr
import rasterio as rio
warnings.simplefilter('ignore')

# 读取文件
a2 = xr.open_dataset('/home/data/lab_d/../ncf/3DIMG_19APR2015_1100_L2B_LST.h5')

# 获取目标变量的图像尺寸（以LST为例，替换为你的实际变量名）
height, width = a2['LST'].shape

# 定义完整的GEOS投影参数
geos_crs = "+proj=geos +lon_0=82 +h=35785831 +x_0=0 +y_0=0 +a=6378137 +b=6356752.3142 +units=m +no_defs"

# 计算GEOS投影下的坐标范围
orbit_radius = 35785831
x_res = (2 * np.pi * orbit_radius) / width
y_res = x_res  # 方形像素
x_min = -width * x_res / 2
x_max = width * x_res / 2
y_min = -height * y_res / 2
y_max = height * y_res / 2

# 创建投影坐标数组（注意y轴从上到下对应纬度递减）
x_coords = np.linspace(x_min, x_max, width)
y_coords = np.linspace(y_max, y_min, height)

# 替换原始的GeoX/GeoY为投影坐标
a2 = a2.assign_coords(x=x_coords, y=y_coords)

# 写入CRS信息
a2 = a2.rio.write_crs(geos_crs)

# 重投影到EPSG:4326，由rioxarray自动计算输出网格
ins_lonlat = a2.rio.reproject(dst_crs="EPSG:4326")

# 查看结果
print(ins_lonlat)

3. 关键说明 #

• 投影参数完整性：必须包含椭球长半轴a和短半轴b，否则重投影计算会出现系统性偏差。
• 坐标转换：原始的GeoX/GeoY是像素索引，不是投影坐标，必须转换为GEOS投影下的米制坐标才能正确重投影。
• 避免手动指定shape：手动指定shape可能导致分辨率不匹配，让rioxarray根据输入数据自动计算输出的网格尺寸和分辨率。

验证方法 #

可以对比重投影后的经纬度范围与文件属性中的left_longitude/right_longitude/upper_latitude/lower_latitude，确认结果是否匹配。

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