diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a548c00..e025d97f 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/MANIFEST.in b/MANIFEST.in index 733eb766..3c53b4f4 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,3 +1,4 @@ recursive-include pyiem *.py +recursive-include pyiem *.pyi recursive-include pyiem *.txt recursive-include pyiem *.npy diff --git a/src/pyiem/grid/nav.py b/src/pyiem/grid/nav.py new file mode 100644 index 00000000..cc729416 --- /dev/null +++ b/src/pyiem/grid/nav.py @@ -0,0 +1,109 @@ +"""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 get_nav(name: str, dom: str) -> CartesianGridNavigation: + """Helper to remove some boilerplate for fetching gridnav.""" + extra = f"_{dom.upper()}" if dom != "" else "" + key = f"{name.upper()}{extra}" + return CartesianGridNavigation(**_GRID_CONFIGS[key]) + + +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}'") diff --git a/src/pyiem/grid/nav.pyi b/src/pyiem/grid/nav.pyi new file mode 100644 index 00000000..c34ea8a0 --- /dev/null +++ b/src/pyiem/grid/nav.pyi @@ -0,0 +1,16 @@ +"""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): ... + +def get_nav(name: str, region: str) -> CartesianGridNavigation: ... diff --git a/src/pyiem/models/gridnav.py b/src/pyiem/models/gridnav.py new file mode 100644 index 00000000..1e17de52 --- /dev/null +++ b/src/pyiem/models/gridnav.py @@ -0,0 +1,164 @@ +"""Grid Navigation Metadata.""" + +from typing import Optional, Union + +import numpy as np +from affine import Affine +from pydantic import BaseModel, ConfigDict, Field, model_validator +from pyproj import CRS, Proj + + +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 Affine(self.dx, 0, self.left_edge, 0, self.dy, self.bottom_edge) + + @property + def affine_image(self): + """Return the transformation associated with upper left origin.""" + return Affine( + self.dx, 0, self.left_edge, 0, 0 - self.dy, self.top_edge + ) + + @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) # skipcq + 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 diff --git a/tests/grid/test_grid_nav.py b/tests/grid/test_grid_nav.py new file mode 100644 index 00000000..bcdbee05 --- /dev/null +++ b/tests/grid/test_grid_nav.py @@ -0,0 +1,22 @@ +"""Test pyiem.grid.nav""" + +from pyiem.grid import nav + + +def test_api(): + """Test basic things.""" + assert nav.IEMRE.bottom == 23.0 + assert nav.IEMRE_CHINA.top == 54.875 + assert nav.IEMRE.affine + assert nav.IEMRE_EUROPE.affine_image + + +def test_prism_calc(): + """Test that PRISM works out to what we expect.""" + assert (nav.PRISM.right - -66.50) < 0.01 + + +def test_get_nav(): + """Test that we can get a helper.""" + assert nav.get_nav("iemRE", "") + assert nav.get_nav("era5LAND", "china") diff --git a/tests/models/test_models_gridnav.py b/tests/models/test_models_gridnav.py new file mode 100644 index 00000000..709ada20 --- /dev/null +++ b/tests/models/test_models_gridnav.py @@ -0,0 +1,79 @@ +"""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_non_even_spacing(): + """Test that this errors.""" + with pytest.raises(ValueError): + CartesianGridNavigation( + left_edge=0, + bottom_edge=0, + top_edge=10, + right_edge=10, + dx=2.2, + dy=0.7, + ) + + +def test_computing_dxdy(): + """Test that we can get a dx and dy.""" + _cgn = CartesianGridNavigation( + left_edge=0, + bottom_edge=0, + nx=10, + ny=10, + right_edge=10, + top_edge=10, + ) + assert _cgn.dx == 1 + assert _cgn.dy == 1 + + +def test_computing_nxny(): + """Test that ny and ny can be computed.""" + _cgn = CartesianGridNavigation( + left_edge=0, + bottom_edge=0, + dx=1, + dy=1, + right_edge=10, + top_edge=10, + ) + assert _cgn.nx == 10 + assert _cgn.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