Grib / Kerchunk recipe?

[rabernat]

The new ECMWF forecast dataset in Azure is a perfect chance for us develop a close but not quite functional type of recipe: indexing and combing a bunch of grib files with Kerchunk.

The current HDFReferenceRecipe class almost does what we need; it just needs to be generalized to permit grib inputs. The kerchunk HRRR example shows how to open grib with kerchunk.

Here is some code to kick things off. cc @cisaacstern, @lsterzinger, @martindurant, @TomAugspurger

from kerchunk.grib2 import scan_grib
from fsspec.implementations.reference import ReferenceFileSystem
import zarr
import xarray as xr

ROOT = 'https://ai4edataeuwest.blob.core.windows.net/ecmwf/'
path = '20220122/18z/0p4-beta/scda/20220122180000-36h-scda-fc.grib2'
url = ROOT + path

storage_options = {}
common = ['time', 'latitude', 'longitude']
%time scan = scan_grib(url, common, storage_options)

rfs = ReferenceFileSystem(scan)
zg = zarr.open_group(rfs.get_mapper(""))
zg.info

%time _ = zg['q'][:]

import xarray as xr
ds = xr.open_zarr(rfs.get_mapper(""), consolidated=False)
print(ds)

<xarray.Dataset>
Dimensions:    (latitude: 451, longitude: 900)
Coordinates:
  * latitude   (latitude) float64 90.0 89.6 89.2 88.8 ... -89.2 -89.6 -90.0
  * longitude  (longitude) float64 -180.0 -179.6 -179.2 ... 178.8 179.2 179.6
Data variables: (12/18)
    d          (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    gh         (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    msl        (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    q          (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    r          (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    skt        (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    ...         ...
    u          (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    u10        (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    unknown    (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    v          (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    v10        (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
    vo         (latitude, longitude) float32 dask.array<chunksize=(451, 900), meta=np.ndarray>
Attributes:
    Conventions:             CF-1.7
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_edition:            2
    GRIB_subCentre:          0
    history:                 2022-01-26T04:08 GRIB to CDM+CF via cfgrib-0.9.9...
    institution:             European Centre for Medium-Range Weather Forecasts

[TomAugspurger]

Just an FYI, that storage container currently has a lifecycle policy to remove files older than 30-ish days. We're still trying to figure out how best to handle these operational forecast datasets. If anyone has thoughts on that, I'd be happy to hear them.

[martindurant]

I would just like to mention that GRIB2 support is not as good as hdf: it still needs cfgrib installed and works by copying pieces of the target file to local temporary files. It should be possible to eliminate that and get direct access to the binary chunks, and so better performance, but I have not had the time to try myself and was hoping someone would emerge with grib expertise who can help ( @DPeterK :) ).

[rabernat]

Perhaps we want to define a recipe that will create a kerchunk-virtual-zarr that always points to the "latest" data, whatever that means. Kind of similar to the ideas explored in this met office blog post. By running this recipe daily, we will ensure that we are never pointing at data that doesn't exist.

[rabernat]

@martindurant - the code above worked quite well out of the box! It took about 20s to index the file, and I didn't even run it from inside Azure. Whether or not a local copies are happening may not be so important here.

[martindurant]

Good to know!

[TomAugspurger]

@martindurant this ECMWF dataset includes "index files". Has anyone checked if those have sufficient information to translate into the index filesystem format needed for a reference filesystem?

https://ai4edataeuwest.blob.core.windows.net/ecmwf/20220126/00z/0p4-beta/enfo/20220126000000-0h-enfo-ef.index is one example, for the file at https://ai4edataeuwest.blob.core.windows.net/ecmwf/20220126/00z/0p4-beta/enfo/20220126000000-0h-enfo-ef.grib2

[rabernat]

Here is what the index file looks like

{"domain": "g", "date": "20220126", "time": "0000", "expver": "0001", "class": "od", "type": "pf", "stream": "enfo", "step": "0", "levtype": "sfc", "number": "29", "param": "2t", "_offset": 0, "_length": 609069}
{"domain": "g", "date": "20220126", "time": "0000", "expver": "0001", "class": "od", "type": "pf", "stream": "enfo", "step": "0", "levtype": "sfc", "number": "27", "param": "10u", "_offset": 609069, "_length": 609069}
{"domain": "g", "date": "20220126", "time": "0000", "expver": "0001", "class": "od", "type": "pf", "stream": "enfo", "step": "0", "levtype": "sfc", "number": "14", "param": "10v", "_offset": 1218138, "_length": 609069}
{"domain": "g", "date": "20220126", "time": "0000", "expver": "0001", "class": "od", "type": "pf", "stream": "enfo", "step": "0", "levtype": "sfc", "number": "2", "param": "2t", "_offset": 1827207, "_length": 609069}
{"domain": "g", "date": "20220126", "time": "0000", "expver": "0001", "class": "od", "type": "pf", "stream": "enfo", "step": "0", "levtype": "sfc", "number": "43", "param": "10u", "_offset": 2436276, "_length": 609069}
...

4182 lines of that.

It definitely looks like the information is there. Sounds like we need a kerchunk scan_grib_index function!

[martindurant]

If those "offsets" are what I think they are (will see), then this might indeed be the case that we don't need cfgrib. Otherwise, the JSON looks just a little cryptic.

[darothen]

This is fantastic... I tried this out with two of NOAA's operational datasets archived on GCS and it worked without any modifications, but was much slower (for calibration, the MWE here took ~40 seconds for scan_grib to run, unclear to me if this is a personal internet bandwidth issue or a GCS issue).

HRRR - https://console.cloud.google.com/storage/browser/high-resolution-rapid-refresh/ - took ~5 minutes for scan_grib to run
GFS - https://console.cloud.google.com/storage/browser/global-forecast-system/ - took ~7 minutes for scan_grib to run, but there was some odd bug in decoding the longitude coordinate and one of the values wound up as a NaN.

@TomAugspurger I was just about to call out the index files, these should be available for most NWP datasets distributed via GRIB (they are for NOAA). There's a good bit of information on leveraging systems that support reading via random access to subset byteranges corresponding to data fields in the GRIB using the offsets present in the index - see https://nomads.ncep.noaa.gov/txt_descriptions/fast_downloading_grib.shtml

[lsterzinger]

From the docs @TomAugspurger linked: "In addition, the keys _offset and _length represent the byte offset and length of the corresponding field. This allows the download of a single field using the HTTP Byte-Range request" so this definitely seems workable

[TomAugspurger]

Thanks all. I'll spend a bit of time this morning seeing if I can wire something up between the index files and fsspec reference filesystem.

[rabernat]

then this might indeed be the case that we don't need cfgrib

even with the index file, it seems like there must be some metadata in the grib file that is not covered by the index file.

[martindurant]

As far as I can tell, these "fields" will still need decoding. Maybe those other details in each JSON line tells you the dtype and compression.

[darothen]

Tag @blaylockbk - his package Herbie has some nascent capability to do partial decoding of HRRR grib files, could be a starting point here.

Somehow the grib_ls tool that ships with eccodes generates a much more detailed catalog of GRIB file contents, but I'm not sure how it works. That's another avenue to potentially explore.