使用Rioxarray将北极立体投影转经纬度致数据集损坏问题

数据老陈

2026-8-5

问题

我尝试将GeoTIFF格式的CMC雪深图从北极立体投影（North Polar Stereo）转换为经纬度（Lat/Lon），该文件的PROJ字符串为：+proj=stere +lat_0=90 +lat_ts=60 +lon_0=10 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs=True。

以下是我的代码：

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling
from rasterio.crs import CRS
from rioxarray.crs import crs_from_user_input
import rioxarray as rxr
#### 读取CMC数据 ####
yrs = np.arange(1980,2020,1)

CMC_dir = '/mnt/data/obs/CMC_snow_depth/'

CMC_fil = ''.join([CMC_dir+'cmc_sdepth_dly_1998_v01.2.tif'])

dates = pd.date_range(start='08-01-1998', end='12-31-1998', freq='D')

### 创建XArray数据集 ###
CMC_dataset = rxr.open_rasterio(CMC_fil,decode_coords="all")
print(CMC_dataset)

# 写入CRS，使用proj4字符串
crs_CMC = CRS.from_string('+proj=stere +lat_0=90 +lat_ts=60 +lon_0=10 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs=True')
CMC_dataset = CMC_dataset.rio.write_crs(crs_CMC)

    
### 从极地立体投影重投影到经纬度 ###
CMC_dataset_latlon = CMC_dataset.rio.reproject("EPSG:4326") #EPSG:4326是经纬度投影
print(CMC_dataset_latlon)

输入数据结构：

<xarray.DataArray (band: 153, y: 706, x: 706)>
[76260708 values with dtype=float64]
Coordinates:
  * band         (band) int64 1 2 3 4 5 6 7 8 ... 147 148 149 150 151 152 153
  * x            (x) float64 -8.394e+06 -8.37e+06 ... 8.37e+06 8.394e+06
  * y            (y) float64 8.394e+06 8.37e+06 ... -8.37e+06 -8.394e+06
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     -1.7e+308
    scale_factor:   1.0
    add_offset:     0.0

但重投影到经纬度后，数据集出现异常，经纬度数值不对：

<xarray.DataArray (band: 153, y: 54, x: 997)>
array([[[ 0.0e+000,  0.0e+000,  0.0e+000, ...,  0.0e+000,  0.0e+000,
          0.0e+000],
        [ 0.0e+000,  0.0e+000,  0.0e+000, ...,  0.0e+000,  0.0e+000,
          0.0e+000],
        [-1.7e+308, -1.7e+308, -1.7e+308, ...,  0.0e+000,  0.0e+000,
         -1.7e+308],
        ...,
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308],
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308],
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308]],

       [[ 0.0e+000,  0.0e+000,  0.0e+000, ...,  0.0e+000,  0.0e+000,
          0.0e+000],
        [ 0.0e+000,  0.0e+000,  0.0e+000, ...,  0.0e+000,  0.0e+000,
          0.0e+000],
        [-1.7e+308, -1.7e+308, -1.7e+308, ...,  0.0e+000,  0.0e+000,
         -1.7e+308],
...
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308],
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308],
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308]],

       [[ 0.0e+000,  0.0e+000,  0.0e+000, ...,  0.0e+000,  0.0e+000,
          0.0e+000],
        [ 0.0e+000,  0.0e+000,  0.0e+000, ...,  0.0e+000,  0.0e+000,
          0.0e+000],
        [-1.7e+308, -1.7e+308, -1.7e+308, ...,  0.0e+000,  0.0e+000,
         -1.7e+308],
        ...,
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308],
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308],
        [-1.7e+308, -1.7e+308, -1.7e+308, ..., -1.7e+308, -1.7e+308,
         -1.7e+308]]])
Coordinates:
  * x            (x) float64 -179.8 -179.5 -179.1 -178.7 ... 179.1 179.4 179.8
  * y            (y) float64 19.3 18.94 18.57 18.21 ... 0.8805 0.5194 0.1583
  * band         (band) int64 1 2 3 4 5 6 7 8 ... 147 148 149 150 151 152 153
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0
    _FillValue:     -1.7e+308

请问这是什么原因导致的？

原因分析与解决方案

核心原因

CRS匹配冲突：手动写入的CRS可能和原始GeoTIFF自带的CRS信息冲突，或者proj字符串的参数细节（如轴方向、椭球定义）与数据实际存储的坐标系统不匹配，导致重投影时地理范围计算错误。
默认重投影的范围/分辨率不合理：rio.reproject自动计算的输出范围没有针对北极区域优化，错误地裁剪掉了高纬度数据，只保留了低纬度区域。
极端填充值干扰：原始数据的填充值-1.7e+308数值过大，重投影时插值逻辑无法正确识别，导致填充值区域扩散覆盖有效数据。

解决步骤

先验证原始数据的CRS
先检查文件是否自带正确CRS，避免手动写入导致冲突：
```
print(CMC_dataset.rio.crs)
```
如果输出已有正确的北极立体投影信息，直接跳过手动写入CRS的步骤。

手动指定重投影的目标范围与分辨率
针对北极区域定义目标经纬度范围（比如北纬60°到90°），并设置合适的分辨率：

# 目标范围：经度-180°到180°，纬度60°N到90°N
target_bounds = (-180, 60, 180, 90)
# 目标分辨率（0.5°，可根据需求调整）
target_res = (0.5, 0.5)

CMC_dataset_latlon = CMC_dataset.rio.reproject(
    "EPSG:4326",
    bounds=target_bounds,
    resolution=target_res,
    resampling=Resampling.nearest  # 雪深数据适合用最近邻插值
)

正确处理填充值
将极端填充值替换为标准NaN，避免插值错误：

# 替换原始填充值为NaN
CMC_dataset = CMC_dataset.where(CMC_dataset != CMC_dataset.rio.nodata, np.nan)
# 重投影时指定NaN为填充值
CMC_dataset_latlon = CMC_dataset.rio.reproject(
    "EPSG:4326",
    bounds=target_bounds,
    resolution=target_res,
    resampling=Resampling.nearest,
    nodata=np.nan
)

验证原始数据的地理范围
用rasterio查看原始数据的实际地理边界，确认CRS是否正确：
```
import rasterio
with rasterio.open(CMC_fil) as src:
    print(src.bounds)
    print(src.crs)
```

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

使用Rioxarray将北极立体投影转经纬度致数据集损坏问题

核心原因 #

解决步骤 #

核心原因

解决步骤