Merge pull request #51 from michaelohanrahan/fix/gridarea-projected-crs

fix: gridarea handles projected-CRS inputs

Author: sibrenloos

docs/changelog.rst

@@ -17,6 +17,13 @@ Changed
 
 Fixed
 -----
+- ``workflows.emissions.gridarea`` now branches on ``crs.is_projected`` and
+  delegates to ``ds.raster.area_grid()`` for projected-CRS inputs. The previous
+  implementation called ``raster_utils._reggrid_area`` unconditionally, which
+  assumes lat/lon in degrees and returns partly-negative cell areas when given
+  projected coordinates in metres — silently corrupting downstream mm->m3/s
+  hydrology forcing conversion and ``setup_emission_raster`` with
+  ``area_division=True``.
 
 v0.4.0 (15 April 2026)
 ========================

hydromt_delwaq/workflows/emissions.py

@@ -20,28 +20,37 @@
 def gridarea(ds):
     """Return a DataArray containing the area in m2 of the reference grid ds.
 
+    Works for both geographic (lat/lon degrees) and projected (metres) CRS.
+    The previous implementation called ``_reggrid_area`` unconditionally,
+    which treats x/y as lat/lon in degrees — fine for global hydromt setups,
+    but for projected CRS (e.g. SVY21 EPSG:3414, Easting/Northing in metres)
+    it computes ``sin(radians(Easting))`` and returns periodic, partly-negative
+    garbage. Delegating to ``ds.raster.area_grid()`` switches on
+    ``crs.is_projected`` and uses ``res_x * res_y * unit_factor**2`` instead.
+
     Parameters
     ----------
-    da : xarray.DataArray or xarray.DataSet
-        DataArray containing reference grid.
+    ds : xarray.DataArray or xarray.DataSet
+        Reference grid; must have a CRS set on the raster accessor.
 
     Returns
     -------
     da_out : xarray.DataArray
-        DataArray containing area in m2 of the reference grid.
-
+        Cell area in m2.
     """
+    crs = ds.raster.crs
+    if crs is not None and crs.is_projected:
+        return ds.raster.area_grid().astype("float32")
+
     realarea = raster_utils._reggrid_area(
         ds.raster.ycoords.values, ds.raster.xcoords.values
     )
-    da_out = xr.DataArray(
+    return xr.DataArray(
         data=realarea.astype("float32"),
         coords=ds.raster.coords,
         dims=ds.raster.dims,
     )
 
-    return da_out
-
 
 def gridlength_gridwidth(ds):
     """Return two DataArray for length and width of the grid."""

tests/test_model_methods.py

@@ -11,6 +11,52 @@
 TESTDATADIR = join(dirname(abspath(__file__)), "data")
 
 
+def _make_grid(crs, x0, y0, res, nx=4, ny=3):
+    x = x0 + (np.arange(nx) + 0.5) * res
+    y = y0 - (np.arange(ny) + 0.5) * res  # N -> S
+    da = xr.DataArray(
+        np.ones((ny, nx), dtype="float32"),
+        coords={"y": y, "x": x},
+        dims=("y", "x"),
+        name="dummy",
+    )
+    da.raster.set_crs(crs)
+    return da
+
+
+def test_gridarea_projected_returns_uniform_positive_area():
+    """Projected CRS (e.g. SVY21 / EPSG:3414) must produce uniform positive areas.
+
+    Regression test for the gridarea bug that fed SVY21 Easting metres into
+    ``_reggrid_area`` (which interprets them as longitude degrees). At
+    Easting=10500-18000 m, ``sin(radians(Easting))`` is a periodic function
+    of large angles and produces partly-negative garbage, which propagated
+    through the mm -> m3/s unit conversion and yielded negative precip values
+    in dynamicdata.nc. The patched gridarea delegates to
+    ``ds.raster.area_grid()`` for projected CRS, which uses res*res*ucf**2.
+    """
+    res = 30.0
+    da = _make_grid("EPSG:3414", x0=10490.0, y0=39280.0, res=res)
+
+    area = gridarea(da)
+
+    assert area.shape == da.shape
+    assert float(area.min()) > 0, "projected cell area must be strictly positive"
+    expected = res * res  # SVY21 linear unit factor is 1.0 (metres)
+    np.testing.assert_allclose(area.values, expected, rtol=1e-6)
+
+
+def test_gridarea_geographic_path_unchanged():
+    """Geographic CRS path must still go through _reggrid_area (spherical cap)."""
+    da = _make_grid("EPSG:4326", x0=4.0, y0=52.0, res=0.01)
+
+    area = gridarea(da)
+
+    assert float(area.min()) > 0
+    # At ~52N, 0.01deg cells are ~700-800 m on a side -> ~5-7e5 m2.
+    assert 4e5 < float(area.mean()) < 9e5
+
+
 def test_setup_grid(example_demission_model):
     # Initialize model and read results
     mod = example_demission_model