Merge pull request #104 from erikmannerfelt/hillshade

Add hillshade function.

tests/test_spatial_tools.py

@@ -4,8 +4,13 @@
     Erik S. Holmlund
 
 """
+import os
+import shutil
+import subprocess
+import tempfile
+import warnings
+
 import geoutils as gu
-import matplotlib.pyplot as plt
 import numpy as np
 import rasterio as rio
 
@@ -71,12 +76,14 @@ def test_stack_rasters(self):
         assert merged_bounds == stacked_dem.bounds
 
         # Check that reference works with input Raster
-        stacked_dem = xdem.spatial_tools.stack_rasters([self.dem1, self.dem2], reference=self.dem)
+        stacked_dem = xdem.spatial_tools.stack_rasters(
+            [self.dem1, self.dem2], reference=self.dem)
         assert self.dem.bounds == stacked_dem.bounds
 
         # Others than int or gu.Raster should raise a ValueError
         try:
-            stacked_dem = xdem.spatial_tools.stack_rasters([self.dem1, self.dem2], reference="a string")
+            stacked_dem = xdem.spatial_tools.stack_rasters(
+                [self.dem1, self.dem2], reference="a string")
         except ValueError as exception:
             if "reference should be" not in str(exception):
                 raise exception
@@ -88,7 +95,8 @@ def test_stack_rasters(self):
         assert stacked_dem.bounds != self.dem.bounds
 
         # This case should preserve original extent
-        stacked_dem2 = xdem.spatial_tools.stack_rasters([self.dem1, self.dem3], reference=self.dem, use_ref_bounds=True)
+        stacked_dem2 = xdem.spatial_tools.stack_rasters(
+            [self.dem1, self.dem3], reference=self.dem, use_ref_bounds=True)
         assert stacked_dem2.bounds == self.dem.bounds
 
     def test_merge_rasters(self):
@@ -103,5 +111,68 @@ def test_merge_rasters(self):
         assert np.abs(np.nanmean(diff)) < 0.0001
 
         # Check that reference works
-        merged_dem2 = xdem.spatial_tools.merge_rasters([self.dem1, self.dem2], reference=self.dem)
+        merged_dem2 = xdem.spatial_tools.merge_rasters(
+            [self.dem1, self.dem2], reference=self.dem)
         assert merged_dem2 == merged_dem
+
+
+def test_hillshade():
+    """Test the hillshade algorithm, partly by comparing it to the GDAL hillshade function."""
+    warnings.simplefilter("error")
+
+    def make_gdal_hillshade(filepath) -> np.ndarray:
+        # rasterio strongly recommends against importing gdal along rio, so this is done here instead.
+        from osgeo import gdal
+        temp_dir = tempfile.TemporaryDirectory()
+        temp_hillshade_path = os.path.join(temp_dir.name, "hillshade.tif")
+        # gdal_commands = ["gdaldem", "hillshade",
+        #                 filepath, temp_hillshade_path,
+        #                 "-az", "315", "-alt", "45"]
+        #subprocess.run(gdal_commands, check=True, stdout=subprocess.PIPE)
+        gdal.DEMProcessing(
+            destName=temp_hillshade_path,
+            srcDS=filepath,
+            processing="hillshade",
+            options=gdal.DEMProcessingOptions(azimuth=315, altitude=45)
+        )
+
+        data = gu.Raster(temp_hillshade_path).data
+        temp_dir.cleanup()
+        return data
+
+    filepath = xdem.examples.FILEPATHS["longyearbyen_ref_dem"]
+    dem = xdem.DEM(filepath)
+
+    xdem_hillshade = xdem.spatial_tools.hillshade(dem.data, resolution=dem.res)
+    gdal_hillshade = make_gdal_hillshade(filepath)
+    diff = gdal_hillshade - xdem_hillshade
+
+    # Check that the xdem and gdal hillshades are relatively similar.
+    assert np.mean(diff) < 5
+    assert xdem.spatial_tools.nmad(diff.filled(np.nan)) < 5
+
+    # Try giving the hillshade invalid arguments.
+    try:
+        xdem.spatial_tools.hillshade(dem.data, dem.res, azimuth=361)
+    except ValueError as exception:
+        if "Azimuth must be a value between 0 and 360" not in str(exception):
+            raise exception
+    try:
+        xdem.spatial_tools.hillshade(dem.data, dem.res, altitude=91)
+    except ValueError as exception:
+        if "Altitude must be a value between 0 and 90" not in str(exception):
+            raise exception
+
+    try:
+        xdem.spatial_tools.hillshade(dem.data, dem.res, z_factor=np.inf)
+    except ValueError as exception:
+        if "z_factor must be a non-negative finite value" not in str(exception):
+            raise exception
+
+    # Introduce some nans
+    dem.data.mask = np.zeros_like(dem.data, dtype=bool)
+    dem.data.mask.ravel()[np.random.choice(
+        dem.data.size, 50000, replace=False)] = True
+
+    # Make sure that this doesn't create weird division warnings.
+    xdem.spatial_tools.hillshade(dem.data, dem.res)

xdem/spatial_tools.py

@@ -316,3 +316,73 @@ def merge_rasters(rasters: list[gu.georaster.Raster], reference: Union[int, gu.R
     )
 
     return merged_raster
+
+
+def hillshade(dem: Union[np.ndarray, np.ma.masked_array], resolution: Union[float, tuple[float, float]],
+              azimuth: float = 315.0, altitude: float = 45.0, z_factor: float = 1.0) -> np.ndarray:
+    """
+    Generate a hillshade from the given DEM.
+
+    :param dem: The input DEM to calculate the hillshade from.
+    :param resolution: One or two values specifying the resolution of the DEM.
+    :param azimuth: The azimuth in degrees (0-360°) going clockwise, starting from north.
+    :param altitude: The altitude in degrees (0-90°). 90° is straight from above.
+    :param z_factor: Vertical exaggeration factor.
+
+    :raises AssertionError: If the given DEM is not a 2D array.
+    :raises ValueError: If invalid argument types or ranges were given.
+
+    :returns: A hillshade with the dtype "float32" with value ranges of 0-255.
+    """
+    # Extract the DEM and mask
+    dem_values, mask = get_array_and_mask(dem.squeeze())
+    # The above is not guaranteed to copy the data, so this needs to be done first.
+    demc = dem_values.copy()
+
+    # Validate the inputs.
+    assert len(demc.shape) == 2, f"Expected a 2D array. Got shape: {dem.shape}"
+    if (azimuth < 0.0) or (azimuth > 360.0):
+        raise ValueError(f"Azimuth must be a value between 0 and 360 degrees (given value: {azimuth})")
+    if (altitude < 0.0) or (altitude > 90):
+        raise ValueError("Altitude must be a value between 0 and 90 degress (given value: {altitude})")
+    if (z_factor < 0.0) or not np.isfinite(z_factor):
+        raise ValueError(f"z_factor must be a non-negative finite value (given value: {z_factor})")
+
+    # Fill the nonfinite values with the median (or maybe just 0?) to not interfere with the gradient analysis.
+    demc[~np.isfinite(demc)] = np.nanmedian(demc)
+
+    # Multiply the DEM with the z_factor to increase the apparent height.
+    demc *= z_factor
+
+    # Parse the resolution argument. If it's subscriptable, it's assumed to be [X, Y] resolution.
+    try:
+        resolution[0]  # type: ignore
+    # If that fails, it's assumed to be the X&Y resolution.
+    except TypeError as exception:
+        if "not subscriptable" not in str(exception):
+            raise exception
+        resolution = (resolution,) * 2  # type: ignore
+
+    # Calculate the gradient of each pixel.
+    x_gradient, y_gradient = np.gradient(demc)
+    # Normalize by the radius of the resolution to make it resolution variant.
+    x_gradient /= resolution[0] * 0.5  # type: ignore
+    y_gradient /= resolution[1] * 0.5  # type: ignore
+
+    azimuth_rad = np.deg2rad(360 - azimuth)
+    altitude_rad = np.deg2rad(altitude)
+
+    # Calculate slope and aspect maps.
+    slope = np.pi / 2.0 - np.arctan(np.sqrt(x_gradient ** 2 + y_gradient ** 2))
+    aspect = np.arctan2(-x_gradient, y_gradient)
+
+    # Create a hillshade from these products.
+    shaded = np.sin(altitude_rad) * np.sin(slope) + np.cos(altitude_rad) * \
+        np.cos(slope) * np.cos((azimuth_rad - np.pi / 2.0) - aspect)
+
+    # Set (potential) masked out values to nan
+    shaded[mask] = np.nan
+
+    # Return the hillshade, scaled to uint8 ranges.
+    # The output is scaled by "(x + 0.6) / 1.84" to make it more similar to GDAL.
+    return np.clip(255 * (shaded + 0.6) / 1.84, 0, 255).astype("float32")