Merge pull request #72 from C2SM-RCM/coli

Coli

docs/source/api/grids.rst

@@ -8,6 +8,8 @@ General Grids
 
 .. autoclass:: emiproc.grids.RegularGrid
 
+.. autoclass:: emiproc.grids.HexGrid
+
 .. autoclass:: emiproc.grids.GeoPandasGrid
 
 

docs/source/tutos/tutorials.rst

@@ -34,3 +34,5 @@ The first tutorial to understand how emiproc works is the edgar_processing tutor
     edgar_processing
     gfed
     icon_oem
+    grids
+    regridding

emiproc/grids.py

@@ -6,6 +6,7 @@
 
 from __future__ import annotations
 
+import logging
 import math
 from pathlib import Path
 import warnings
@@ -30,6 +31,8 @@
 R_EARTH = 6371000  # m
 
 
+logger = logging.getLogger(__name__)
+
 # Type alias
 # minx, miny, maxx, maxy
 BoundingBox = tuple[float, float, float, float]
@@ -366,6 +369,141 @@ def from_centers(
         return grid
 
 
+class HexGrid(Grid):
+    """A grid with hexagonal cells.
+
+    The grid is similar to a regular grid, but the cells are hexagons.
+    This implies that the lines of the grid are not all parallel.
+
+    For the arguments, look at :py:class:`RegularGrid`.
+    Note that `dx and dy` have been replaced by `spacing` parameter,
+
+    :arg spacing: The distance between the centers of two adjacent hexagons
+        on a row. This is also the diameter of the circle inscribed in the hexagon.
+    :arg oriented_north: If True, the hexagons are oriented with the top
+        and bottom sides parallel on the y-axis. If False, the hexagons
+        are oriented with the left and right sides parallel on the x-axis.
+
+    """
+
+    def __init__(
+        self,
+        xmin: float,
+        ymin: float,
+        xmax: float | None = None,
+        ymax: float | None = None,
+        nx: int | None = None,
+        ny: int | None = None,
+        spacing: float | None = None,
+        oriented_north: bool = True,
+        name: str | None = None,
+        crs: int | str = WGS84,
+    ):
+
+        # Correct the delta in case the orientation is on the other ax
+        correct = lambda x: x * np.sqrt(3) / 2
+
+        self.xmin, self.ymin = xmin, ymin
+
+        # Check if they did not specify all the optional parameters
+        if all((p is not None for p in [xmax, ymax, nx, ny, spacing])):
+            raise ValueError(
+                "Specified too many parameters. "
+                "Specify only 2 of the following: "
+                "(xmax, ymax), (nx, ny), (spacing)"
+            )
+
+        if spacing is None and xmax is None and ymax is None:
+            raise ValueError(
+                "Cannot create grid with only nx and ny. "
+                "Specify at least spacing or xmax, ymax."
+            )
+
+        if spacing is not None:
+            # Guess the nx and ny values, override the max
+            dx = spacing if oriented_north else correct(spacing)
+            dy = correct(spacing) if oriented_north else spacing
+
+        # Calclate the number of cells if not specified
+        if nx is None and ny is None:
+            if spacing is None:
+                raise ValueError(
+                    "Either nx and ny or dx and dy must be specified. "
+                    f"Received: {nx=}, {ny=}, {spacing=}"
+                )
+            if xmax is None or ymax is None:
+                raise ValueError(
+                    "When using only dx dy, xmax and ymax must be specified. "
+                    f"Received: {xmax=}, {ymax=}"
+                )
+
+            nx = math.ceil((xmax - xmin) / dx)
+            ny = math.ceil((ymax - ymin) / dy)
+
+        if xmax is None and ymax is None:
+            xmax = xmin + nx * dx
+            ymax = ymin + ny * dy
+        self.xmax, self.ymax = xmax, ymax
+
+        # Calculate all grid parameters
+        self.nx, self.ny = nx, ny
+        self.dx, self.dy = (xmax - xmin) / nx, (ymax - ymin) / ny
+
+        self.lon_range = np.arange(xmin, xmax, self.dx) + self.dx / 2
+        self.lat_range = np.arange(ymin, ymax, self.dy) + self.dy / 2
+
+        assert len(self.lon_range) == nx, f"{len(self.lon_range)=} != {nx=}"
+        assert len(self.lat_range) == ny, f"{len(self.lat_range)=} != {ny=}"
+
+        if name is None:
+            name = f"HexGrid x({xmin},{xmax})_y({ymin},{ymax})_nx({nx})_ny({ny})"
+
+        self.oriented_north = oriented_north
+
+        super().__init__(name, crs)
+
+    @cached_property
+    def cells_as_polylist(self) -> list[Polygon]:
+        """Return all the cells as a list of polygons."""
+
+        x_centers, y_centers = np.meshgrid(self.lon_range, self.lat_range)
+        # Shift the odd rows
+        if self.oriented_north:
+            x_centers[1::2] += self.dx / 2
+        else:
+            y_centers[:, 1::2] += self.dy / 2
+
+        # Calculat the corners of the hexagons
+        corners = []
+        flatten = lambda x: x.flatten(order="F")
+        x_centers = flatten(x_centers)
+        y_centers = flatten(y_centers)
+
+        # half_offset = np.sqrt(3) / 2
+        # half_offset = np.sqrt(3) / 8
+        half_offset = 1 / (np.sqrt(3))
+        offsets_x = [0, 1, 1, 0, -1, -1]
+        offsets_y = [
+            2 - half_offset,
+            half_offset,
+            -half_offset,
+            -2 + half_offset,
+            -half_offset,
+            half_offset,
+        ]
+        if not self.oriented_north:
+            offsets_x, offsets_y = offsets_y, offsets_x
+
+        for off_x, off_y in zip(offsets_x, offsets_y):
+            x = x_centers + self.dx / 2 * off_x
+            y = y_centers + self.dy / 2 * off_y
+            corners.append(np.stack([x, y], axis=-1))
+        # Reshape to 1D (order set for backward compatibility)
+
+        coords = np.stack(corners, axis=1)
+        return polygons(coords)
+
+
 class TNOGrid(RegularGrid):
     """Contains the grid from the TNO emission inventory
     This grid is defined as a standard lat/lon coordinate system.

emiproc/tests_utils/test_grids.py

@@ -2,7 +2,7 @@
 
 from shapely.geometry import Point, Polygon
 
-from emiproc.grids import RegularGrid, GeoPandasGrid
+from emiproc.grids import HexGrid, RegularGrid, GeoPandasGrid
 
 # Regular grid for testing
 # Note that this grid is big enough to include the inventories from the
@@ -11,6 +11,8 @@
     xmin=-1, xmax=5, ymin=-2, ymax=3, nx=10, ny=15, name="Test Regular Grid"
 )
 
+hex_grid = HexGrid(xmin=-1, xmax=5, ymin=-2, ymax=3, nx=10, ny=15, name="Test Hex Grid")
+
 basic_serie = gpd.GeoSeries(
     [
         Polygon(((0, 0), (0, 1), (1, 1), (1, 0))),

tests/test_hex_grid.py

@@ -0,0 +1,91 @@
+# %%
+import pytest
+import geopandas as gpd
+from emiproc.grids import HexGrid
+from emiproc.tests_utils.test_grids import hex_grid
+
+
+def test_creation():
+    HexGrid(xmin=-1, xmax=5, ymin=-2, ymax=3, nx=10, ny=15, name="Test Regular Grid")
+
+
+def test_creation_orientation():
+    HexGrid(
+        xmin=-1,
+        xmax=5,
+        ymin=-2,
+        ymax=3,
+        nx=10,
+        ny=15,
+        name="Test Regular Grid",
+        oriented_north=False,
+    )
+
+
+def test_creation_n_d():
+    HexGrid(xmin=-1, ymin=-2, spacing=0.1, nx=10, ny=15, name="Test Regular Grid")
+
+
+def test_creation_max_d():
+    HexGrid(xmin=-1, ymin=-2, xmax=5, ymax=3, spacing=0.5, name="Test Regular Grid")
+
+
+def test_creation_fails_max_d_n():
+    pytest.raises(
+        ValueError,
+        HexGrid,
+        xmin=-1,
+        ymin=-2,
+        xmax=5,
+        ymax=3,
+        spacing=0.2,
+        nx=10,
+        ny=15,
+        name="Test Regular Grid",
+    )
+
+
+def test_creation_fails_not_enough():
+    pytest.raises(
+        ValueError,
+        HexGrid,
+        xmin=-1,
+        ymin=-2,
+        xmax=5,
+        ymax=3,
+        name="Test Regular Grid",
+    )
+    pytest.raises(
+        ValueError,
+        HexGrid,
+        xmin=-1,
+        ymin=-2,
+        spacing=0.1,
+        name="Test Hex  Grid",
+    )
+    pytest.raises(
+        ValueError,
+        HexGrid,
+        xmin=-1,
+        ymin=-2,
+        nx=10,
+        ny=15,
+        name="Test Hex  Grid",
+    )
+
+
+def test_polylist():
+    hex_grid.cells_as_polylist
+
+    # Test how we iterate over the cells in the polygon list
+
+
+def test_area():
+    hex_grid.cell_areas
+
+
+def test_shape():
+    nx, ny = hex_grid.shape
+
+    assert nx == hex_grid.nx
+    assert ny == hex_grid.ny