使用 rioxarray 将 xarray 转换为 3D 并导出为 tif 文件

transform xarray to 3D and export as tif file with rioxarray

提问人:Giacomo 提问时间:7/11/2023 更新时间:7/11/2023 访问量:82

问:

我有一个多维 xarray。我想为调用的维度选择一个特定值,然后删除它,只保留 、 和 维度并将 xarray 导出为变量的 tif 文件。numbersteplatitudelongitude

这是我的起始xarray:

<xarray.Dataset>
Dimensions:            (number: 2, step: 2, latitude: 45, longitude: 45)
Coordinates:
  * number             (number) int32 1 2
    time               datetime64[ns] ...
  * step               (step) timedelta64[ns] 06:00:00 12:00:00
    heightAboveGround  float64 ...
  * latitude           (latitude) float64 28.8 28.4 28.0 27.6 ... 12.0 11.6 11.2
  * longitude          (longitude) float64 111.6 112.0 112.4 ... 128.8 129.2
    valid_time         (step) datetime64[ns] dask.array<chunksize=(2,), meta=np.ndarray>
Data variables:
    v10                (number, step, latitude, longitude) float32 dask.array<chunksize=(2, 2, 15, 15), meta=np.ndarray>
    u10                (number, step, latitude, longitude) float32 dask.array<chunksize=(2, 2, 15, 15), meta=np.ndarray>
Attributes:
    GRIB_edition:            2
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2023-07-11T10:38 GRIB to CDM+CF via cfgrib-0.9.1...

在这里,我屏蔽了 xarray,以便只返回一个,并且我试图压缩维度,但它仍然存在。number

ensemble = 1
mask = (xarr.number.values == ensemble)
xarr1 = xarr.sel(number=mask).squeeze(drop=True)
xarr1.rio.write_crs("epsg:4326", inplace=True) # set wgs84 proj to xarray
<xarray.Dataset>
Dimensions:            (number: 0, step: 2, latitude: 45, longitude: 45)
Coordinates:
  * number             (number) int32 
    time               datetime64[ns] ...
  * step               (step) timedelta64[ns] 06:00:00 12:00:00
    heightAboveGround  float64 ...
  * latitude           (latitude) float64 28.8 28.4 28.0 27.6 ... 12.0 11.6 11.2
  * longitude          (longitude) float64 111.6 112.0 112.4 ... 128.8 129.2
    valid_time         (step) datetime64[ns] dask.array<chunksize=(2,), meta=np.ndarray>
    spatial_ref        int32 0
Data variables:
    v10                (number, step, latitude, longitude) float32 dask.array<chunksize=(0, 2, 15, 15), meta=np.ndarray>
    u10                (number, step, latitude, longitude) float32 dask.array<chunksize=(0, 2, 15, 15), meta=np.ndarray>
Attributes:
    GRIB_edition:            2
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2023-07-11T10:56 GRIB to CDM+CF via cfgrib-0.9.1...

我计算了风速,并将其作为新变量添加到xarray中ws

ds2 = xarr1.assign(ws=(xarr1["u10"] ** 2 + xarr1["v10"] ** 2) ** 0.5)

我创建了一个仅包含风速变量的新 xarray

xarr_ws = ds2.ws
<xarray.DataArray 'ws' (number: 0, step: 2, latitude: 45, longitude: 45)>
dask.array<pow, shape=(0, 2, 45, 45), dtype=float32, chunksize=(0, 2, 24, 24), chunktype=numpy.ndarray>
Coordinates:
  * number             (number) int32 
    time               datetime64[ns] ...
  * step               (step) timedelta64[ns] 06:00:00 12:00:00
    heightAboveGround  float64 ...
  * latitude           (latitude) float64 28.8 28.4 28.0 27.6 ... 12.0 11.6 11.2
  * longitude          (longitude) float64 111.6 112.0 112.4 ... 128.8 129.2
    valid_time         (step) datetime64[ns] dask.array<chunksize=(2,), meta=np.ndarray>
    spatial_ref        int32 0

使用 rioxarray 导出到 tif 文件,但这失败了,因为我有 4 个维度而不是 3 个维度。

# save tif with wind speed with all timesteps
xarr_ws.rio.to_raster(f'{storm_name[1:]}_ens{ensemble}_wind_speed_by_timestep.tif')

错误信息:

in _check_dimensions raise TooManyDimensions(rioxarray.exceptions.TooManyDimensions: Only 2D and 3D data arrays supported. Data variable: ws

如何将尺寸从 减小到 ?(number: 0, step: 2, latitude: 45, longitude: 45)(step: 2, latitude: 45, longitude: 45)

我试过了,但没有运气.squeeze('number').drop(dim='number').drop_dims(number)

多维数组 python-xarray 维度

评论

1赞 cyril 7/11/2023
在屏蔽操作之后,拍摄了 0 个维度编号的切片,我相信这不是您想要的数据集。
0赞 Giacomo 7/11/2023
谢谢你@cyril很好。我的代码实际上是正确的,问题是我使用的值不是维度的一部分。ensemblenumber

答: 暂无答案