fix: transform points to the DEM CRS before sampling the DEM

The DEM was sampled with model CRS coordinates. A DEM in another CRS
returned no data for every point, so all points got an elevation of 0.

Author: lachlangrose

loopstructural/main/data_manager.py

@@ -355,6 +355,26 @@ def _refresh_dem_function(self):
         `dem_layer`, even when `use_dem` is True).
         """
         if self.use_dem and self.dem_layer is not None:
+            # Callers pass (x, y) in the model CRS, but `sample()` expects the
+            # DEM's own CRS. Cache the transform and rebuild it only when the
+            # model or DEM CRS changes, because this is called once per point.
+            cache = {'model_crs': None, 'dem_crs': None, 'transform': None}
+
+            def to_dem_crs(x, y):
+                model_crs = self.get_model_crs()
+                dem_crs = self.dem_layer.crs()
+                if (
+                    model_crs is None
+                    or not model_crs.isValid()
+                    or not dem_crs.isValid()
+                    or model_crs == dem_crs
+                ):
+                    return QgsPointXY(x, y)
+                if cache['model_crs'] != model_crs or cache['dem_crs'] != dem_crs:
+                    cache['model_crs'] = model_crs
+                    cache['dem_crs'] = dem_crs
+                    cache['transform'] = QgsCoordinateTransform(model_crs, dem_crs, self.project)
+                return cache['transform'].transform(QgsPointXY(x, y))
 
             def dem_function(x, y):
                 if not self.dem_layer.isValid():
@@ -363,7 +383,7 @@ def dem_function(x, y):
                         log_level=2,
                     )
                     return 0.0
-                value = self.dem_layer.dataProvider().sample(QgsPointXY(x, y), 1)[0]
+                value = self.dem_layer.dataProvider().sample(to_dem_crs(x, y), 1)[0]
                 # `sample()` returns NaN (not None) for points outside the
                 # raster's extent or on nodata cells, depending on provider.
                 if value is None or np.isnan(value):

tests/qgis/test_dem_function.py

@@ -0,0 +1,66 @@
+from unittest.mock import Mock
+
+import numpy as np
+import pytest
+from osgeo import gdal, osr
+from qgis.core import (
+    QgsCoordinateReferenceSystem,
+    QgsCoordinateTransform,
+    QgsPointXY,
+    QgsProject,
+    QgsRasterLayer,
+)
+
+from loopstructural.main.data_manager import ModellingDataManager
+
+DEM_VALUE = 460.0
+MODEL_CRS = "EPSG:32755"
+
+
+def _write_geographic_dem(path):
+    """Write a small constant-value GeoTIFF in EPSG:4326 over (147E, 42S)."""
+    driver = gdal.GetDriverByName("GTiff")
+    dataset = driver.Create(str(path), 10, 10, 1, gdal.GDT_Float32)
+    dataset.SetGeoTransform((146.9, 0.02, 0.0, -41.9, 0.0, -0.02))
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(4326)
+    dataset.SetProjection(srs.ExportToWkt())
+    dataset.GetRasterBand(1).WriteArray(np.full((10, 10), DEM_VALUE, dtype=np.float32))
+    dataset.FlushCache()
+    dataset = None
+
+
+@pytest.fixture
+def data_manager(tmp_path):
+    path = tmp_path / "dem.tif"
+    _write_geographic_dem(path)
+    dem_layer = QgsRasterLayer(str(path), "dem")
+    assert dem_layer.isValid()
+    manager = ModellingDataManager(
+        project=QgsProject.instance(), mapCanvas=Mock(), logger=Mock()
+    )
+    manager.set_model_manager(Mock())
+    manager.set_model_crs(QgsCoordinateReferenceSystem(MODEL_CRS), use_project_crs=False)
+    manager.set_dem_layer(dem_layer)
+    manager.set_use_dem(True)
+    return manager
+
+
+def _dem_centre_in_model_crs():
+    transform = QgsCoordinateTransform(
+        QgsCoordinateReferenceSystem("EPSG:4326"),
+        QgsCoordinateReferenceSystem(MODEL_CRS),
+        QgsProject.instance(),
+    )
+    return transform.transform(QgsPointXY(147.0, -42.0))
+
+
+def test_dem_in_another_crs_is_sampled_at_model_coordinates(data_manager):
+    """The DEM is sampled after the model (x, y) is transformed to the DEM CRS."""
+    centre = _dem_centre_in_model_crs()
+    assert data_manager.dem_function(centre.x(), centre.y()) == pytest.approx(DEM_VALUE)
+
+
+def test_point_outside_the_dem_gives_zero(data_manager):
+    centre = _dem_centre_in_model_crs()
+    assert data_manager.dem_function(centre.x() + 1e6, centre.y()) == 0.0