feat: OmGrid interface based on crs_wkt

examples/plot_map.py

@@ -3,9 +3,7 @@
 # /// script
 # requires-python = ">=3.12"
 # dependencies = [
-#     "omfiles==1.0.1",
-#     "fsspec>=2025.7.0",
-#     "s3fs",
+#     "omfiles[fsspec,grids]>=1.1.0",
 #     "matplotlib",
 #     "cartopy",
 # ]
@@ -17,69 +15,61 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from omfiles import OmFileReader
+from omfiles.grids import OmGrid
 
+MODEL_DOMAIN = "dmi_harmonie_arome_europe"
+VARIABLE = "relative_humidity_2m"
 # Example: URI for a spatial data file in the `data_spatial` S3 bucket
 # See data organization details: https://github.com/open-meteo/open-data?tab=readme-ov-file#data-organization
 # Note: Spatial data is only retained for 7 days. The example file below may no longer exist.
 # Please update the URI to match a currently available file.
-s3_uri = "s3://openmeteo/data_spatial/dwd_icon/2025/09/23/0000Z/2025-09-30T0000.om"
-
-# The following two incantations are equivalent
-#
-# from fsspec.implementations.cached import CachingFileSystem
-# from s3fs import S3FileSystem
-# s3_fs = S3FileSystem(anon=True, default_block_size=65536, default_cache_type="none")
-# backend = CachingFileSystem(
-#     fs=s3_fs, cache_check=3600, block_size=65536, cache_storage="cache", check_files=False, same_names=True
-# )
-# with OmFileReader.from_fsspec(backend, s3_uri) as reader:
+S3_URI = f"s3://openmeteo/data_spatial/{MODEL_DOMAIN}/2026/01/10/0000Z/2026-01-12T0000.om"
 
 backend = fsspec.open(
-    f"blockcache::{s3_uri}",
+    f"blockcache::{S3_URI}",
     mode="rb",
     s3={"anon": True, "default_block_size": 65536},
     blockcache={"cache_storage": "cache"},
 )
 with OmFileReader(backend) as reader:
     print("reader.is_group", reader.is_group)
 
-    child = reader.get_child_by_name("relative_humidity_2m")
+    child = reader.get_child_by_name(VARIABLE)
     print("child.name", child.name)
-
-    # Get the full data array
     print("child.shape", child.shape)
     print("child.chunks", child.chunks)
+    # Get the full data array
     data = child[:]
     print(f"Data shape: {data.shape}")
     print(f"Data range: {np.nanmin(data)} to {np.nanmax(data)}")
 
     # Create the plot
     fig = plt.figure(figsize=(12, 8))
-    ax = plt.axes(projection=ccrs.PlateCarree())  # use PlateCarree projection
+    ax = plt.axes(projection=ccrs.PlateCarree())
 
     # Add map features
-    ax.add_feature(cfeature.COASTLINE)
-    ax.add_feature(cfeature.BORDERS)
+    ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
+    ax.add_feature(cfeature.BORDERS, linewidth=0.5)
     ax.add_feature(cfeature.OCEAN, alpha=0.3)
     ax.add_feature(cfeature.LAND, alpha=0.3)
 
     # Create coordinate arrays
-    # Currently, the files don't contain any information about the spatial coordinates,
-    # so you need to provide these coordinate arrays manually.
-    height, width = data.shape
-    lon = np.linspace(-180, 180, width)  # Adjust these bounds
-    lat = np.linspace(-90, 90, height)  # Adjust these bounds
-    lon_grid, lat_grid = np.meshgrid(lon, lat)
+    num_y, num_x = data.shape
+    grid = OmGrid(reader.get_child_by_name("crs_wkt").read_scalar(), (num_y, num_x))
+    lon_grid, lat_grid = grid.get_meshgrid()
+    crs = grid.crs
+    assert crs is not None, "CRS is None, this should only happen for gaussian grids"
 
     # Plot the data
-    im = ax.contourf(lon_grid, lat_grid, data, levels=20, transform=ccrs.PlateCarree(), cmap="viridis")
-    plt.colorbar(im, ax=ax, shrink=0.6, label=child.name)
+    im = ax.contourf(lon_grid, lat_grid, data, cmap="coolwarm")
     ax.gridlines(draw_labels=True, alpha=0.3)
