Advice on how to append multiple dimensions when to_zarr

[songhan89]

Hi,

I have a ECMWF Ensemble GRIB dataset that I wish to save as .zarr format.

I can create a zarr store but found that when i try to append_dim (ingest new grib data everyday and want to append to the same zarr store), it is only possible to append to one dimension.

Ideally, i would want to set append_dim = ['time', 'valid_time']

The xarray dataset structure is as followed.

Would like to check with folks here what is your advice on this ?

• Songhan

[jhamman]

Hi @songhan89 - nice to see you around here.

In your example, are time and valid_time actually distinct dimensions? Can you also share what the data variable shapes you expect to have before and after the append step?

[songhan89]

Hi @jhamman ,

Sorry for the slow response. The use case is basically to compare different model forecast runs for the same time of the day (e.g comparing +24 hr from 00Z run vs +12 hr from 12Z run)

I experimented using earthkit.data to open the ECMWF grib data and turned out that it returned nicely in what I envisioned. Not sure what ecmwf has built in the preprocessing. In this case, I only need to append_dim = time.

valid_time (model run time + forecast step hour) actually need not to be another separate distinct dimension, since it can be composed from time (model run time) and step (forecast step hour)

with earthkit + xarray

import earthkit
import xarray as xr
import zarr
import adlfs
from glob import glob

files = glob("./samples/N1S*")
earth_ds = earthkit.data.from_source_lazily('file', files)

args = {'engine': 'cfgrib',
 'filter_by_keys': {'dataType': 'fc',
  'typeOfLevel': ['surface', 'isobaricInhPa']}}
  
ds1 = earth_ds.to_xarray(**args)
ds1

xarray

Previously was only using xarray to open the grib file and this would require me to do additional steps to expand_dims, and that i would have set valid_time to be a distinct dimension (which is not necessary).

def preprocess(ds):
    return ds.expand_dims(['time', 'step', 'valid_time'])

xr.open_mfdataset(
    files,
    preprocess=preprocess,
    parallel=True,
    **args
)

Sample files

[dcherian]

A few points:

1. time and step are the orthogonal dimensions here, valid_time = time + step . You should not need to add it to the expand_dims call.
2. append_dim only works along a single dimension.

IIUC you should only need to append along time. Here's how:

1. construct a schema dataset for the update with chunking that matches what's on disk. This doesn't need to contain real values and can be done using dask.array.ones or similar. The shape of this dataset should match the size of latitude ,longitude, step dimensions on disk. The size of the time dimension is 1 assuming you are ingesting a single forecast run. Make sure the time coordinate has the appropriate timestamp. https://github.com/earth-mover/forecast-datacube-demo/blob/4e0c7bce3fc49203453d74b83460d4cf30b9b919/modal_app.py#L170-L172
2. Append that schema to the dataset on-disk schematic_update.to_zarr(store, mode="a-", append_dim="time", compute=False). The compute=False is important, this will not write any of the data variables but simply setup the metadata correctly and write the new timestamp to disk. (https://github.com/earth-mover/forecast-datacube-demo/blob/4e0c7bce3fc49203453d74b83460d4cf30b9b919/modal_app.py#L170-L195)
3. Now ingest the available grib files for that run. For this read each grib in to Xarray and write to disk using the region kwarg. Note the convenient region="auto" functionality (https://github.com/earth-mover/forecast-datacube-demo/blob/4e0c7bce3fc49203453d74b83460d4cf30b9b919/modal_app.py#L188-L195). ±This step can be parallelized as long as you respect the chunking along step appropriately.

At a conceptual level, you first "expand" the dataset on-disk using an "append". Then you write real values to the expanded portion using the region kwarg.