如何将具有异常尺寸的 netCDF 转换为标准的 netCDF（ltime、lat、lon）（python）

Question

我有多个最终要合并的 netCDF 文件。一个示例 netCDF 如下所示。

import xarray as xr
import numpy as np
import cftime

Rain_nc = xr.open_dataset('filepath.nc', decode_times=False)
print(Rain_nc)

<xarray.Dataset>
Dimensions: (land: 67209, tstep:248)
Dimensions without coordinates: land, tstep
Data variables:
    lon    (land) float32...
    lat    (land) float32...
    timestp(tstep) int32...
    time   (tstep) int32...
    Rainf  (tstep, land) float32...

维度“land”是从 1 到 67209 的计数，“tstep”是从 1 到 248 的计数。

变量 'lat' 和 'lon' 是纬度和经度值，形状为 (67209,)

变量“时间”是自月初以来的秒数（netcdf 为一个月）

接下来我将维度从“tstep”交换为“time”，将其转换为以后合并并将坐标设置为 lon、lat 和 time。

rain_nc = rain_nc.swap_dims({'tstep':'time'})
rain_nc = rain_nc.set_coords(['lon', 'lat', 'time'])

rain_nc['time'] = cftime.num2date(rain_nc['time'], units='seconds since 2016-01-01 00:00:00', calendar = 'standard')
rain_nc['time'] = cftime.date2num(rain_nc['time'], units='seconds since 1970-01-01 00:00:00', calendar = 'standard')

这给我留下了以下数据集：

print(rain_nc)

<xarray.Dataset>
Dimensions: (land: 67209, time: 248)
Coordinates:
    lon        (land)float32
    lat        (land)float32
  * time       (time)float64
Dimensions without coordinates: land
Data variables:
    timestp   (time)int32
    Rainf     (time, land)


print(rain_nc['land'])
<xarray.DataArray 'land' (land: 67209)>
array([    0,    1,    2,..., 67206, 67207, 67208])
Coordinates:
    lon     (land) float32 ...
    lat     (land) float32 ...
Dimensions without coordinates: land

我感兴趣的Rainf变量如下：

<xarray.DataArray 'Rainf' (time: 248, land: 67209)>
[16667832 values with dtype=float32]
Coordinates:
    lon      (land) float32 -179.75 -179.75 -179.75 ... 179.75 179.75 
179.75
    lat      (land) float32 71.25 70.75 68.75 68.25 ... -16.25 -16.75 
-19.25
  * time     (time) float64 1.452e+09 1.452e+09 ... 1.454e+09 1.454e+09
Dimensions without coordinates: land
Attributes:
    title:       Rainf
    units:       kg/m2s
    long_name:   Mean rainfall rate over the \nprevious 3 hours
    actual_max:  0.008114143
    actual_min:  0.0
    Fill_value:  1e+20

从这里我想创建一个带有维度（时间、纬度、经度）和变量 Rainf 的 netCDF。

我尝试创建一个新的 netCDF（或更改这个），但是当我尝试传递 Rainf 变量时不起作用，因为它的形状为 (248, 67209) 并且需要形状为 (248, 67209, 67209 ）。尽管“Rainf”的当前“土地”维度具有纬度和经度坐标。是否可以将此变量重塑为具有时间、纬度和经度维度？

Answer 1

最后，您似乎想要将"land" 尺寸重塑为("lat", "lon") 尺寸。

所以，你有一些类似这样的 DataArray：

# Setting sizes and coordinates
lon_size, lat_size = 50, 80                                                                                                                                                                           
lon, lat = [arr.flatten() for arr in np.meshgrid(range(lon_size), range(lat_size))]                                                                                                                   
land_size = lon_size * lat_size                                                                                                                                                                       
time_size = 100 

da = xr.DataArray( 
    dims=("time", "land"), 
    data=np.random.randn(time_size, land_size), 
    coords=dict( 
        time=np.arange(time_size), 
        lon=("land", lon), 
        lat=("land", lat), 
    ) 
)  

看起来像这样：

>>> da
<xarray.DataArray (time: 100, land: 4000)>
array([[...]])
Coordinates:
  * time     (time) int64 0 1 ... 98 99
    lon      (land) int64 0 1 ... 48 49
    lat      (land) int64 0 0 ... 79 79
Dimensions without coordinates: land

首先，我们将使用.set_index() 方法告诉xarray "land" 索引应该由"lon" 和"lat" 坐标表示：

>>> da.set_index(land=("lon", "lat"))                                                                                                                                                                    
<xarray.DataArray (time: 100, land: 4000)>
array([[...]])
Coordinates:
  * time     (time) int64 0 1 ... 98 99
  * land     (land) MultiIndex
  - lon      (land) int64 0 1 ... 48 49
  - lat      (land) int64 0 0 ... 79 79

尺寸仍然是("time", "land")，但现在"land" 是MultiIndex。

请注意，如果您此时尝试写入 NETCDF，您将遇到以下错误：

>>> da.set_index(land=("lon", "lat")).to_netcdf("data.nc")   
NotImplementedError: variable 'land' is a MultiIndex, which cannot yet be serialized to netCDF files (https://github.com/pydata/xarray/issues/1077). Use reset_index() to convert MultiIndex levels into coordinate variables instead.

它告诉你使用.reset_index() 方法。但这不是你想要的，因为它只会回到原来的da 状态。

你现在想要的是使用.unstack()方法：

>>> da.set_index(land=("lon", "lat")).unstack("land")                                                                                                                                                    
<xarray.DataArray (time: 100, lon: 50, lat: 80)>
array([[[...]]])
Coordinates:
  * time     (time) int64 0 1 ... 98 99
  * lon      (lon) int64 0 1 ... 48 49
  * lat      (lat) int64 0 1 ... 78 79

它有效地杀死了"land" 维度并提供了所需的输出。