-    plt.title(f"2D Map: {child.name}")
-    ax.set_global()
+    ax.set_aspect("equal")
+    # ax.set_global()
+    plt.colorbar(im, ax=ax, orientation="vertical", pad=0.05, aspect=40, shrink=0.55, label=VARIABLE)
+    plt.title(f"{MODEL_DOMAIN} {VARIABLE} Map\nCRS: {crs.name}", fontsize=12, fontweight="bold", pad=16)
     plt.tight_layout()
 
-    output_filename = f"map_{child.name.replace('/', '_')}.png"
+    output_filename = f"map_{VARIABLE}.png"
     plt.savefig(output_filename, dpi=300, bbox_inches="tight")
     print(f"Plot saved as: {output_filename}")
     plt.close()

examples/readme_example.py

@@ -3,25 +3,25 @@
 # /// script
 # requires-python = ">=3.12"
 # dependencies = [
-#     "omfiles==1.0.1",
-#     "fsspec>=2025.7.0",
-#     "s3fs",
+#     "omfiles[fsspec]>=1.1.0",
 # ]
 # ///
 
 import fsspec
 import numpy as np
 from omfiles import OmFileReader
 
+MODEL_DOMAIN = "dwd_icon"
+VARIABLE = "temperature_2m"
 # Example: URI for a spatial data file in the `data_spatial` S3 bucket
 # See data organization details: https://github.com/open-meteo/open-data?tab=readme-ov-file#data-organization
 # Note: Spatial data is only retained for 7 days. The example file below may no longer exist.
 # Please update the URI to match a currently available file.
