🚀 Introduce pyiem.grid.nav with a fancy interface

CHANGELOG.md

@@ -17,6 +17,7 @@ All notable changes to this library are documented in this file.
 ### New Features
 
 - Add `allowed_as_list` option to `iemapp()` helper to stop lists.
+- Add `pyiem.grid.nav` with IEM grid information in a fancy form.
 - Add `MapPlot.imshow` with some optimized panel plotting.
 - Add maximum risk threshold within SPC outlook message (#969).
 - Add `pyiem.era5land` with IEM grid reference information.

MANIFEST.in

@@ -1,3 +1,4 @@
 recursive-include pyiem *.py
+recursive-include pyiem *.pyi
 recursive-include pyiem *.txt
 recursive-include pyiem *.npy

src/pyiem/grid/nav.py

@@ -0,0 +1,102 @@
+"""Build Grid Navigation Metadata from the NetCDF file."""
+
+from pyiem.models.gridnav import CartesianGridNavigation
+
+_GRID_CONFIGS = {
+    "IEMRE": {
+        "left_edge": -126.0625,
+        "bottom_edge": 22.9375,
+        "dx": 0.125,
+        "dy": 0.125,
+        "nx": 488,
+        "ny": 216,
+    },
+    "IEMRE_CHINA": {
+        "left_edge": 69.9375,
+        "bottom_edge": 14.9375,
+        "dx": 0.125,
+        "dy": 0.125,
+        "nx": 560,
+        "ny": 320,
+    },
+    "IEMRE_EUROPE": {
+        "left_edge": -10.0625,
+        "bottom_edge": 34.9375,
+        "dx": 0.125,
+        "dy": 0.125,
+        "nx": 400,
+        "ny": 280,
+    },
+    # Lamely hardcoded for now
+    "ERA5LAND": {
+        "left_edge": -126.05,
+        "bottom_edge": 22.95,
+        "dx": 0.1,
+        "dy": 0.1,
+        "nx": 610,
+        "ny": 270,
+    },
+    "ERA5LAND_CHINA": {
+        "left_edge": 69.95,
+        "bottom_edge": 14.95,
+        "dx": 0.1,
+        "dy": 0.1,
+        "nx": 700,
+        "ny": 400,
+    },
+    "ERA5LAND_EUROPE": {
+        "left_edge": -10.05,
+        "bottom_edge": 34.95,
+        "dx": 0.1,
+        "dy": 0.1,
+        "nx": 500,
+        "ny": 350,
+    },
+    "STAGE4": {
+        "crs": (
+            "+proj=stere +a=6371200 +b=6371200 +lat_0=90 "
+            "+lon_0=-105, +lat_ts=60"
+        ),
+        "left_edge": -1_904_912.924,
+        "bottom_edge": -7_619_986.180,
+        "dx": 4_762.5,
+        "dy": 4_762.5,
+        "nx": 1121,
+        "ny": 881,
+    },
+    "STAGE4_PRE2002": {
+        "crs": (
+            "+proj=stere +a=6371200 +b=6371200 +lat_0=90 "
+            "+lon_0=-105, +lat_ts=60"
+        ),
+        "left_edge": -2_097_827.439,
+        "bottom_edge": -7_622_315.608,
+        "dx": 4_763.0,
+        "dy": 4_763.0,
+        "nx": 1160,
+        "ny": 880,
+    },
+    "MRMS_IEMRE": {  # Specific to the IEM and not in general
+        "left_edge": -126.0,
+        "bottom_edge": 23.0,
+        "dx": 0.01,
+        "dy": 0.01,
+        "nx": 6100,
+        "ny": 2700,
+    },
+    "PRISM": {
+        "left_edge": -125.0 - (1 / 24.0) / 2.0,
+        "bottom_edge": 24.083333 - (1 / 24.0) / 2.0,
+        "dx": 1 / 24.0,
+        "dy": 1 / 24.0,
+        "nx": 1405,
+        "ny": 621,
+    },
+}
+
+
+def __getattr__(name: str):
+    """Build stuff on the fly."""
+    if name in _GRID_CONFIGS:
+        return CartesianGridNavigation(**_GRID_CONFIGS[name])
+    raise AttributeError(f"module '{__name__}' has no attribute '{name}'")

src/pyiem/grid/nav.pyi

@@ -0,0 +1,14 @@
+"""Stub file."""
+
+from pyiem.models.gridnav import CartesianGridNavigation
+
+class IEMRE(CartesianGridNavigation): ...
+class IEMRE_CHINA(CartesianGridNavigation): ...
+class IEMRE_EUROPE(CartesianGridNavigation): ...
+class ERA5LAND(CartesianGridNavigation): ...
+class ERA5LAND_CHINA(CartesianGridNavigation): ...
+class ERA5LAND_EUROPE(CartesianGridNavigation): ...
+class STAGE4(CartesianGridNavigation): ...
+class STAGE4_PRE2002(CartesianGridNavigation): ...
+class MRMS_IEMRE(CartesianGridNavigation): ...
+class PRISM(CartesianGridNavigation): ...

src/pyiem/models/gridnav.py

@@ -0,0 +1,162 @@
+"""Grid Navigation Metadata."""
+
+from typing import Optional, Union
+
+import numpy as np
+from pydantic import BaseModel, ConfigDict, Field, model_validator
+from pyproj import CRS, Proj
+from rasterio.transform import from_origin
+
+
+class CartesianGridNavigation(BaseModel):
+    """Navigation for cartesian grid with (0,0) in lower left.
+
+    The `left_edge` and `bottom_edge` are the only required fields. The
+    rest are optional, but you need to have enough information to define
+    the grid, ie provide `dx` and `dy` or `nx` and `ny`.
+    """
+
+    model_config = ConfigDict(arbitrary_types_allowed=True)
+
+    crs: Union[str, CRS] = Field(
+        default="EPSG:4326",
+        description="The coordinate reference system of the grid",
+    )
+    left_edge: float = Field(
+        default=...,
+        description="The left edge of the grid in projection units",
+    )
+    bottom_edge: float = Field(
+        default=...,
+        description="The bottom edge of the grid in projection units",
+    )
+    top_edge: float = Field(
+        default=None,
+        description="The top edge of the grid in projection units",
+    )
+    right_edge: float = Field(
+        default=None,
+        description="The right edge of the grid in projection units",
+    )
+    dx: float = Field(
+        default=None,
+        description="The grid cell width in projection units",
+        gt=0,
+    )
+    dy: float = Field(
+        default=None,
+        description="The grid cell height in projection units",
+        gt=0,
+    )
+    nx: int = Field(
+        default=None,
+        description="The number of grid cells in the x direction",
+        gt=0,
+    )
+    ny: int = Field(
+        default=None,
+        description="The number of grid cells in the y direction",
+        gt=0,
+    )
+
+    @property
+    def x_points(self) -> np.ndarray:
+        """These are the centers of the cells in the x direction."""
+        return np.arange(self.nx) * self.dx + self.left
+
+    @property
+    def y_points(self) -> np.ndarray:
+        """These are the centers of the cells in the y direction."""
+        return np.arange(self.ny) * self.dy + self.bottom
+
+    @property
+    def x_edges(self) -> np.ndarray:
+        """These are the edges of the x cells (n=NX + 1)."""
+        return np.arange(self.nx + 1) * self.dx + self.left_edge
+
+    @property
+    def y_edges(self) -> np.ndarray:
+        """These are the edges of the y cells (n=NY + 1)."""
+        return np.arange(self.ny + 1) * self.dy + self.bottom_edge
+
+    @property
+    def left(self) -> float:
+        """The centroid of the left most grid cell."""
+        return self.left_edge + (self.dx / 2.0)
+
+    @property
+    def right(self) -> float:
+        """The centroid of the right most grid cell."""
+        return self.right_edge - (self.dx / 2.0)
+
+    @property
+    def bottom(self) -> float:
+        """The centroid of the bottom most grid cell."""
+        return self.bottom_edge + (self.dy / 2.0)
+
+    @property
+    def top(self) -> float:
+        """The centroid of the top most grid cell."""
+        return self.top_edge - (self.dy / 2.0)
+
+    @property
+    def affine(self):
+        """Return the affine transformation."""
+        return from_origin(self.left_edge, self.bottom_edge, self.dx, self.dy)
+
+    @property
+    def affine_image(self):
+        """Return the transformation associated with upper left origin."""
+        return from_origin(self.left_edge, self.top_edge, self.dx, -self.dy)
+
+    @model_validator(mode="before")
+    def complete_definition(cls, values):
+        """Use information that was provided to compute other fields."""
+        # We have required fields left_edge, bottom_edge
+        # Require that either dx/dy is provided or nx/ny is provided
+        if values.get("top_edge") is None:
+            values["top_edge"] = values["bottom_edge"] + (
+                values["ny"] * values["dy"]
+            )
+        if values.get("right_edge") is None:
+            values["right_edge"] = values["left_edge"] + (
+                values["nx"] * values["dx"]
+            )
+        if values.get("dx") is None:
+            values["dx"] = (values["right_edge"] - values["left_edge"]) / (
+                values["nx"]
+            )
+        if values.get("dy") is None:
+            values["dy"] = (values["top_edge"] - values["bottom_edge"]) / (
+                values["ny"]
+            )
+        # Be a bit more careful here that our grid generates a near integer
+        for key, spacing, edges in [
+            ("nx", "dx", ["left_edge", "right_edge"]),
+            ("ny", "dy", ["bottom_edge", "top_edge"]),
+        ]:
+            if values.get(key) is not None:
+                continue
+            nn = (values[edges[1]] - values[edges[0]]) / values[spacing]
+            if abs(nn - int(nn)) > 0.01:
+                msg = f"Computed {key} is not approximately an integer"
+                raise ValueError(msg)
+            values[key] = int(nn)
+
+        return values
+
+    def find_ij(
+        self, lon: float, lat: float
+    ) -> tuple[Optional[int], Optional[int]]:
+        """Find the grid cell that contains the given lon/lat (EPSG: 4326)."""
+        x, y = Proj(self.crs)(lon, lat)
+        if (
+            x < self.left_edge
+            or x >= self.right_edge
+            or y < self.bottom_edge
+            or y >= self.top_edge
+        ):
+            return None, None
+        i = int((x - self.left_edge) / self.dx)
+        j = int((y - self.bottom_edge) / self.dy)
+        return i, j

tests/grid/test_grid_nav.py

@@ -0,0 +1,13 @@
+"""Test pyiem.grid.nav"""
+
+from pyiem.grid import nav
+
+
+def test_api():
+    """Test basic things."""
+    assert nav.IEMRE.bottom == 23.0
+
+
+def test_prism_calc():
+    """Test that PRISM works out to what we expect."""
+    assert (nav.PRISM.right - -66.50) < 0.01

tests/models/test_models_gridnav.py

@@ -0,0 +1,38 @@
+"""Test the gridnav module."""
+
+import pytest
+
+from pyiem.models.gridnav import CartesianGridNavigation
+
+
+@pytest.fixture
+def cgn() -> CartesianGridNavigation:
+    """Return a basic CartesianGridNavigation."""
+    return CartesianGridNavigation(
+        left_edge=0,
+        bottom_edge=0,
+        dx=1,
+        dy=1,
+        nx=10,
+        ny=10,
+    )
+
+
+def test_api(cgn):
+    """Test basic things."""
+    assert cgn.bottom == 0.5
+    assert len(cgn.x_points) == 10
+    assert len(cgn.y_points) == 10
+    assert len(cgn.x_edges) == 11
+    assert len(cgn.y_edges) == 11
+    assert cgn.right == cgn.x_points[-1]
+
+
+def test_find_ij(cgn):
+    """See if we can get the right cell."""
+    i, j = cgn.find_ij(0.5, 0.5)
+    assert i == 0
+    assert j == 0
+    i, j = cgn.find_ij(0.5, 10.5)
+    assert i is None
+    assert j is None