E4DRR/gik-gefs-par

Without GIK:  Download ~24 GB of GRIB files per day  →  then process
With GIK:     Read ~10 MB of parquet refs + stream    →  direct analysis

catalog.parquet                                    # 791 KB, 57,780-row index
run_par_gefs_agg/
  monthly_agg/
    {YYYY}/{MM}_00z.parquet                        # one optimised parquet per month

import pandas as pd

# Read only the rows for one date + one member — pyarrow skips non-matching row groups
df = pd.read_parquet(
    "hf://datasets/E4DRR/gik-gefs-par/run_par_gefs_agg/monthly_agg/2024/06_00z.parquet",
    filters=[("date", "=", "20240615"), ("member", "=", "gep01")],
)
print(f"{len(df)} reference rows for 2024-06-15 gep01")
print(df[["key", "member", "date"]].head())

df = pd.read_parquet(
    "hf://datasets/E4DRR/gik-gefs-par/run_par_gefs_agg/monthly_agg/2024/06_00z.parquet",
    filters=[("date", "=", "20240615")],
)
print(f"{df.member.nunique()} members, {len(df)} rows total")

import ast, base64
import dask, dask.array as da
import fsspec, gribberish
import numpy as np, pandas as pd, xarray as xr

GEFS_GRID = (721, 1440)
GEFS_LATS = np.linspace(90, -90, 721)
GEFS_LONS = np.linspace(0, 359.75, 1440)            # GEFS uses 0..360 longitude
STEPS = [3, 6, 12, 24, 48, 72, 120, 168, 240]       # forecast hours

# 1. Read all 30 members for one date in one HF call
#    (~30 MB pulled via parquet filter pushdown + HF range reads)
df = pd.read_parquet(
    "hf://datasets/E4DRR/gik-gefs-par/run_par_gefs_agg/monthly_agg/2024/06_00z.parquet",
    filters=[("date", "=", "20240615")],
)

def member_zstore(sub):
    blob = sub[sub["key"] == "refs"].iloc[0]["value"]
    if isinstance(blob, bytes):
        blob = blob.decode("utf-8")
    return ast.literal_eval(blob)

# 2. Lazy chunk fetcher — runs only when xarray pulls data
s3 = fsspec.filesystem("s3", anon=True)

def lazy_chunk(ref):
    @dask.delayed
    def _fetch():
        url, off, ln = ref[0], ref[1], ref[2]
        with s3.open(url, "rb") as f:
            f.seek(off)
            raw = f.read(ln)
        return gribberish.parse_grib_array(raw, 0).reshape(GEFS_GRID).astype(np.float32)
    return da.from_delayed(_fetch(), shape=GEFS_GRID, dtype=np.float32)

# 3. Build the (member, step) lazy stack — no S3 reads here, just metadata
member_arrays, member_names = [], []
for member, sub in df.groupby("member"):
    zs = member_zstore(sub)
    val = zs.get("tp/accum/surface/step/0", "")
    if isinstance(val, str) and val.startswith("base64:"):
        step_hours = np.frombuffer(base64.b64decode(val[7:]), dtype="<f8")
    else:
        step_hours = np.arange(0, 243, 3, dtype=float)
    tp_refs = {
        int(k.rsplit("/", 1)[1].split(".")[0]): v
        for k, v in zs.items()
        if k.startswith("tp/accum/surface/tp/") and isinstance(v, list)
    }
    chunks = []
    for h in STEPS:
        idx = int(np.argmin(np.abs(step_hours - h)))
        ref = tp_refs.get(idx)
        chunks.append(lazy_chunk(ref) if ref else
                      da.full(GEFS_GRID, np.nan, dtype=np.float32))
    member_arrays.append(da.stack(chunks, axis=0))
    member_names.append(member)

ds = xr.Dataset(
    {"tp": (["member", "step", "latitude", "longitude"],
            da.stack(member_arrays, axis=0))},
    coords={"member": member_names, "step": STEPS,
            "latitude": GEFS_LATS, "longitude": GEFS_LONS},
)
print(ds)
# <xarray.Dataset>  Dimensions: member=30, step=9, latitude=721, longitude=1440
# tp: dask.array<...>  ← zero bytes in memory, all delayed

# 4. Fetch only what you need — each .load() triggers parallel S3 byte-range reads
step48 = ds.tp.sel(step=48).load()                       # all 30 members, T+48h
ea = ds.tp.sel(step=48,
               latitude=slice(15, -12),                  # GEFS lat is 90..-90
               longitude=slice(25, 52)).load()           # East Africa subset (~5 MB total)