-s3_uri = "s3://openmeteo/data_spatial/dwd_icon/2025/09/23/0000Z/2025-09-30T0000.om"
+S3_URI = f"s3://openmeteo/data_spatial/{MODEL_DOMAIN}/2026/01/10/0000Z/2026-01-12T0000.om"
 
 # Create and open filesystem, wrapping it in a blockcache
 backend = fsspec.open(
-    f"blockcache::{s3_uri}",
+    f"blockcache::{S3_URI}",
     mode="rb",
     s3={"anon": True, "default_block_size": 65536},  # s3 settings
     blockcache={"cache_storage": "cache"},  # blockcache settings
@@ -35,7 +35,7 @@
     print(f"root.is_scalar: {root.is_scalar}")  # False
     print(f"root.is_group: {root.is_group}")  # True
 
-    temperature_reader = root.get_child_by_name("temperature_2m")
+    temperature_reader = root.get_child_by_name(VARIABLE)
     print(f"temperature_reader.is_array: {temperature_reader.is_array}")  # True
     print(f"temperature_reader.is_scalar: {temperature_reader.is_scalar}")  # False
     print(f"temperature_reader.is_group: {temperature_reader.is_group}")  # False

examples/regrid_ecmwf_ifs_hres_gaussian_O1280_to_0.1_degree.py

@@ -3,9 +3,7 @@
 # /// script
 # requires-python = ">=3.12"
 # dependencies = [
-#     "omfiles>=1.0.1",
-#     "fsspec>=2025.7.0",
-#     "s3fs",
+#     "omfiles[fsspec]>=1.1.0",
 #     "matplotlib",
 #     "cartopy",
 #     "earthkit-regrid==0.5.0",
@@ -20,22 +18,24 @@
 from earthkit.regrid import interpolate
 from omfiles import OmFileReader
 
+MODEL_DOMAIN = "ecmwf_ifs"
+VARIABLE = "temperature_2m"
 # Example: URI for a spatial data file in the `data_spatial` S3 bucket
 # See data organization details: https://github.com/open-meteo/open-data?tab=readme-ov-file#data-organization
 # Note: Spatial data is only retained for 7 days. The example file below may no longer exist.
 # Please update the URI to match a currently available file.
-s3_ifs_spatial_uri = f"s3://openmeteo/data_spatial/ecmwf_ifs/2025/10/01/0000Z/2025-10-01T0000.om"
+S3_URI = f"s3://openmeteo/data_spatial/{MODEL_DOMAIN}/2026/01/10/0000Z/2026-01-12T0000.om"
 
 backend = fsspec.open(
-    f"blockcache::{s3_ifs_spatial_uri}",
+    f"blockcache::{S3_URI}",
     mode="rb",
     s3={"anon": True, "default_block_size": 65536},
     blockcache={"cache_storage": "cache"},
 )
 with OmFileReader(backend) as reader:
     print("reader.is_group", reader.is_group)
 
-    child = reader.get_child_by_name("temperature_2m")
+    child = reader.get_child_by_name(VARIABLE)
     print("child.name", child.name)
 
     # Get the full data array
@@ -52,11 +52,11 @@
 
     # Create plot
     fig = plt.figure(figsize=(12, 8))
-    ax = plt.axes(projection=ccrs.PlateCarree())  # use PlateCarree projection
+    ax = plt.axes(projection=ccrs.PlateCarree())
 
     # Add map features
-    ax.add_feature(cfeature.COASTLINE)
-    ax.add_feature(cfeature.BORDERS)
+    ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
+    ax.add_feature(cfeature.BORDERS, linewidth=0.5)
     ax.add_feature(cfeature.OCEAN, alpha=0.3)
     ax.add_feature(cfeature.LAND, alpha=0.3)
 
@@ -68,14 +68,14 @@
     lon_grid, lat_grid = np.meshgrid(lon, lat)
 
     # Plot the data
-    im = ax.contourf(lon_grid, lat_grid, regridded, levels=20, transform=ccrs.PlateCarree(), cmap="viridis")
-    plt.colorbar(im, ax=ax, shrink=0.6, label=child.name)
-    ax.gridlines(draw_labels=True, alpha=0.3)
-    plt.title(f"2D Map: {child.name}")
+    im = ax.contourf(lon_grid, lat_grid, regridded, levels=20, cmap="coolwarm")
     ax.set_global()
+    ax.gridlines(draw_labels=True, alpha=0.3)
+    plt.colorbar(im, ax=ax, orientation="vertical", pad=0.05, aspect=40, shrink=0.55, label=VARIABLE)
+    plt.title(f"{MODEL_DOMAIN} {VARIABLE} Regridded to 0.1° Map", fontsize=12, fontweight="bold", pad=16)
     plt.tight_layout()
 
-    output_filename = f"map_ifs_{child.name.replace('/', '_')}.png"
+    output_filename = f"map_ifs_{VARIABLE}.png"
     plt.savefig(output_filename, dpi=300, bbox_inches="tight")
     print(f"Plot saved as: {output_filename}")
     plt.close()

examples/spatial_xarray.py

@@ -3,10 +3,7 @@
 # /// script
 # requires-python = ">=3.12"
 # dependencies = [
-#     "omfiles==1.0.1",
-#     "fsspec>=2025.7.0",
-#     "s3fs",
-#     "xarray",
+#     "omfiles[fsspec,grids,xarray]>=1.1.0",
 #     "matplotlib",
 #     "cartopy",
 # ]
@@ -18,63 +15,57 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import xarray as xr
-from omfiles import OmFileReader
+from omfiles.grids import OmGrid
 
-PLOT_VARIABLE = "temperature_2m"
-MODEL = "dwd_icon"
+MODEL_DOMAIN = "dwd_icon"
+VARIABLE = "temperature_2m"
 
 # Example: URI for a spatial data file in the `data_spatial` S3 bucket
 # See data organization details: https://github.com/open-meteo/open-data?tab=readme-ov-file#data-organization
 # Note: Spatial data is only retained for 7 days. The example file below may no longer exist.
 # Please update the URI to match a currently available file.
-s3_spatial_uri = f"s3://openmeteo/data_spatial/{MODEL}/2025/09/23/0000Z/2025-09-30T0000.om"
+S3_URI = f"s3://openmeteo/data_spatial/{MODEL_DOMAIN}/2026/01/10/0000Z/2026-01-12T0000.om"
 
 backend = fsspec.open(
-    f"blockcache::{s3_spatial_uri}",
+    f"blockcache::{S3_URI}",
     mode="rb",
     s3={"anon": True, "default_block_size": 65536},
     blockcache={"cache_storage": "cache"},
 )
 
 ds = xr.open_dataset(backend, engine="om")  # type: ignore
+print(ds.attrs)
 print(ds.variables.keys())  # any of these keys can be used for plotting
-ds = ds.rename_dims({"dim0": "lat", "dim1": "lon"})
-# You need to know what exactly the dimensions of the specified domain is referring to.
-# Icon is using a global regular grid:
-# https://github.com/open-meteo/open-meteo/blob/a4cdae1ad139f9dfa6dd2552c0636c7e572dcb52/Sources/App/Icon/Icon.swift#L146
-ds["lat"] = np.linspace(-90, 90, ds.sizes["lat"], endpoint=True)
-ds["lon"] = np.linspace(-180, 180, num=ds.sizes["lon"], endpoint=False)
 
 fig = plt.figure(figsize=(12, 8))
-ax = ax = plt.axes(projection=ccrs.PlateCarree())
+ax = plt.axes(projection=ccrs.PlateCarree())
 
 ax.add_feature(cfeature.COASTLINE, linewidth=0.8)
 ax.add_feature(cfeature.BORDERS, linewidth=0.5)
 ax.add_feature(cfeature.OCEAN, alpha=0.3)
 ax.add_feature(cfeature.LAND, alpha=0.3)
-ax.add_feature(cfeature.LAKES, alpha=0.3)
-ax.add_feature(cfeature.RIVERS, alpha=0.3)
 
-plot_data = ds[PLOT_VARIABLE]  # shape: (lat, lon)
-lon2d, lat2d = np.meshgrid(ds["lon"].values, ds["lat"].values)
+data = ds[VARIABLE]  # shape: (lat, lon)
+# Use OmGrid with the crs_wkt attribute to get the lat/lon grid
+grid = OmGrid(ds.attrs["crs_wkt"], ds[VARIABLE].shape)
+lon2d, lat2d = grid.get_meshgrid()
 
-min = int(plot_data.min().values)
-max = int(plot_data.max().values)
+min = int(data.min().values)
+max = int(data.max().values)
 stepsize = int((max - min) / 30)
 
-c = ax.contourf(
+im = ax.contourf(
     lon2d,
     lat2d,
-    plot_data,
+    data,
     levels=np.arange(min, max, stepsize),
-    cmap="Spectral_r",  # or "RdYlBu_r"
+    cmap="Spectral_r",
     vmin=min,
     vmax=max,
-    transform=ccrs.PlateCarree(),
     extend="both",
 )
-cb = plt.colorbar(c, ax=ax, orientation="vertical", pad=0.02, aspect=40, shrink=0.8)
-cb.set_label(PLOT_VARIABLE, fontsize=14)
-plt.title(f"{MODEL} {PLOT_VARIABLE}", fontsize=14, fontweight="bold", pad=20)
+ax.gridlines(draw_labels=True, alpha=0.3)
+plt.colorbar(im, ax=ax, orientation="vertical", pad=0.05, aspect=40, shrink=0.55, label=VARIABLE)
+plt.title(f"{MODEL_DOMAIN} {VARIABLE}", fontsize=12, fontweight="bold", pad=16)
 plt.tight_layout()
 plt.savefig("xarray_map.png", dpi=300, bbox_inches="tight")

pyproject.toml

@@ -21,7 +21,7 @@ maintainers = [
     { name = "Terraputix", email = "terraputix@mailbox.org" },
 ]
 dependencies = [
-    "numpy>=1.20.0",
+    "numpy>=1.21.0",
 ]
 
 [project.urls]
@@ -36,14 +36,16 @@ codec = [
     "numcodecs>=0.12.1",
 ]
 xarray = ["xarray>=2023.1.0"]
-fsspec = ["fsspec>=2023.10.0", "s3fs>=2023.1.0"]
+fsspec = ["fsspec>=2023.1.0", "s3fs>=2023.1.0"]
+grids = ["pyproj>=3.1.0"]
 all = [
     "zarr>=2.18.2",
     "numcodecs>=0.12.1",
     "xarray>=2023.1.0",
     "fsspec>=2023.10.0",
-    "s3fs>=2023.1.0"
-    ]
+    "s3fs>=2023.1.0",
+    "pyproj>=3.1.0"
+]
 
 [dependency-groups]
 documentation = [

python/omfiles/grids/__init__.py

@@ -0,0 +1,11 @@
+"""Provides classes for representing grids in OM files."""
+
+from .gaussian import GaussianGrid
+from .om_grid import OmGrid
+from .regular import RegularGrid
+
+__all__ = [
+    "GaussianGrid",
+    "OmGrid",
+    "RegularGrid",
+]