Remove RichDEM from source, moved utility functions to conftest.py (#612

)

CONTRIBUTING.md

@@ -50,6 +50,8 @@ pytest
 
 Running `pytest` will trigger a script that automatically downloads test data from [https://github.com/GlacioHack/xdem-data](https://github.com/GlacioHack/xdem-data) used to run all tests.
 
+RichDEM should only be used for testing purposes within the xDEM project. The functionality of xDEM must not depend on RichDEM.
+
 ### Documentation
 
 If your changes need to be reflected in the documentation, update the related pages located in `doc/source/`. The documentation is written in MyST markdown syntax, similar to GitHub's default Markdown (see [MyST-NB](https://myst-nb.readthedocs.io/en/latest/authoring/text-notebooks.html) for details).

dev-environment.yml

@@ -20,7 +20,6 @@ dependencies:
 
   # Optional dependencies
   - pytransform3d
-  - richdem
 
   # Test dependencies
   - gdal  # To test against GDAL
@@ -29,6 +28,7 @@ dependencies:
   - pyyaml
   - flake8
   - pylint
+  - richdem
 
   # Doc dependencies
   - sphinx

doc/source/terrain.md

@@ -10,9 +10,7 @@ and tested for consistency against [gdaldem](https://gdal.org/programs/gdaldem.h
 ## Quick use
 
 Terrain attribute methods can either be called directly from a {class}`~xdem.DEM` (e.g., {func}`xdem.DEM.slope`) or
-through the {class}`~xdem.terrain` module which allows array input. If computational performance
-is key, xDEM can rely on [RichDEM](https://richdem.readthedocs.io/) by specifying `use_richdem=True` for speed-up
-of its supported attributes (slope, aspect, curvature).
+through the {class}`~xdem.terrain` module which allows array input.
 
 ## Slope
 

setup.cfg

@@ -52,14 +52,14 @@ opt =
     opencv
     openh264
     pytransform3d
-    richdem
     noisyopt
 test =
     pytest
     pytest-xdist
     pyyaml
     flake8
     pylint
+    richdem
 doc =
     sphinx
     sphinx-book-theme

tests/conftest.py

@@ -0,0 +1,142 @@
+from typing import Callable, List, Union
+
+import geoutils as gu
+import numpy as np
+import pytest
+import richdem as rd
+from geoutils.raster import RasterType
+
+from xdem._typing import NDArrayf
+
+
+@pytest.fixture(scope="session")  # type: ignore
+def raster_to_rda() -> Callable[[RasterType], rd.rdarray]:
+    def _raster_to_rda(rst: RasterType) -> rd.rdarray:
+        """
+        Convert geoutils.Raster to richDEM rdarray.
+        """
+        arr = rst.data.filled(rst.nodata).squeeze()
+        rda = rd.rdarray(arr, no_data=rst.nodata)
+        rda.geotransform = rst.transform.to_gdal()
+        return rda
+
+    return _raster_to_rda
+
+
+@pytest.fixture(scope="session")  # type: ignore
+def get_terrainattr_richdem(raster_to_rda: Callable[[RasterType], rd.rdarray]) -> Callable[[RasterType, str], NDArrayf]:
+    def _get_terrainattr_richdem(rst: RasterType, attribute: str = "slope_radians") -> NDArrayf:
+        """
+        Derive terrain attribute for DEM opened with geoutils.Raster using RichDEM.
+        """
+        rda = raster_to_rda(rst)
+        terrattr = rd.TerrainAttribute(rda, attrib=attribute)
+        terrattr[terrattr == terrattr.no_data] = np.nan
+        return np.array(terrattr)
+
+    return _get_terrainattr_richdem
+
+
+@pytest.fixture(scope="session")  # type: ignore
+def get_terrain_attribute_richdem(
+    get_terrainattr_richdem: Callable[[RasterType, str], NDArrayf]
+) -> Callable[[RasterType, Union[str, list[str]], bool, float, float, float], Union[RasterType, list[RasterType]]]:
+    def _get_terrain_attribute_richdem(
+        dem: RasterType,
+        attribute: Union[str, List[str]],
+        degrees: bool = True,
+        hillshade_altitude: float = 45.0,
+        hillshade_azimuth: float = 315.0,
+        hillshade_z_factor: float = 1.0,
+    ) -> Union[RasterType, List[RasterType]]:
+        """
+        Derive one or multiple terrain attributes from a DEM using RichDEM.
+        """
+        if isinstance(attribute, str):
+            attribute = [attribute]
+
+        if not isinstance(dem, gu.Raster):
+            raise ValueError("DEM must be a geoutils.Raster object.")
+
+        terrain_attributes = {}
+
+        # Check which products should be made to optimize the processing
+        make_aspect = any(attr in attribute for attr in ["aspect", "hillshade"])
+        make_slope = any(
+            attr in attribute
+            for attr in [
+                "slope",
+                "hillshade",
+                "planform_curvature",
+                "aspect",
+                "profile_curvature",
+                "maximum_curvature",
+            ]
+        )
+        make_hillshade = "hillshade" in attribute
+        make_curvature = "curvature" in attribute
+        make_planform_curvature = "planform_curvature" in attribute or "maximum_curvature" in attribute
+        make_profile_curvature = "profile_curvature" in attribute or "maximum_curvature" in attribute
+
+        if make_slope:
+            terrain_attributes["slope"] = get_terrainattr_richdem(dem, "slope_radians")
+
+        if make_aspect:
+            # The aspect of RichDEM is returned in degrees, we convert to radians to match the others
+            terrain_attributes["aspect"] = np.deg2rad(get_terrainattr_richdem(dem, "aspect"))
+            # For flat slopes, RichDEM returns a 90° aspect by default, while GDAL return a 180° aspect
+            # We stay consistent with GDAL
+            slope_tmp = get_terrainattr_richdem(dem, "slope_radians")
+            terrain_attributes["aspect"][slope_tmp == 0] = np.pi
+
+        if make_hillshade:
+            # If a different z-factor was given, slopemap with exaggerated gradients.
+            if hillshade_z_factor != 1.0:
+                slopemap = np.arctan(np.tan(terrain_attributes["slope"]) * hillshade_z_factor)
+            else:
+                slopemap = terrain_attributes["slope"]
+
+            azimuth_rad = np.deg2rad(360 - hillshade_azimuth)
+            altitude_rad = np.deg2rad(hillshade_altitude)
+
+            # The operation below yielded the closest hillshade to GDAL (multiplying by 255 did not work)
+            # As 0 is generally no data for this uint8, we add 1 and then 0.5 for the rounding to occur between
+            # 1 and 255
+            terrain_attributes["hillshade"] = np.clip(
+                1.5
+                + 254
+                * (
+                    np.sin(altitude_rad) * np.cos(slopemap)
+                    + np.cos(altitude_rad) * np.sin(slopemap) * np.sin(azimuth_rad - terrain_attributes["aspect"])
+                ),
+                0,
+                255,
+            ).astype("float32")
+
+        if make_curvature:
+            terrain_attributes["curvature"] = get_terrainattr_richdem(dem, "curvature")
+
+        if make_planform_curvature:
+            terrain_attributes["planform_curvature"] = get_terrainattr_richdem(dem, "planform_curvature")
+
+        if make_profile_curvature:
+            terrain_attributes["profile_curvature"] = get_terrainattr_richdem(dem, "profile_curvature")
+
+        # Convert the unit if wanted.
+        if degrees:
+            for attr in ["slope", "aspect"]:
+                if attr not in terrain_attributes:
+                    continue
+                terrain_attributes[attr] = np.rad2deg(terrain_attributes[attr])
+
+        output_attributes = [terrain_attributes[key].reshape(dem.shape) for key in attribute]
+
+        if isinstance(dem, gu.Raster):
+            output_attributes = [
+                gu.Raster.from_array(attr, transform=dem.transform, crs=dem.crs, nodata=-99999)
+                for attr in output_attributes
+            ]
+
+        return output_attributes if len(output_attributes) > 1 else output_attributes[0]
+
+    return _get_terrain_attribute_richdem

tests/test_terrain.py

@@ -183,7 +183,7 @@ def test_attribute_functions_against_gdaldem(self, attribute: str) -> None:
         "attribute",
         ["slope_Horn", "aspect_Horn", "hillshade_Horn", "curvature", "profile_curvature", "planform_curvature"],
     )  # type: ignore
-    def test_attribute_functions_against_richdem(self, attribute: str) -> None:
+    def test_attribute_functions_against_richdem(self, attribute: str, get_terrain_attribute_richdem) -> None:
         """
         Test that all attribute functions give the same results as those of RichDEM within a small tolerance.
 
@@ -202,12 +202,14 @@ def test_attribute_functions_against_richdem(self, attribute: str) -> None:
 
         # Functions for RichDEM wrapper methods
         functions_richdem = {
-            "slope_Horn": lambda dem: xdem.terrain.slope(dem, degrees=True, use_richdem=True),
-            "aspect_Horn": lambda dem: xdem.terrain.aspect(dem, degrees=True, use_richdem=True),
-            "hillshade_Horn": lambda dem: xdem.terrain.hillshade(dem, use_richdem=True),
-            "curvature": lambda dem: xdem.terrain.curvature(dem, use_richdem=True),
-            "profile_curvature": lambda dem: xdem.terrain.profile_curvature(dem, use_richdem=True),
-            "planform_curvature": lambda dem: xdem.terrain.planform_curvature(dem, use_richdem=True),
+            "slope_Horn": lambda dem: get_terrain_attribute_richdem(dem, attribute="slope", degrees=True),
+            "aspect_Horn": lambda dem: get_terrain_attribute_richdem(dem, attribute="aspect", degrees=True),
+            "hillshade_Horn": lambda dem: get_terrain_attribute_richdem(dem, attribute="hillshade"),
+            "curvature": lambda dem: get_terrain_attribute_richdem(dem, attribute="curvature"),
+            "profile_curvature": lambda dem: get_terrain_attribute_richdem(dem, attribute="profile_curvature"),
+            "planform_curvature": lambda dem: get_terrain_attribute_richdem(
+                dem, attribute="planform_curvature", degrees=True
+            ),
         }
 
         # Copy the DEM to ensure that the inter-test state is unchanged, and because the mask will be modified.
@@ -355,15 +357,6 @@ def test_get_terrain_attribute_errors(self) -> None:
         """Test the get_terrain_attribute function raises appropriate errors."""
 
         # Below, re.escape() is needed to match expressions that have special characters (e.g., parenthesis, bracket)
-        with pytest.raises(
-            ValueError,
-            match=re.escape("RichDEM can only compute the slope and aspect using the " "default method of Horn (1981)"),
-        ):
-            xdem.terrain.slope(self.dem, method="ZevenbergThorne", use_richdem=True)
-
-        with pytest.raises(ValueError, match="To derive RichDEM attributes, the DEM passed must be a Raster object"):
-            xdem.terrain.slope(self.dem.data, resolution=self.dem.res, use_richdem=True)
-
         with pytest.raises(
             ValueError,
             match=re.escape(

xdem/dem.py

@@ -375,69 +375,44 @@ def to_vcrs(
             )
 
     @copy_doc(terrain, remove_dem_res_params=True)
-    def slope(
-        self,
-        method: str = "Horn",
-        degrees: bool = True,
-        use_richdem: bool = False,
-    ) -> RasterType:
-        return terrain.slope(self, method=method, degrees=degrees, use_richdem=use_richdem)
+    def slope(self, method: str = "Horn", degrees: bool = True) -> RasterType:
+        return terrain.slope(self, method=method, degrees=degrees)
 
     @copy_doc(terrain, remove_dem_res_params=True)
     def aspect(
         self,
         method: str = "Horn",
         degrees: bool = True,
-        use_richdem: bool = False,
     ) -> RasterType:
 
-        return terrain.aspect(self, method=method, degrees=degrees, use_richdem=use_richdem)
+        return terrain.aspect(self, method=method, degrees=degrees)
 
     @copy_doc(terrain, remove_dem_res_params=True)
     def hillshade(
-        self,
-        method: str = "Horn",
-        azimuth: float = 315.0,
-        altitude: float = 45.0,
-        z_factor: float = 1.0,
-        use_richdem: bool = False,
+        self, method: str = "Horn", azimuth: float = 315.0, altitude: float = 45.0, z_factor: float = 1.0
     ) -> RasterType:
 
-        return terrain.hillshade(
-            self, method=method, azimuth=azimuth, altitude=altitude, z_factor=z_factor, use_richdem=use_richdem
-        )
+        return terrain.hillshade(self, method=method, azimuth=azimuth, altitude=altitude, z_factor=z_factor)
 
     @copy_doc(terrain, remove_dem_res_params=True)
-    def curvature(
-        self,
-        use_richdem: bool = False,
-    ) -> RasterType:
+    def curvature(self) -> RasterType:
 
-        return terrain.curvature(self, use_richdem=use_richdem)
+        return terrain.curvature(self)
 
     @copy_doc(terrain, remove_dem_res_params=True)
-    def planform_curvature(
-        self,
-        use_richdem: bool = False,
-    ) -> RasterType:
+    def planform_curvature(self) -> RasterType:
 
-        return terrain.planform_curvature(self, use_richdem=use_richdem)
+        return terrain.planform_curvature(self)
 
     @copy_doc(terrain, remove_dem_res_params=True)
-    def profile_curvature(
-        self,
-        use_richdem: bool = False,
-    ) -> RasterType:
+    def profile_curvature(self) -> RasterType:
 
-        return terrain.profile_curvature(self, use_richdem=use_richdem)
+        return terrain.profile_curvature(self)
 
     @copy_doc(terrain, remove_dem_res_params=True)
-    def maximum_curvature(
-        self,
-        use_richdem: bool = False,
-    ) -> RasterType:
+    def maximum_curvature(self) -> RasterType:
 
-        return terrain.maximum_curvature(self, use_richdem=use_richdem)
+        return terrain.maximum_curvature(self)
 
     @copy_doc(terrain, remove_dem_res_params=True)
     def topographic_position_index(self, window_size: int = 3) -> RasterType: