Geo API

Geo Package

Generic geospatial logic and data for the NGED substation forecast project: H3 spatial indexing and the Great Britain boundary the numerical weather prediction (NWP) grid is clipped to. H3 is a grid system that tiles the globe in hexagons at nested resolutions, with 12 pentagons where hexagons alone cannot close the sphere, and each cell is identified by an index.

Map of Great Britain using H3 resolution 5 hexagons

Purpose

The geo package decouples generic geospatial operations from dataset-specific ingestion logic, such as the processing of European Centre for Medium-Range Weather Forecasts (ECMWF) data in dynamical_data. Any package in the workspace can therefore perform a spatial transformation — mapping a latitude/longitude grid to H3 hexagons, for example — without depending on heavy or unrelated packages.

compute_h3_grid_weights_for_boundary accepts any boundary polygon, not only the Great Britain shape this package ships. Accepting any boundary polygon adds no complexity here and means a new region plugs into the same H3 gridding rather than forking the gridding code. Keeping the function boundary-agnostic is design principle 5 applied to this package.

Two neighbouring jobs are deliberately not here. The per-substation H3 index (h3_res_5 on TimeSeriesMetadata) is computed by nged_data straight from each substation's coordinates. The spatial aggregation that consumes the grid weights computed here happens in dynamical_data at ECMWF ingest. The weather grid is square and the H3 grid is hexagonal, so a grid weight is the fraction of one hexagon that one square grid point covers, and each hexagon takes a weighted share of the grid points overlapping it.

Contents

• h3 — compute_h3_grid_weights_for_boundary() and compute_h3_grid_weights(), which map the H3 grid onto the regular lat/lon NWP grid. The sampling and snapping method is documented on the functions themselves.
• great_britain.load — load_gb_boundary(), which loads the Great Britain boundary polygon that the NWP grid is clipped to, from the packaged GeoJSON file.

geo.h3

H3-related utilities for geospatial operations.

Classes

Functions:

compute_h3_grid_weights_for_boundary(boundary, nwp_grid_size_degrees, h3_res, child_h3_res=None)

Computes the H3 grid weights for a geospatial boundary.

Generates the spatial mapping between the hexagonal H3 grid and the regular lat/lon NWP grid.

The mapping is calculated by sampling each H3 cell with finer-resolution child cells and determining which regular grid cell each child falls into. The nwp_grid_size_degrees parameter is used to snap high-resolution H3 cells to the nearest regular NWP grid points.

Parameters:

Name	Type	Description	Default
boundary	BaseGeometry	The geospatial boundary to find H3 cells for.	required
nwp_grid_size_degrees	float	The size of the regular lat/lon grid in degrees.	required
h3_res	int	The H3 resolution to use for the grid.	required
child_h3_res	int | None	The H3 resolution to use for the underlying points. If None, it defaults to h3_res + 2.	None

Returns:

Type	Description
DataFrame[H3GridWeights]	One row per (H3 cell, NWP grid point) pair that overlap within boundary. proportion
DataFrame[H3GridWeights]	holds the fraction of that H3 cell's child cells falling inside the grid point's box. This
DataFrame[H3GridWeights]	function resolves boundary to its covering h3_index list and then delegates to
DataFrame[H3GridWeights]	compute_h3_grid_weights — see that function for the rest.

Source code in packages/geo/src/geo/h3.py

def compute_h3_grid_weights_for_boundary(
    boundary: BaseGeometry,
    nwp_grid_size_degrees: float,
    h3_res: int,
    child_h3_res: int | None = None,
) -> pt.DataFrame[H3GridWeights]:
    """Computes the H3 grid weights for a geospatial boundary.

    Generates the spatial mapping between the hexagonal H3 grid and the regular lat/lon NWP grid.

    The mapping is calculated by sampling each H3 cell with finer-resolution child cells and
    determining which regular grid cell each child falls into. The `nwp_grid_size_degrees`
    parameter is used to snap high-resolution H3 cells to the nearest regular NWP grid points.

    Args:
        boundary: The geospatial boundary to find H3 cells for.
        nwp_grid_size_degrees: The size of the regular lat/lon grid in degrees.
        h3_res: The H3 resolution to use for the grid.
        child_h3_res: The H3 resolution to use for the underlying points. If None,
            it defaults to h3_res + 2.

    Returns:
        One row per (H3 cell, NWP grid point) pair that overlap within `boundary`. `proportion`
        holds the fraction of that H3 cell's child cells falling inside the grid point's box. This
        function resolves `boundary` to its covering `h3_index` list and then delegates to
        `compute_h3_grid_weights` — see that function for the rest.
    """
    _LOG.info(f"Generating H3 cells at resolution {h3_res}...")

    h3_index = h3.geo_to_cells(geo=boundary, res=h3_res)
    if not h3_index:
        raise ValueError(
            f"No H3 cells found for the given boundary at resolution {h3_res}."
            " Check if the boundary geometry is valid and covers the expected area."
        )

    return compute_h3_grid_weights(
        nwp_grid_size_degrees=nwp_grid_size_degrees, h3_index=h3_index, child_h3_res=child_h3_res
    )

compute_h3_grid_weights(nwp_grid_size_degrees, h3_index, child_h3_res=None)

Computes the proportion mapping for H3 grid cells to a regular lat/lng grid.

This function takes a list of H3 indices. Each H3 resolution subdivides the one above it, so a cell at resolution n contains a set of smaller child cells at resolution n+1. For each index, this function counts how many of its child H3 cells at the finer resolution child_h3_res fall into each cell of a regular lat/lng grid of size nwp_grid_size_degrees.

The regular grid is assumed to be perfectly aligned to 0.0 (e.g., 0.0, nwp_grid_size_degrees, nwp_grid_size_degrees * 2).

Parameters:

Name	Type	Description	Default
nwp_grid_size_degrees	float	The size of the regular NWP lat/lng grid in degrees (e.g., 0.25).	required
h3_index	list[int]	List of 64-bit H3 discrete spatial indices.	required
child_h3_res	int | None	The H3 resolution to use for the underlying points. Must be strictly greater than the resolution of the input h3_index list. If None, it defaults to that resolution + 2.	None

Returns:

Type	Description
DataFrame[H3GridWeights]	One row per (H3 cell, NWP grid point) pair that overlap, with proportion holding the
DataFrame[H3GridWeights]	fraction of that H3 cell's child cells falling inside the grid point's box. The
DataFrame[H3GridWeights]	proportion values for any one h3_index sum to 1.

Source code in packages/geo/src/geo/h3.py

def compute_h3_grid_weights(
    nwp_grid_size_degrees: float, h3_index: list[int], child_h3_res: int | None = None
) -> pt.DataFrame[H3GridWeights]:
    """Computes the proportion mapping for H3 grid cells to a regular lat/lng grid.

    This function takes a list of H3 indices. Each H3 resolution subdivides the one above it, so a
    cell at resolution *n* contains a set of smaller child cells at resolution *n+1*. For each
    index, this function counts how many of its child H3 cells at the finer resolution
    `child_h3_res` fall into each cell of a regular lat/lng grid of size `nwp_grid_size_degrees`.

    The regular grid is assumed to be perfectly aligned to 0.0 (e.g., 0.0,
    `nwp_grid_size_degrees`, `nwp_grid_size_degrees * 2`).

    Args:
        nwp_grid_size_degrees: The size of the regular NWP lat/lng grid in degrees (e.g., 0.25).
        h3_index: List of 64-bit H3 discrete spatial indices.
        child_h3_res: The H3 resolution to use for the underlying points. Must be
            strictly greater than the resolution of the input `h3_index` list. If None,
            it defaults to that resolution + 2.

    Returns:
        One row per (H3 cell, NWP grid point) pair that overlap, with `proportion` holding the
        fraction of that H3 cell's child cells falling inside the grid point's box. The
        `proportion` values for any one `h3_index` sum to 1.
    """
    # Reusable-package input validation, not a reachable production state: the only non-test caller
    # is `compute_h3_grid_weights_for_boundary` above, which has already raised on an empty cell
    # list. This guards direct callers of the public function; `test_h3.py` exercises it.
    if len(h3_index) == 0:
        raise ValueError("h3_index is empty.")

    # Ensure nwp_grid_size_degrees is strictly positive to avoid division by zero or
    # nonsensical snapping.
    if nwp_grid_size_degrees <= 0:
        raise ValueError(
            f"nwp_grid_size_degrees must be strictly positive, not {nwp_grid_size_degrees}."
        )

    df = pl.DataFrame({"h3_index": h3_index}, schema={"h3_index": pl.UInt64}).sort("h3_index")

    # Check resolution of all H3 indices to ensure consistency.
    h3_res_unique = df.select(plh3.get_resolution("h3_index")).unique()
    if h3_res_unique.height > 1:
        raise ValueError(f"All H3 indices must have the same resolution. {h3_res_unique=}")
    h3_res = h3_res_unique.item()

    # One H3 cell has about seven children at the next finer resolution, so stepping the resolution
    # up by two gives roughly 7^2 children. The `+2` heuristic therefore provides ~49 sample points
    # per H3 cell. That sample count balances spatial precision, for area-weighting against a
    # 0.25-degree grid, against computation time and memory overhead. Increasing the resolution gap
    # further could cause an exponential explosion in the number of child cells and potentially
    # trigger OOM errors.
    child_h3_res = h3_res + 2 if child_h3_res is None else child_h3_res

    if child_h3_res <= h3_res:
        raise ValueError(f"{child_h3_res=} must be strictly greater than {h3_res=}.")

    # Logged here rather than on entry so that child_h3_res is the resolution actually used: on the
    # default path the argument is None until the line above resolves it.
    _LOG.info(
        f"Computing H3 grid weights for grid size {nwp_grid_size_degrees} from {len(h3_index)}"
        f" H3 cells at resolution {h3_res}, sampled at child resolution {child_h3_res}..."
    )

    half_grid_size = nwp_grid_size_degrees / 2

    def _snap_to_grid(lat_or_lon: pl.Expr) -> pl.Expr:
        """Snap a latitude or longitude to the closest NWP grid line.

        The half-grid offset binning below ensures that points are snapped to the *closest* grid
        center rather than the bottom-left corner of the grid cell. Adding `half_grid_size`
        before flooring shifts the bin boundaries so that the grid points (0, 0.25, 0.5, etc.)
        are at the center of each bin.
        """
        return (
            (lat_or_lon + half_grid_size) / nwp_grid_size_degrees
        ).floor() * nwp_grid_size_degrees

    weights_df = (
        df.with_columns(child_h3=plh3.cell_to_children("h3_index", child_h3_res))
        # empty_as_null=False on the .explode("child_h3") below matches the Polars 2.0 default and
        # silences the deprecation warning. The setting has no effect on output today.
        # cell_to_children always returns a non-empty list here: child_h3_res > h3_res is enforced
        # above, and every H3 cell has children at any finer resolution. That holds for a pentagon
        # as well as a hexagon, and an H3 grid holds 12 pentagons because a sphere cannot be tiled
        # by hexagons alone. The empty-list branch the two settings disagree on is therefore
        # unreachable. An *invalid* index yields a null list, not an empty list. A null list
        # explodes to a single null row under both settings, and H3GridWeights.validate() below
        # rejects that row regardless. Neither setting therefore risks silent data loss.
        .explode("child_h3", empty_as_null=False)
        # nwp_lat takes cell_to_lat and nwp_lon takes cell_to_lng: swapping the two column
        # assignments below is the one orientation bug this expression can introduce. The swap
        # survives the row-count and weight-sum assertions in `test_grid_weights_invariant`, because
        # the output stays well-formed and only its geography is wrong. Two tests in
        # `geo/tests/test_h3.py` do catch the swap, and mutation testing verified both: the swap was
        # written into the code deliberately, and each test was watched failing.
        # `test_grid_weights_snap_to_nearest_grid_centre` builds its expected points geographically,
        # in its oracle. `test_grid_weights_preserve_geographic_orientation` pins two well-separated
        # landmarks in Great Britain. The orientation-coverage table covers the wider
        # set of orientation bugs this mapping can carry, and which test catches each:
        # <https://openclimatefix.github.io/nged-substation-forecast/architecture/testing/#nwp-grid-h3-orientation-coverage>
        .with_columns(
            nwp_lat=_snap_to_grid(plh3.cell_to_lat("child_h3")),
            nwp_lon=_snap_to_grid(plh3.cell_to_lng("child_h3")),
        )
        .group_by(["h3_index", "nwp_lat", "nwp_lon"])
        .len(name="h3_children_in_nwp_grid_box")
        .with_columns(
            h3_children_per_h3_parent=pl.col("h3_children_in_nwp_grid_box").sum().over("h3_index"),
        )
        .with_columns(
            proportion=pl.col("h3_children_in_nwp_grid_box") / pl.col("h3_children_per_h3_parent")
        )
        .sort(["h3_index", "nwp_lat", "nwp_lon"])
    )

    return pt.DataFrame(weights_df).set_model(H3GridWeights).drop().cast().validate()

geo.great_britain.load

Loading the Great Britain boundary polygon that the NWP H3 grid is clipped to.

Functions:

load_gb_boundary()

Loads the boundary geometry for Great Britain from a local GeoJSON file.

The boundary is buffered to ensure that coastal substations and nearby islands are included in the resulting H3 grid without spatial distortion.

Source code in packages/geo/src/geo/great_britain/load.py

def load_gb_boundary() -> BaseGeometry:
    """Loads the boundary geometry for Great Britain from a local GeoJSON file.

    The boundary is buffered to ensure that coastal substations and nearby islands are included
    in the resulting H3 grid without spatial distortion.
    """
    geojson_path = Path(__file__).parent / "england_scotland_wales.geojson"
    _LOG.info(f"Loading GB boundary from {geojson_path}")
    file_contents = geojson_path.read_text()
    shape: BaseGeometry = shapely.from_geojson(file_contents)
    return shape.buffer(distance=0.25)  # This takes about 30 seconds.