Add IFS humidity data set

Author: treigerm

README.md

@@ -19,6 +19,7 @@ To download all the data used for the benchmark run the following commands:
 uv run python -m climatebenchpress.data_loader.datasets.esa_biomass_cci
 uv run python -m climatebenchpress.data_loader.datasets.cams
 uv run python -m climatebenchpress.data_loader.datasets.ifs_uncompressed
+uv run python -m climatebenchpress.data_loader.datasets.ifs_humidity
 uv run python -m climatebenchpress.data_loader.datasets.nextgems
 uv run python -m climatebenchpress.data_loader.datasets.cmip6.access_ta
 uv run python -m climatebenchpress.data_loader.datasets.cmip6.access_tos

pyproject.toml

@@ -8,15 +8,9 @@ dependencies = [
     "cf-xarray~=0.10.0",
     "cftime~=1.6.0",
     "dask>=2024.12.0,<2025.4",
-<<<<<<< HEAD
     "earthkit-regrid~=0.5.0",
     "fsspec>=2024.10.0,<2025.4",
     "gribscan~=0.0.14",
-=======
-    "earthkit-regrid>=0.5.0",
-    "fsspec>=2024.10.0,<2025.4",
-    "gribscan>=0.0.14",
->>>>>>> 868ed98 (Add IFS uncompressed data)
     "healpy~=1.18.0",
     # These versions need to be pinned to be compatible with the NextGEMS
     # catalog at https://data.nextgems-h2020.eu/online.yaml.

src/climatebenchpress/data_loader/datasets/all.py

@@ -4,5 +4,6 @@
 from .cmip6.all import *
 from .era5 import *
 from .esa_biomass_cci import *
+from .ifs_humidity import *
 from .ifs_uncompressed import *
 from .nextgems import *

src/climatebenchpress/data_loader/datasets/ifs_humidity.py

@@ -0,0 +1,74 @@
+__all__ = ["IFSHumidityDataset"]
+
+import argparse
+from pathlib import Path
+
+import xarray as xr
+
+from .. import (
+    monitor,
+    open_downloaded_canonicalized_dataset,
+    open_downloaded_tiny_canonicalized_dataset,
+)
+from .abc import Dataset
+from .ifs_uncompressed import load_hplp_data, regrid_to_regular
+
+BASE_URL = "https://object-store.os-api.cci1.ecmwf.int/esiwacebucket"
+
+
+class IFSHumidityDataset(Dataset):
+    """Dataset for the humidity field of the uncompressed IFS data.
+
+    Contains data from the [hplp](https://apps.ecmwf.int/ifs-experiments/rd/hplp/)
+    experiment from the Integrated Forecasting System (IFS) model. Crucially,
+    this dataset contains uncompressed 64-bit floating point data.
+    """
+
+    name = "ifs-humidity"
+
+    @staticmethod
+    def download(download_path: Path, progress: bool = True):
+        donefile = download_path / "download.done"
+        if donefile.exists():
+            return
+
+        ds = load_hplp_data(leveltype="ml", gridtype="reduced_gg", step=0)
+        ds = ds[["q"]]
+        ds_regridded = regrid_to_regular(
+            ds,
+            in_grid={"grid": "O400"},
+            out_grid={"grid": [0.25, 0.25]},
+        )
+        downloadfile = download_path / "ifs_humidity.zarr"
+        with monitor.progress_bar(progress):
+            ds_regridded.to_zarr(
+                downloadfile, mode="w", encoding=dict(), compute=False
+            ).compute()
+
+    @staticmethod
+    def open(download_path: Path) -> xr.Dataset:
+        ds = xr.open_dataset(download_path / "ifs_humidity.zarr")
+        ds = ds.isel(time=slice(0, 1)).chunk(-1)
+
+        # Needed to make the dataset CF-compliant.
+        ds.longitude.attrs["axis"] = "X"
+        ds.latitude.attrs["axis"] = "Y"
+        ds.level.attrs["axis"] = "Z"
+        ds.time.attrs["standard_name"] = "time"
+        return ds
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--basepath", type=Path, default=Path())
+    args = parser.parse_args()
+
+    ds = open_downloaded_canonicalized_dataset(
+        IFSHumidityDataset, basepath=args.basepath
+    )
+    open_downloaded_tiny_canonicalized_dataset(
+        IFSHumidityDataset, basepath=args.basepath
+    )
+
+    for v, da in ds.items():
+        print(f"- {v}: {da.dims}")

src/climatebenchpress/data_loader/datasets/ifs_uncompressed.py

@@ -126,33 +126,49 @@ def regrid_to_regular(ds, in_grid, out_grid):
     """
     out_data = {var: [] for var in ds.data_vars}
     for var in ds.data_vars:
+        var_has_level = "level" in ds[var].dims
         for time in ds.time:
-            r = earthkit.regrid.interpolate(
-                ds[var].sel(time=time).values,
-                in_grid=in_grid,
-                out_grid=out_grid,
-                method="linear",
-            )
-            out_data[var].append(r)
+            if var_has_level:
+                level_data = []
+                for level in ds[var].level:
+                    r = earthkit.regrid.interpolate(
+                        ds[var].sel(time=time, level=level).values,
+                        in_grid=in_grid,
+                        out_grid=out_grid,
+                        method="linear",
+                    )
+                    level_data.append(r)
+                out_data[var].append(level_data)
+            else:
+                r = earthkit.regrid.interpolate(
+                    ds[var].sel(time=time).values,
+                    in_grid=in_grid,
+                    out_grid=out_grid,
+                    method="linear",
+                )
+                out_data[var].append(r)
 
     dx = out_grid["grid"][0]
-    assert out_grid["grid"][0] == out_grid["grid"][1], (
-        "Only grids with equal latitude and longitude spacing are supported."
-    )
+    assert (
+        out_grid["grid"][0] == out_grid["grid"][1]
+    ), "Only grids with equal latitude and longitude spacing are supported."
     lats = np.linspace(90, -90, int(180 / dx) + 1)
     lons = np.linspace(0, 360 - dx, int(360 / dx))
     coords = {
         "time": ds.time,
         "latitude": lats,
         "longitude": lons,
     }
-    out_ds = xr.Dataset(
-        {
-            var: (("time", "latitude", "longitude"), out_data[var])
-            for var in ds.data_vars
-        },
-        coords=coords,
-    )
+
+    data_vars = {}
+    for var in ds.data_vars:
+        if "level" in ds[var].dims:
+            coords["level"] = ds[var].level
+            data_vars[var] = (("time", "level", "latitude", "longitude"), out_data[var])
+        else:
+            data_vars[var] = (("time", "latitude", "longitude"), out_data[var])
+
+    out_ds = xr.Dataset(data_vars, coords=coords)
     return out_ds