Feat/186 raster georeference

CHANGELOG.md

@@ -6,6 +6,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ---
 
+## [Unreleased]
+
+### Changed
+- `load_geotiff` sets `RasterBackend.transform` and `RasterBackend.crs` from
+  the file (they were only in the returned `meta` dict), and `save_geotiff`
+  falls back to them: `save_geotiff(backend, uri)` writes a georeferenced
+  GeoTIFF without a `meta` dict. Explicit `crs`/`transform` arguments still
+  override `meta`, and `meta` overrides the backend.
+
+### Added
+- `RasterBackend.cell_area()`: the area of each cell from `transform` and
+  `crs` — square metres on the ellipsoid for a geographic CRS (a 1/12° cell
+  is ~86 km² at the equator and ~72 km² at 33° S), the pixel area for a
+  projected one.
+
 ## [0.6.5] — 2026-09-22
 
 ### Fixed

dissmodel/geo/raster/backend.py

@@ -498,6 +498,40 @@ def from_xarray(cls, ds, nodata_value: float | None = None) -> RasterBackend:
 
     # ── spatial operations ────────────────────────────────────────────────────
 
+    def cell_area(self) -> np.ndarray:
+        """
+        Area of each cell, shape ``(rows, cols)``, from ``transform`` and ``crs``.
+
+        On a geographic CRS the result is in **square metres** on the CRS's
+        ellipsoid (``pyproj.Geod``): a 1/12° cell is ~86 km² at the equator and
+        ~72 km² at 34° S, so a fraction map's areas are fraction × cell area,
+        not fraction × a nominal size. On a projected CRS it is the planar
+        pixel area, in the CRS's units squared.
+
+        Requires a north-up ``transform`` (no rotation) and, for a geographic
+        CRS, ``pyproj``.
+        """
+        if self.transform is None or self.crs is None:
+            raise ValueError("cell_area needs the backend's transform and crs")
+        t = self.transform
+        if t.b != 0 or t.d != 0:
+            raise ValueError("cell_area supports north-up transforms only (no rotation)")
+        rows, cols = self.shape
+        from pyproj import CRS, Geod
+
+        crs = CRS.from_user_input(self.crs)
+        if not crs.is_geographic:
+            return np.full((rows, cols), abs(t.a * t.e))
+        geod = crs.get_geod() or Geod(ellps="WGS84")
+        west, east = t.c, t.c + t.a
+        area = np.empty(rows)
+        for r in range(rows):
+            north, south = t.f + r * t.e, t.f + (r + 1) * t.e
+            a, _ = geod.polygon_area_perimeter([west, east, east, west], [north, north, south, south])
+            area[r] = abs(a)
+        return np.repeat(area[:, None], cols, axis=1)
+
+
     @staticmethod
     def shift2d(arr: np.ndarray, dr: int, dc: int) -> np.ndarray:
         """

dissmodel/io/raster.py

@@ -80,7 +80,7 @@ def _read_geotiff(
 
     with rasterio.open(path_str) as ds:
         rows, cols = ds.height, ds.width
-        backend    = RasterBackend(shape=(rows, cols))
+        backend    = RasterBackend(shape=(rows, cols), transform=ds.transform, crs=ds.crs)
 
         if band_spec:
             for i, (name, dtype, nodata) in enumerate(band_spec, start=1):
@@ -126,17 +126,20 @@ def save_geotiff(
 
     Parameters
     ----------
-    data : (RasterBackend, dict)
-        Backend and metadata dict (as returned by load_geotiff).
+    data : (RasterBackend, dict) or RasterBackend
+        Backend and metadata dict (as returned by load_geotiff), or the
+        backend alone.
     uri : str
         Destination URI. Local path or s3://bucket/key.
     band_spec : list of (name, dtype, nodata) or None
         Bands to write in order. Missing bands are filled with nodata.
         If None, all arrays in the backend are written.
     crs : str or None
-        CRS string (e.g. "EPSG:31984"). Overrides meta["crs"].
+        CRS string (e.g. "EPSG:31984"). Overrides meta["crs"], which
+        overrides ``backend.crs``.
     transform : Affine or None
-        Affine geotransform. Overrides meta["transform"].
+        Affine geotransform. Overrides meta["transform"], which overrides
+        ``backend.transform``.
     compress : str
         Compression algorithm. Default: "deflate".
 
@@ -150,9 +153,9 @@ def save_geotiff(
     if not HAS_RASTERIO:
         raise ImportError("rasterio is required — pip install rasterio")
 
-    backend, meta = data
-    resolved_crs       = crs       or (meta.get("crs") if meta else None)
-    resolved_transform = transform or (meta.get("transform") if meta else None)
+    backend, meta = data if isinstance(data, tuple) else (data, None)
+    resolved_crs       = crs       or (meta.get("crs") if meta else None) or backend.crs
+    resolved_transform = transform or (meta.get("transform") if meta else None) or backend.transform
 
     with tempfile.NamedTemporaryFile(suffix=".tif", delete=False) as f:
         tmp = f.name

paper.md

@@ -179,7 +179,7 @@ year-by-year reference outputs generated in a containerised TerraME
 that perturbing the regression coefficients breaks the tolerance criterion.
 
 The discrete CLUE-S-like allocation in `disslucc`, with a logistic-regression
-potential, reproduces the Lab15 case study (Moju municipality, 5,914 cells, 6 steps)
+potential, reproduces the Lab15 case study (Mojui, Pará 5,914 cells, 6 steps)
 from the reference LuccME implementation [@LuccME] cell for cell — zero quantity and
 zero allocation disagreement [@PontiusMillones2011] — at 10.3 ms/step. A shipped
 discriminance test shows the final map is also reproduced by a trivial static ranking,

tests/io/test_raster_io.py

@@ -120,6 +120,22 @@ def test_s3_upload(self, backend):
         assert len(payload) > 0
         assert isinstance(checksum, str) and len(checksum) == 64
 
+    def test_backend_georeference_is_the_fallback(self, backend, tmp_path):
+        path = tmp_path / "own.tif"
+        backend.transform = rasterio.transform.from_bounds(10, 20, 510, 420, 5, 4)
+        backend.crs = "EPSG:31984"
+        save_geotiff(backend, str(path))  # the backend alone, no meta
+        with rasterio.open(str(path)) as ds:
+            assert "31984" in ds.crs.to_string()
+            assert ds.transform == backend.transform
+
+    def test_meta_overrides_the_backend_georeference(self, backend, tmp_path):
+        path = tmp_path / "meta_wins.tif"
+        backend.crs = "EPSG:31984"
+        save_geotiff((backend, {"crs": "EPSG:4326"}), str(path))
+        with rasterio.open(str(path)) as ds:
+            assert ds.crs.to_string() == "EPSG:4326"
+
     def test_mixed_dtype_bands_are_not_truncated(self, backend, tmp_path):
         """Regression: int32 + float32 bands must promote to a common
         dtype instead of truncating floats to the first band's dtype."""
