Using the Zarr library to read HDF5

[rsignell-usgs]

The USGS contracted the HDFGroup to do a test:
Could we make HDF5 format as performant on the Cloud as Zarr format by writing the HDF5 chunk locations into .zmetadata and then having the Zarr library read from those chunks instead of Zarr format chunks?

From our first test the answer appears to be YES: https://gist.github.com/rsignell-usgs/3cbe15670bc2be05980dec7c5947b540

We modified both the zarr and xarray libraries to make that notebook possible, adding the FileChunkStore concept. The modified libraries are: https://github.com/rsignell-usgs/hurricane-ike-water-levels/blob/zarr-hdf5/binder/environment.yml#L20-L21

Feel free to try running the notebook yourself:
(If you run into a 'stream is closed` error computing the max of the zarr data, just run the cell again.
I'm trying to figure out why that error occurs sometimes)

[Cadair]

I like the sound of the FileChunkStore. The ability to delegate reading chunks (and also metadata) to other libraries would be super useful for instance for asdf-format/asdf#718

[rabernat]

This recent blog post by @RPrudden gives some very pertinent suggestions:
https://medium.com/informatics-lab/arrays-on-the-fly-9509b6b2e46a

[alimanfoo]

Hi @rsignell-usgs, thanks a lot for posting, nice proof of concept.

[rsignell-usgs]

Thanks @alimanfoo! As @RPrudden pointed out to me, probably the neatest thing here is a mechanism for allowing existing community formats to be used efficiently via the Zarr library.

I've got a draft Medium post, in case folks are interested in commenting.

[djhoese]

Moving a question from gitter here: What would it take to do this type of operation for HDF5 files without generating the external .zmetadata file? For HDF5, could the C library be modified to read some byte range for the header, parse out the chunk locations for variable data, and then use that to load the chunks of data? I'm not sure if the HDF5 format has a static length header size which would probably be required to do this reliably.

[rsignell-usgs]

@ajelenak can provide a better answer, but my understanding is that it would be pretty complicated to make the HDF5 library have the functionality we demonstrated here. That's what made using the Zarr library very convenient!

I know there are several annoying steps in the current workflow that could be improved upon. We could imagine computing the augmented .zmetadata on the fly, or perhaps invent some JSON convention "virtual zarr file" that points to both (1) the augmented .zmetadata and (2) the binary file containing the chunks.

Do folks have other ideas how this could be made more user friendly?

[djhoese]

We could imagine computing the augmented .zmetadata on the fly

When the binary file is first created? Or when the binary data is being read?

I'm just hoping I can read GOES-R ABI NetCDF4 data from Amazon/Google without having to convince NOAA to add .zmetadata files for everything. 😉

[ajelenak]

The HDF5 library has the S3 Virtual File Driver (released in v1.10.6, I think) that enables access to HDF5 files in S3. You also can use h5py with the library without this virtual file driver. I have an example notebook how to do that. Both of these methods are not optimized and need to make frequent requests for small amounts of file content to figure out chunk file locations.

To answer @rsignell-usgs, the HDF5 library would need some form of the very similar information as in .zmetadata so a post-creation step to generate that information would also be needed. Of course, the library could be extended to generate the "consolidated file content" on every file close operation.

One solution to this problem now could be that someone sets up a Lambda function that will generate .zmetadata each time NOAA uploads a new GOES file in their S3 bucket.

[rsignell-usgs]

@djhoese , you don't have to convince NOAA to create the augmented .zmetadata! You can create them yourself and start working effectively with that data right now!

[djhoese]

you don't have to convince NOAA to create the augmented .zmetadata! You can create them yourself and start working effectively with that data right now!

I thought I needed the entire input file to properly generate the .zmetadata with the script mentioned in the medium post?

[rsignell-usgs]

@djhoese, you just need to extract the metadata from the existing GOES NetCDF4 files and stick it in another file. Then you reference both files when you read, as in our example notebook

[djhoese]

you just need to extract the metadata from the existing GOES NetCDF4 files

I just read the notebook @ajelenak linked to. This makes it more clear. When the Python file-like object from fsspec is passed to h5py.File it doesn't read the entire file, but knows to only parse specific byte ranges to get all the metadata it needs. Even though it makes a ton of requests, It won't download the entire file which is what I was worried about. So theoretically you don't need to make the .zmetadata file, you could generate that information on the fly from a h5py.File object, but for the best performance and least amount of HTTP requests (as @ajelenak pointed out) a .zmetadata file should be created before processing. Correct me if I'm wrong.

[ajelenak]

You are correct.

[rsignell-usgs]

@djhoese, I'm afraid if you are after efficient reading of the GOES data, we would have some more work to do. I downloaded a sample file and the chunks are tiny:

import xarray as xr
file = '/users/rsignell/downloads/OR_ABI-L2-SSTF-M6_G16_s20200412100059_e20200412159366_c20200412206173.nc'
ds = xr.open_dataset(file)
print(ds['SST'].encoding)

produces:

{'zlib': True,
 'shuffle': False,
 'complevel': 1,
 'fletcher32': False,
 'contiguous': False,
 'chunksizes': (226, 226),
 'source': 'c:\\users\\rsignell\\downloads\\OR_ABI-L2-SSTF-M6_G16_s20200412100059_e20200412159366_c20200412206173.nc',
 'original_shape': (5424, 5424),
 'dtype': dtype('int16'),
 '_Unsigned': 'true',
 '_FillValue': 65535,
 'scale_factor': 0.00244163,
 'add_offset': 180.0,
 'coordinates': 'retrieval_local_zenith_angle quantitative_local_zenith_angle retrieval_solar_zenith_angle t y x'}

so each chunk is only 216x216 Int16 values, so 2162162/1e6=0.1MB and that's before compression!

@ajelenak had some ideas that when we encounter tiny chunks we might create meta chunks to make the S3 byte-range requests and dask jobs bigger. But that would require more thought and effort...

[rsignell-usgs]

For those interested, the Medium Blog Post on this work.