首页 > 解决方案 > xarray concat中的空切片

问题描述

我有一堆 HDF5 文件。基本上,每个给定时间月都有一个分辨率约为 1 度的全球地图。我想创建一个涵盖整个时期的 xarray 数据集。但是,一些 HDF5 文件是空的(例如,给定月份没有数据)。

最好的方法是什么?我认为这xr.open_mfdataset会起作用,但这不起作用(可能是因为 HDF5 有纬度/经度的“假”变量)。使用 rasterio 作为我想要的工作的子数据集:

def read_month(fname):
    ds = xr.open_rasterio(f'HDF5:"{fname}"://layer', chunks={"band":1})
    return ds
    
da = xr.concat([read_month(f) for f in fnames], 
               dim="band", coords={"band":times}
              )

fnames上面完全错过了空文件,并且代码以正确的结构读取,但读者没有解释地理参考或时间信息。

<xarray.DataArray (band: 19, y: 180, x: 360)>
dask.array<concatenate, shape=(19, 180, 360), dtype=float64, chunksize=(1, 180, 360), chunktype=numpy.ndarray>
Coordinates:
  * band     (band) int64 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  * y        (y) float64 0.5 1.5 2.5 3.5 4.5 ... 175.5 176.5 177.5 178.5 179.5
  * x        (x) float64 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
Attributes:
    transform:   (1.0, 0.0, 0.0, 0.0, 1.0, 0.0)
    res:         (1.0, -1.0)
    is_tiled:    0
    nodatavals:  (0.0,)
    scales:      (1.0,)
    offsets:     (0.0,)

更新

我安装了 hdf5netcdf,并尝试如下读取数据(基本上,不要尝试解码 CF 数据,并处理“假暗淡”)。然后我需要重命名虚假维度并翻转纬度:

fnames = sorted([f for f in Path("/somewhere/over/the/rainbow").glob("nightmare*h5") if f.stat().st_size > 10500])
dates = pd.to_datetime([f.name.split(".")[-2] for f in fnames], format="%Y%j")
for fname in fnames:    
    xds = xr.open_dataset(fname, chunks={},
                          engine="h5netcdf", 
                          backend_kwargs = {"phony_dims":"access"}, 
                          decode_cf=False
                         )

    xds = xds.swap_dims({"phony_dim_0":"lat", "phony_dim_1":"lon"}
                    ).assign_coords({"lat":np.linspace(90, -90, 180)})
    ss.append(xds)

da=xr.concat(ss, dates)
da = da.swap_dims({"concat_dim":"time"})

据我所知,这工作正常,但它很hacky并且需要很长时间才能完成。我想我现在可以将它保存为 netcdf 文件,但如果能够看看这是否是正确的方法,那就太好了。

标签: python-xarray

解决方案


推荐阅读