Help with Parquet storage

[norlandrhagen]

Hi there Martin,

I'm trying to learn how to use the parquet storage in the hope of creating a section in the Pythia cookbook to demonstrate this. I'm going off of this test and the code block in the Kerchunk docs.

My initial attempt is below, but when I try to read the parquet I get a FileNotFoundError:

FileNotFoundError: [Errno 2] No such file or directory: '[/Users/nrhagen/Documents/carbonplan/pythia/kerchunk-cookbook/notebooks/foundations/combined.parq/TSLB/refs.0.parq](https://file+.vscode-resource.vscode-cdn.net/Users/nrhagen/Documents/carbonplan/pythia/kerchunk-cookbook/notebooks/foundations/combined.parq/TSLB/refs.0.parq)'

The only file in the combined.parq directory is ['.zmetadata']. So it seems like I'm not writing the combined reference correctly.

Also, I'm not clear on this line from the test LazyReferenceMapper.create(10, temp_dir, fs). Is the first arg of LazyReferenceMapper.create the length of input files or?

Thanks again for the help! It would be great to figure out how to use the parquet functionality.

from tempfile import TemporaryDirectory
import xarray as xr
from kerchunk.combine import MultiZarrToZarr
from kerchunk.hdf import SingleHdf5ToZarr
import os 
from kerchunk import hdf, combine, df
from fsspec.implementations.reference import LazyReferenceMapper, ReferenceFileSystem
import fsspec 

file_pattern = [
  's3://wrf-se-ak-ar5/ccsm/rcp85/daily/2060/WRFDS_2060-01-01.nc',
  's3://wrf-se-ak-ar5/ccsm/rcp85/daily/2060/WRFDS_2060-01-02.nc'
]

single_ref_sets = [hdf.SingleHdf5ToZarr(_).translate() for _ in file_pattern]

fs = fsspec.filesystem("file")
td = TemporaryDirectory()
temp_dir = td.name
temp_dir = str(temp_dir)
out = LazyReferenceMapper.create(10, temp_dir, fs)

mzz = MultiZarrToZarr(
    single_ref_sets,
    remote_protocol="memory",
    concat_dims=["Time"],identical_dims=['south_north', 'west_east', 'interp_levels', 'soil_layers_stag'],
   out=out,
).translate()

if not os.path.exists("combined.parq"):
    os.makedirs("combined.parq")
df.refs_to_dataframe(mzz, "combined.parq")

fs = ReferenceFileSystem(
    "combined.parq", lazy=True)
ds = xr.open_dataset(
    fs.get_mapper(), engine="zarr",
    backend_kwargs={"consolidated": False}
)

cc @rsignell-usgs if you've already figured this out!

[martindurant]

Following translate(), you need

out.flush()   # write all non-full reference batches
fs = ReferenceFileSystem(temp_dir,  remote_protocol="s3", target_protocol="file")
ds = xr.open_dataset(
    fs.get_mapper(), engine="zarr",
    backend_kwargs={"consolidated": False}
)

The value of 10 is how many references are in each component parquet file. The higher the number, the more efficient the storage and fewer files; the lower the number, the better random access. A value of 10000 is fine for reach use - it is 10 in testing to make sure we have multiple batches.

Not providing a remote_protocol exposes a bug, btw, but that's not important right now.

[dwest77a]

