使用rioxarray将NASA IMERG数据重投影至GK-2A兰伯特投影

核心问题排查方向 #

重投影后可视化和cartopy预期不符，大概率是这几个环节出了问题：

• GK-2A投影参数定义不精准
• 重投影时目标网格的范围、分辨率和GK-2A数据不匹配
• rioxarray重采样方法选错
• 可视化时坐标系映射错误

步骤1：确认GK-2A兰伯特投影的精准参数 #

GK-2A的兰伯特正形投影参数必须完全对应，不能有偏差，标准参数如下（替换成你手里的官方参数）：

gk2a_proj = "+proj=lcc +lat_1=30 +lat_2=60 +lat_0=38 +lon_0=127 +x_0=2100000 +y_0=600000 +ellps=WGS84 +units=m +no_defs"

重点：要和GK-2A数据元数据里的投影参数完全一致，包括标准纬线顺序、中央经线、假东/北坐标、椭球模型

步骤2：修正rioxarray重投影代码 #

手动定义网格容易出错，建议直接匹配GK-2A的网格，或者精准定义目标范围：

import rioxarray
import xarray as xr

# 读IMERG数据，确保坐标系正确
imerg_ds = xr.open_dataset("your_imerg_data.nc")
imerg_ds = imerg_ds.rio.set_crs("EPSG:4326")  # IMERG默认是WGS84地理坐标系

# 替换成你的GK-2A投影参数
gk2a_crs = "+proj=lcc +lat_1=30 +lat_2=60 +lat_0=38 +lon_0=127 +x_0=2100000 +y_0=600000 +ellps=WGS84 +units=m +no_defs"

# 方法1：直接匹配GK-2A数据的网格（推荐，避免手动误差）
# gk2a_ds = xr.open_dataset("your_gk2a_data.nc")
# imerg_reproj = imerg_ds.rio.reproject_match(gk2a_ds, resampling="bilinear")

# 方法2：手动定义目标网格（如果没有GK-2A数据）
x_min, x_max = 1500000, 2700000  # 你的感兴趣区域x范围（米）
y_min, y_max = 0, 1200000        # 你的感兴趣区域y范围（米）
resolution = 2000                # GK-2A通常是2km分辨率，按需调整

# 创建目标网格模板，注意y轴是北到南递减
target_ds = xr.Dataset(
    {
        "x": xr.DataArray(
            range(int(x_min), int(x_max)+resolution, resolution),
            dims=["x"], attrs={"units": "meters"}
        ),
        "y": xr.DataArray(
            range(int(y_max), int(y_min)-resolution, -resolution),
            dims=["y"], attrs={"units": "meters"}
        )
    }
).rio.set_crs(gk2a_crs)

# 重投影，降水数据用bilinear或nearest重采样
imerg_reproj = imerg_ds.rio.reproject_match(target_ds, resampling="bilinear")

# 保存结果
imerg_reproj.to_netcdf("imerg_gk2a_reproj.nc")

步骤3：修正可视化代码 #

可视化时最容易踩的坑是坐标系映射错误，正确的做法是：

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

# 用cartopy定义和GK-2A完全一致的投影
gk2a_cartopy_crs = ccrs.LambertConformal(
    central_longitude=127,
    central_latitude=38,
    standard_parallels=(30, 60),
    false_easting=2100000,
    false_northing=600000,
    globe=ccrs.Globe(ellipse='WGS84')
)

fig, ax = plt.subplots(subplot_kw={"projection": gk2a_cartopy_crs})
# 重投影后的IMERG已经是GK-2A坐标系，所以transform设为目标投影
imerg_reproj.precipitation.plot(ax=ax, transform=gk2a_cartopy_crs)
ax.coastlines()
plt.show()

注意：如果误把transform设为EPSG:4326，数据会严重偏移，这是最常见的可视化错误

常见错误修正总结 #

1. 投影参数不匹配：核对GK-2A数据的元数据，确保proj字符串完全一致
2. 网格方向错误：兰伯特投影的y轴通常北高南低，手动创建网格时要倒序生成y值
3. 重采样方法错误：降水数据不要用sum、mean这类会改变降水量级的方法，用bilinear或nearest
4. 可视化transform错误：重投影后的数据是目标坐标系，必须对应cartopy的投影对象

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