GDAL中旋转经纬度网格转EPSG:4326投影失败求助

处理带旋转经纬度网格的NetCDF文件时，投影转换后在QGIS中无法与意大利等地物正确对齐，尝试多种GDAL命令及Python（xarray、numpy、pyproj）方法均未解决，具体细节如下：

数据摘要 #

Dimensions:                         (time: 24, rlat: 652, rlon: 570, bnds: 2)
Coordinates:
  * time                            (time) datetime64[ns] 192B 2006-01-01 ......
  * rlon                            (rlon) float64 5kB -5.913 -5.893 ... 5.467
  * rlat                            (rlat) float64 5kB -6.408 -6.388 ... 6.612
    lat                             (rlat, rlon) float32 1MB ...
    lon                             (rlat, rlon) float32 1MB ...
Dimensions without coordinates: bnds
Data variables:
    TOT_PREC                        (time, rlat, rlon) float32 36MB ...
    crs_rotated_latitude_longitude  int64 8B ...
    time_bnds                       (time, bnds) datetime64[ns] 384B ...
Attributes: (12/13)
    CDI:             Climate Data Interface version 1.9.8 (https://mpimet.mpg...
    Conventions:     CF-1.4
    history:         Generated by CMCC DDS version 0.9.0 2024-07-29 08:33:07....
    source:          COSMO
    institution:     CMCC (Euro-Mediterranean Center on Climate Change), REMH...
    title:           cclm-sp_2.4_terra_urb_2.3.1 simulation (0.02 Deg) forced...
    ...              ...
    experiment_id:   Historical with Urban Parametrization, DT=20sec
    conventionsURL:  http://www.cfconventions.org/
    contact:         Mario Raffa (mario.raffa@cmcc.it)
    references:      http//www.clm-community.eu, http://www.cmcc.it
    creation_date:   2021-12-31 12:05:48
    CDO:             Climate Data Operators version 1.9.8 (https://mpimet.mpg...

尝试的方法及错误 #

1. GDAL命令指定SRS #

执行命令：

gdal_translate -a_srs "EPSG:4326" precipitation.nc precipitation_with_srs.nc

错误提示：NetCDF: Index exceeds dimension bound.

2. Python代码转换 #

执行代码：

import xarray as xr
import numpy as np
from pyproj import Transformer
import os

# Path to the original NetCDF file
input_file_path = r'C:\Users\HUNDESA\Desktop\hun - Copy\precipitation_amount_era5-downscaled-over-italy_hourly_418082.nc'

# Path for the new NetCDF file
output_file_path = r'C:\Users\HUNDESA\Desktop\hun - Copy\transformed_file.nc'

# Open the original NetCDF file using xarray
ds = xr.open_dataset(input_file_path)

# Rotation parameters
north_pole_lat = 47
north_pole_lon = -168

# Define the coordinate system transformation
transformer = Transformer.from_proj(
    proj_from=f"+proj=aeqd +lat_0={north_pole_lat} +lon_0={north_pole_lon}",
    proj_to="epsg:4326"  # WGS84
)

# Create a copy of the dataset
new_ds = ds.copy()

# Transform the coordinates if necessary


# Save the new dataset to a new NetCDF file
encoding = {
    'TOT_PREC': {'dtype': 'float32'},
    'latitude': {'dtype': 'float32'},
    'longitude': {'dtype': 'float32'},
}

# Convert DataArray to DataArray with proper dtype if needed
for var_name in encoding.keys():
    if var_name in new_ds.variables:
        new_ds[var_name] = new_ds[var_name].astype(encoding[var_name]['dtype'])

# Save the dataset
new_ds.to_netcdf(output_file_path, encoding=encoding)

问题：无法将数据正确投影至EPSG:4326，可视化结果扭曲。

内容的提问来源于stack exchange，提问作者Hundesa Siraj Mohammed