I have two problems, one of which I was able to bypass although it may result in issues further on.
Firstly the metadata in .zmetadata appears to be missing details required by whatever parser is in ReferenceFileSystem (namely the zarr_format property is missing in a lot of cases. My process is pretty much the same as above:

pq = 'combined.parq'
files = ['file1.nc','file2.nc'...]
single_ref_sets = []
for url in files:
    print('[INFO]', url)
    single_ref_sets.append(hdf.SingleHdf5ToZarr(url).translate())
    
out = LazyReferenceMapper.create(10, pq, fs = fsspec.filesystem("file"))

out_dict = combine.MultiZarrToZarr(
 single_ref_sets,
 remote_protocol="memory",
 concat_dims=["time"],
 out=out).translate()

out.flush()

df.refs_to_dataframe(out_dict, pq)

fs = fsspec.implementations.reference.ReferenceFileSystem(
    pq, remote_protocol='s3', target_protocol="file", lazy=True)
ds = xr.open_dataset(
    fs.get_mapper(), engine="zarr",
    backend_kwargs={"consolidated": False, "decode_times": False}
)

The second issue is after I've 'fixed' the metadata by basically scraping all the .zarray/.zattrs from the kerchunk saved file for this set of data and dumping them into .zmetadata. This adds all the correct properties. But then I get different errors, always with the last steps of opening the parquet file. At the moment I get
GroupNotFoundError: group not found at path ''
which I think is because the kwarg 'group' is not provided in xr.open_dataset call?

[dwest77a]

I also get TypeError: argument of type 'float' is not iterable if I manage to bypass the group error, at the point of plotting the dataset because there's a chunk reference for time that ends up as ['nan'] instead of the [path, size, offset] structure, I think possibly because the time reference gets base64 encoded at the MultiZarrToZarr stage.

[martindurant]

Since you have an fs instance, the easiest way to debug would be to poke into this filesystem to see what's there, e.g. fs.ls(""). Is there an .zmetadata or .zgroup? Can you get their contents?

 remote_protocol="memory",

is this right?

[dwest77a]

Yes I can get the .zmetadata contents, that's how I know the metadata needs correcting as a lot of the .zarray parts are missing zarr_format.
I've worked out the time reference is base64 encoded which I think is causing the later issue, but if I just try this with inline_threshold added I get a different problem.

pq = 'combined.parq'
files = ['file1.nc','file2.nc'...]
single_ref_sets = []
for url in files:
    print('[INFO]', url)
    single_ref_sets.append(hdf.SingleHdf5ToZarr(url, inline_threshold=-1).translate())
    
out = LazyReferenceMapper.create(10, pq, fs = fsspec.filesystem("file"))

out_dict = combine.MultiZarrToZarr(
 single_ref_sets,
 remote_protocol="memory",
 concat_dims=["time"],
 out=out).translate()

I just get a keyerror in fsspec/implementations/memory:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/.local/lib/python3.10/site-packages/fsspec/implementations/memory.py:216, in MemoryFileSystem.cat_file(self, path, start, end, **kwargs)
    215 try:
--> 216     return bytes(self.store[path].getbuffer()[start:end])
    217 except KeyError:

KeyError: '/neodc/esacci/biomass/data/agb/maps/v4.0/netcdf/ESACCI-BIOMASS-L4-AGB-MERGED-100m-2010-fv4.0.nc'

The filepath definitely exists

[martindurant]

Note that this is in memory.py:

File ~/.local/lib/python3.10/site-packages/fsspec/implementations/memory.py:216, in MemoryFileSystem.cat_file(self, path, start, end, **kwargs)
    215 try:
--> 216     return bytes(self.store[path].getbuffer()[start:end])
    217 except KeyError:

because you specified remote_protocol="memory". Should be "file"?

[dwest77a]

Ok I've changed the remote_protocol to 'file' but still getting one of the references as [nan] instead of the normal format which is giving issues. I've taken a closer look and found it's because there are missing chunks in the variables I'm looking at, but when the parquet file is read in, the chunk reference is presumably set to 'nan'. I'll try and find out where this is happening.

Does parquet use any kind of templating to generate kerchunk references [path, size, offset] that it might just fill with [nan] when reconstructing from the metadata and parquet refs?

[martindurant]

Are there really missing chunks? It should come back with path=None (not nan and certainly not "nan"), size=offset=0.

[dwest77a]

I'll dig into the LazyReferenceMapper to find out where this is happening, but in the normal kerchunk json file for this dataset there are missing chunks where there is a land/sea mask in place

[martindurant]

hanks for digging. Probably me tests have defined/existing chunks everywhere.

In LazyReferenceMapper._load_one_key we have

        if selection[0] is None:
            raise KeyError("This reference has been deleted")

to catch missing keys, but apparently there is a nan here. It would be good to find out the dtype of the columns in refs and whether this is a np.nan, pd.NaN or "nan".

[martindurant]

Any luck making progress here?

[NikosAlexandris]

I got

..
│   200 │   │   # check metadata format version                                                    │
│   201 │   │   zarr_format = meta.get("zarr_format", None)                                        │
│   202 │   │   if zarr_format != cls.ZARR_FORMAT:                                                 │
│ ❱ 203 │   │   │   raise MetadataError("unsupported zarr format: %s" % zarr_format)               │
│   204 │   │                                                                                      │
│   205 │   │   meta = dict(zarr_format=zarr_format)                                               │
│   206 │   │   return meta                                                                        │
│                                                                                                  │
│ ╭─────────────────────────────────────────── locals ───────────────────────────────────────────╮ │
│ │         cls = <class 'zarr.meta.Metadata2'>                                                  │ │
│ │        meta = {                                                                              │ │
│ │               │   '_ARRAY_DIMENSIONS': ['time', 'lat', 'lon'],                               │ │
│ │               │   'add_offset': 0.0,                                                         │ │
│ │               │   'cell_methods': 'time: sum',                                               │ │
│ │               │   'long_name': 'Sunshine Duration',                                          │ │
│ │               │   'missing_value': -999,                                                     │ │
│ │               │   'scale_factor': 0.0010000000474974513,                                     │ │
│ │               │   'standard_name': 'duration_of_sunshine',                                   │ │
│ │               │   'units': 'h'                                                               │ │
│ │               }                                                                              │ │
│ │           s = b'{"_ARRAY_DIMENSIONS":["time","lat","lon"],"add_offset":0.0,"cell_methods":"… │ │
│ │ zarr_format = None                                                                           │ │
│ ╰──────────────────────────────────────────────────────────────────────────────────────────────╯ │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
MetadataError: unsupported zarr format: None

using

def create_kerchunk_reference(
    source_directory: Annotated[Path, typer_argument_source_directory],
    output_directory: Annotated[Path, typer_argument_output_directory],
    pattern: Annotated[str, typer_option_filename_pattern] = '*.nc',
    workers: Annotated[int, typer_option_number_of_workers] = 4,
    dry_run: Annotated[bool, typer.Option("--dry-run", help="Run the command without making any changes.")] = False,
    verbose: Annotated[int, typer_option_verbose] = VERBOSE_LEVEL_DEFAULT,
):
    """Reference local NetCDF files using Kerchunk"""
    file_paths = list(source_directory.glob(pattern))
    if not file_paths:
        logger.info(
            "No files found in the source directory matching the pattern."
        )
        return
    if dry_run:
        print(f"[bold]Dry run[/bold] of [bold]operations that would be performed[/bold]:")
        print(f"> Reading files in [code]{source_directory}[/code] matching the pattern [code]{pattern}[/code]")
        print(f"> Number of files matched: {len(file_paths)}")
        print(f"> Creating single reference files to [code]{output_directory}[/code]")
        return  # Exit for a dry run
    output_directory.mkdir(parents=True, exist_ok=True)

    # Map verbosity level to display mode
    mode = DisplayMode(verbose)
    with display_context[mode]:
        with multiprocessing.Pool(processes=workers) as pool:
            from functools import partial

            partial_create_single_reference = partial(
                create_single_reference, output_directory=output_directory
            )
            results = pool.map(partial_create_single_reference, file_paths)

and

def combine_kerchunk_references_to_parquet(
    source_directory: Annotated[Path, typer_argument_source_directory],
    pattern: Annotated[str, typer_option_filename_pattern] = "*.json",
    combined_reference: Annotated[Path, typer_argument_kerchunk_combined_reference] = "combined_kerchunk.parq",
    dry_run: Annotated[bool, typer.Option("--dry-run", help="Run the command without making any changes.")] = False,
    verbose: Annotated[int, typer_option_verbose] = VERBOSE_LEVEL_DEFAULT,
):
    """Combine multiple kerchunked datasets into a single logical parquet aggregate
    dataset using fsspec.implementations.reference.ReferenceFileSystem"""
    mode = DisplayMode(verbose)
    with display_context[mode]:
        source_directory = Path(source_directory)
        reference_file_paths = list(source_directory.glob(pattern))
        reference_file_paths = list(map(str, reference_file_paths))

        if dry_run:
            print(f"[bold]Dry run[/bold] of [bold]operations that would be performed[/bold]:")
            print(f"> Reading files in [code]{source_directory}[/code] matching the pattern [code]{pattern}[/code]")
            print(f"> Number of files matched: {len(reference_file_paths)}")
            print(f"> Writing combined reference file to [code]{combined_reference}[/code]")
            return  # Exit for a dry run

        # Create LazyReferenceMapper to pass to MultiZarrToZarr
        fs = fsspec.filesystem("file")
        import os

        combined_reference.mkdir(parents=True, exist_ok=True)
        from fsspec.implementations.reference import LazyReferenceMapper
        output_lazy = LazyReferenceMapper.create(1000, str(combined_reference), fs)

        from kerchunk.combine import MultiZarrToZarr

        # Combine single references
        mzz = MultiZarrToZarr(
            reference_file_paths,
            remote_protocol="file",
            concat_dims=["time"],
            identical_dims=["lat", "lon"],
            out=output_lazy,
        )
        multifile_kerchunk = mzz.translate()

        output_lazy.flush()  # Write all non-full reference batches

        # Read from the Parquet storage
        kerchunk.df.refs_to_dataframe(multifile_kerchunk, str(combined_reference))

        filesystem = fsspec.implementations.reference.ReferenceFileSystem(
            fo=str(combined_reference),
            remote_protocol='file',
            target_protocol="file",
            lazy=True
        )

        # test...
        ds = xr.open_dataset(
            filesystem.get_mapper(),
            engine="zarr",
            backend_kwargs={"consolidated": False},
        )
        print(ds)

[NikosAlexandris]

If I install

pip install \
git+https://github.com/fsspec/kerchunk \
git+https://github.com/martindurant/filesystem_spec.git@ref_mode \
fastparquet \
--force

and remove the .create() from output_lazy = LazyReferenceMapper.create(1000, str(combined_reference), fs)

I get

│   181 │   def _get_chunk_sizes(self, field):                                                     │
│   182 │   │   if field not in self.chunk_sizes:                                                  │
│   183 │   │   │   zarray = self.zmetadata[f"{field}/.zarray"]                                    │
│ ❱ 184 │   │   │   size_ratio = self.np.array(zarray["shape"]) / self.np.array(                   │
│   185 │   │   │   │   zarray["chunks"]                                                           │
│   186 │   │   │   )                                                                              │
│   187 │   │   │   self.chunk_sizes[field] = self.np.ceil(size_ratio).astype(int)                 │
│                                                                                                  │
│ ╭───────────────────────────────────────── locals ─────────────────────────────────────────╮     │
│ │  field = 'lon'                                                                           │     │
│ │   self = <fsspec.implementations.reference.LazyReferenceMapper object at 0x7f630eb73ad0> │     │
│ │ zarray = {'zarr_format': 2}                                                              │     │
│ ╰──────────────────────────────────────────────────────────────────────────────────────────╯     │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
KeyError: 'shape'

[martindurant]

@NikosAlexandris - can you make this a reproducible piece of code I can run? Are the files you are scanning public?

[NikosAlexandris]

@NikosAlexandris - can you make this a reproducible piece of code I can run? Are the files you are scanning public?

@martindurant The data are public, but you need to order + download. Update : Here 2 files. If you'd like more, where can I upload/share some (2 files ~350MB or 9 files ~1.6GB or a month of data ~5.4GB) ?

Following :

kerchunk_to_parquet.py works after

pip install kerchunk fsspec fastparquet h5py xarray typer rich 

for example

python kerchunk_to_parquet.py SIS/ references.json --source-pattern 'SISin2020010*.nc' 

and ends with

│ /software/kerchunk_to_parquet/.rekx_virtual_environment/lib/python3.11/site-packages/zarr/ │
│ meta.py:203 in decode_group_metadata                                                       │
│                                                                                            │
│   200 │   │   # check metadata format version                                              │
│   201 │   │   zarr_format = meta.get("zarr_format", None)                                  │
│   202 │   │   if zarr_format != cls.ZARR_FORMAT:                                           │
│ ❱ 203 │   │   │   raise MetadataError("unsupported zarr format: %s" % zarr_format)         │
│   204 │   │                                                                                │
│   205 │   │   meta = dict(zarr_format=zarr_format)                                         │
│   206 │   │   return meta                                                                  │
│                                                                                            │
│ ╭──────────────────────────────────────── locals ────────────────────────────────────────╮ │
│ │         cls = <class 'zarr.meta.Metadata2'>                                            │ │
│ │        meta = {                                                                        │ │
│ │               │   '_ARRAY_DIMENSIONS': ['lat'],                                        │ │
│ │               │   'axis': 'Y',                                                         │ │
│ │               │   'bounds': 'lat_bnds',                                                │ │
│ │               │   'long_name': 'latitude',                                             │ │
│ │               │   'standard_name': 'latitude',                                         │ │
│ │               │   'units': 'degrees_north'                                             │ │
│ │               }                                                                        │ │
│ │           s = b'{"_ARRAY_DIMENSIONS":["lat"],"axis":"Y","bounds":"lat_bnds","long_nam… │ │
│ │ zarr_format = None                                                                     │ │
│ ╰────────────────────────────────────────────────────────────────────────────────────────╯ │
╰────────────────────────────────────────────────────────────────────────────────────────────╯
MetadataError: unsupported zarr format: None

kerchunk_to_parquet_alternative.py works after

pip install git+https://github.com/fsspec/kerchunk git+https://github.com/martindurant/filesystem_spec.git@ref_mode fastparquet h5py xarray typer rich

and

python kerchunk_to_parquet_alternative.py SIS/ references.json --source-pattern 'SISin2020010*.nc'

ends with

│ /software/kerchunk_to_parquet/.rekx_virtual_environment/lib/python3.11/site-packages/fsspe │
│ c/implementations/reference.py:184 in _get_chunk_sizes                                     │
│                                                                                            │
│   181 │   def _get_chunk_sizes(self, field):                                               │
│   182 │   │   if field not in self.chunk_sizes:                                            │
│   183 │   │   │   zarray = self.zmetadata[f"{field}/.zarray"]                              │
│ ❱ 184 │   │   │   size_ratio = self.np.array(zarray["shape"]) / self.np.array(             │
│   185 │   │   │   │   zarray["chunks"]                                                     │
│   186 │   │   │   )                                                                        │
│   187 │   │   │   self.chunk_sizes[field] = self.np.ceil(size_ratio).astype(int)           │
│                                                                                            │
│ ╭──────────────────────────────────────── locals ────────────────────────────────────────╮ │
│ │  field = 'time'                                                                        │ │
│ │   self = <fsspec.implementations.reference.LazyReferenceMapper object at               │ │
│ │          0x7f1b28ffa750>                                                               │ │
│ │ zarray = {                                                                             │ │
│ │          │   '_ARRAY_DIMENSIONS': ['lon'],                                             │ │
│ │          │   'axis': 'X',                                                              │ │
│ │          │   'bounds': 'lon_bnds',                                                     │ │
│ │          │   'long_name': 'longitude',                                                 │ │
│ │          │   'standard_name': 'longitude',                                             │ │
│ │          │   'units': 'degrees_east'                                                   │ │
│ │          }                                                                             │ │
│ ╰────────────────────────────────────────────────────────────────────────────────────────╯ │
╰────────────────────────────────────────────────────────────────────────────────────────────╯
KeyError: 'shape'

ps- One can also do python kerchunk_to_parquet.py --help for usage details.