Xarray将含二维经纬度的UKCP18 NetCDF4文件转TIFF方法

问题背景 #

• 待处理数据为空间分辨率2.2km的UKCP18气候数据，采用专用旋转极地理参考网格，目标为转换为WGS84（EPSG:4326）坐标系的TIFF格式文件
• 数据原生空间维度为grid_latitude、grid_longitude，文件内同时存储二维非维度坐标latitude、longitude，二者均为原生网格维度的二维映射值
• 要求避免使用gdalwarp执行不必要的重投影操作，直接基于文件自带的二维经纬度信息完成空间参考配置，最终得到可直接使用的栅格文件

当前数据结构 #

<xarray.DataArray 'tasmin' (grid_latitude: 606, grid_longitude: 484)>
dask.array<mean_agg-aggregate, shape=(606, 484), dtype=float32, chunksize=(606, 484), chunktype=numpy.ndarray>
Coordinates:
  * grid_latitude   (grid_latitude) float64 -4.683 -4.647 -4.611 ... 8.027 8.063
  * grid_longitude  (grid_longitude) float64 353.9 354.0 354.0 ... 364.3 364.3
    latitude        (grid_latitude, grid_longitude) float64 dask.array<chunksize=(606, 484), meta=np.ndarray>
    longitude       (grid_latitude, grid_longitude) float64 dask.array<chunksize=(606, 484), meta=np.ndarray>

目标 #

• 优先方案：无重采样、无精度损失导出带正确地理参考的TIFF文件
• 可选方案：转换为以一维latitude、longitude为维度的规则网格结构，结构示例如下：

<xarray.DataArray 'tasmin' (latitude: 606, longitude: 484)>

数据有效性验证：执行da.plot(x='longitude', y='latitude')可正常绘制英国区域不规则网格上的最低气温分布图，证明文件自带经纬度坐标有效。

实现步骤 #

依赖安装 #

提前安装所需依赖库：

pip install xarray rioxarray rasterio dask numpy xesmf

方案一：零重采样直接导出曲线网格TIFF（推荐） #

该方案完全基于原始数据自带的二维经纬度写入地理参考，无任何插值、重投影操作，精度无损失，不需要调用gdalwarp。

import xarray as xr
import rioxarray

# 读取目标数据，替换为实际文件路径与变量名
ds = xr.open_dataset("your_ukcp18_data.nc", chunks="auto")
da = ds["tasmin"]

# 预处理：将0~360范围经度转换为-180~180范围，避免坐标错位
da["longitude"] = xr.where(da["longitude"] > 180, da["longitude"] - 360, da["longitude"])

# 写入坐标系信息，绑定二维经纬度为空间坐标
da = da.rio.write_crs("EPSG:4326")
da = da.rio.set_spatial_dims(x_dim="grid_longitude", y_dim="grid_latitude")
da = da.assign_coords({
    "x": da.longitude,
    "y": da.latitude
})

# 导出为地理配准TIFF
da.rio.to_raster(
    "ukcp18_tasmin_wgs84_native.tif",
    recalc_transform=False,
    tiled=True,
    compress="deflate"
)

导出的TIFF为曲线栅格格式，可直接在QGIS、ArcGIS Pro等主流GIS工具中正常加载使用，坐标与原始数据完全一致。

方案二：重采样为规则经纬网格（满足维度替换需求） #

如果必须得到以一维经纬度为维度的规则网格结构，可使用xesmf执行最邻近重采样，重采样误差极小，远优于gdalwarp默认重投影效果。

import numpy as np
import xesmf as xe

# 基于原始数据经纬度范围生成目标规则网格，0.02度分辨率匹配原始2.2km空间精度
lon_min = float(da.longitude.min().values)
lon_max = float(da.longitude.max().values)
lat_min = float(da.latitude.min().values)
lat_max = float(da.latitude.max().values)

target_lon = np.arange(lon_min, lon_max, 0.02)
target_lat = np.arange(lat_min, lat_max, 0.02)
target_grid = xr.Dataset({
    "longitude": (["longitude"], target_lon),
    "latitude": (["latitude"], target_lat)
})

# 执行最邻近重采样，无平滑模糊问题
regridder = xe.Regridder(da, target_grid, method="nearest_s2d", periodic=False)
da_regular = regridder(da)

# 导出规则网格TIFF
da_regular = da_regular.rio.write_crs("EPSG:4326")
da_regular.rio.to_raster(
    "ukcp18_tasmin_wgs84_regular.tif",
    tiled=True,
    compress="deflate"
)

执行完成后da_regular即为(latitude: xxx, longitude: xxx)结构的规则网格DataArray，完全符合目标结构要求。

注意事项 #

• 原始数据自带的二维经纬度已经是EPSG:4326坐标系下的真实坐标，不需要额外做投影参数转换
• 若不需要规则经纬网格，优先选择方案一，完全规避重采样带来的精度损失
• 处理大体积NetCDF文件时，open_dataset添加chunks="auto"参数启用dask分块加载，避免内存溢出

内容的提问来源于stack exchange，提问作者Kieran Sam Bajpai