@@ -153,6 +169,20 @@ def test_roundtrip_recovers_crs_and_transform(self, saved_tif):
         assert "31984" in meta["crs"].to_string()
         assert meta["transform"] is not None
 
+    def test_backend_carries_crs_and_transform(self, saved_tif):
+        path, band_spec, _ = saved_tif
+        (loaded, meta), _ = load_geotiff(str(path), band_spec=band_spec)
+        assert loaded.crs == meta["crs"]
+        assert loaded.transform == meta["transform"]
+
+    def test_load_then_save_the_backend_keeps_the_georeference(self, saved_tif, tmp_path):
+        path, _, _ = saved_tif
+        (loaded, _), _ = load_geotiff(str(path))
+        copy = tmp_path / "copy.tif"
+        save_geotiff(loaded, str(copy))
+        with rasterio.open(str(path)) as a, rasterio.open(str(copy)) as b:
+            assert a.crs == b.crs and a.transform == b.transform
+
     def test_without_band_spec_recovers_tag_names(self, saved_tif):
         path, _, _ = saved_tif
         (loaded, _), _ = load_geotiff(str(path))

tests/raster/test_backend.py

@@ -230,3 +230,38 @@ def test_moore_includes_diagonals(self):
     def test_von_neumann_no_diagonals(self):
         diagonals = [(dr, dc) for dr, dc in DIRS_VON_NEUMANN if dr != 0 and dc != 0]
         assert diagonals == []
+
+
+class TestCellArea:
+    """cell_area(): true areas on a geographic grid, pixel area on a projected one."""
+
+    def _geographic(self, res: float, top: float, rows: int = 3, cols: int = 2) -> RasterBackend:
+        pytest.importorskip("pyproj")
+        from affine import Affine
+
+        return RasterBackend(shape=(rows, cols), transform=Affine(res, 0, -50.0, 0, -res, top), crs="EPSG:4326")
+
+    def test_one_degree_at_the_equator(self):
+        b = self._geographic(1.0, 0.5, rows=1, cols=1)
+        # 1° × 1° centred on the equator, WGS 84: 110.57 km × 111.32 km
+        assert b.cell_area()[0, 0] == pytest.approx(12_309e6, rel=1e-3)
+
+    def test_cells_shrink_away_from_the_equator_and_rows_repeat(self):
+        b = self._geographic(1 / 12, 0.0, rows=400)
+        area = b.cell_area()
+        assert area.shape == (400, 2)
+        assert (np.diff(area[:, 0]) < 0).all()
+        assert (area[:, 0] == area[:, 1]).all()
+        assert area[0, 0] / 1e6 == pytest.approx(85.5, abs=0.1)     # 1/12° at the equator
+        assert area[-1, 0] / 1e6 == pytest.approx(71.7, abs=0.1)    # 1/12° at 33.3° S
+
+    def test_projected_grid_is_the_pixel_area(self):
+        pytest.importorskip("pyproj")
+        from affine import Affine
+
+        b = RasterBackend(shape=(2, 3), transform=Affine(30.0, 0, 0, 0, -30.0, 0), crs="EPSG:31983")
+        np.testing.assert_array_equal(b.cell_area(), np.full((2, 3), 900.0))
+
+    def test_needs_the_georeference(self):
+        with pytest.raises(ValueError, match="transform and crs"):
+            RasterBackend(shape=(2, 2)).cell_area()