adapt terrain scripts for flow-aligned hillshade

src/pdemtools/_accessor.py

@@ -257,27 +257,29 @@ def terrain(
                 # calculate unidirectional hillshade
                 if hillshade_z_factor == 1:
                     hillshade_arr = hillshade(
-                        slope_arr,
-                        aspect_arr,
+                        np.rad2deg(slope_arr),
+                        np.rad2deg(aspect_arr),
                         hillshade_altitude,
                         hillshade_azimuth,
                         norm=True,
                     )
                 if hillshade_z_factor != 1:
                     hillshade_arr = hillshade(
-                        slope(
-                            p(
-                                self._obj.values * hillshade_z_factor,
-                                resolution,
-                                method,
-                            ),
-                            q(
-                                self._obj.values * hillshade_z_factor,
-                                resolution,
-                                method,
-                            ),
+                        np.rad2deg(
+                            slope(
+                                p(
+                                    self._obj.values * hillshade_z_factor,
+                                    resolution,
+                                    method,
+                                ),
+                                q(
+                                    self._obj.values * hillshade_z_factor,
+                                    resolution,
+                                    method,
+                                ),
+                            )
                         ),
-                        aspect_arr,
+                        np.rad2deg(aspect_arr),
                         hillshade_altitude,
                         hillshade_azimuth,
                         norm=True,
@@ -310,18 +312,45 @@ def terrain(
                     )
 
                 hs = (
-                    w_225 * hillshade(slope_input, aspect_arr, 60, 225, norm=False)
-                    + w_270 * hillshade(slope_input, aspect_arr, 60, 270, norm=False)
-                    + w_315 * hillshade(slope_input, aspect_arr, 60, 315, norm=False)
-                    + w_360 * hillshade(slope_input, aspect_arr, 60, 360, norm=False)
+                    w_225
+                    * hillshade(
+                        np.rad2deg(slope_input),
+                        np.rad2deg(aspect_arr),
+                        60,
+                        225,
+                        norm=False,
+                    )
+                    + w_270
+                    * hillshade(
+                        np.rad2deg(slope_input),
+                        np.rad2deg(aspect_arr),
+                        60,
+                        270,
+                        norm=False,
+                    )
+                    + w_315
+                    * hillshade(
+                        np.rad2deg(slope_input),
+                        np.rad2deg(aspect_arr),
+                        60,
+                        315,
+                        norm=False,
+                    )
+                    + w_360
+                    * hillshade(
+                        np.rad2deg(slope_input),
+                        np.rad2deg(aspect_arr),
+                        60,
+                        360,
+                        norm=False,
+                    )
                 ) / 2
 
                 del w_225, w_270, w_315, w_360
 
                 # get normalised hillshade with darkest points as 0
-                hillshade_arr = 1 - (hs - np.nanmin(hs)) / (
-                    np.nanmax(hs) - np.nanmin(hs)
-                )
+                # hillshade_arr = 1 - (hs - np.nanmin(hs)) / (
+                hillshade_arr = (hs - np.nanmin(hs)) / (np.nanmax(hs) - np.nanmin(hs))
                 del hs
 
         if "horizontal_curvature" in attribute:

src/pdemtools/_geomorphometry.py

@@ -273,23 +273,48 @@ def t(z, w, method):
 
 
 def slope(p_arr, q_arr):
+    """outputs in radians"""
     return np.arctan(np.sqrt(p_arr**2 + q_arr**2))
 
 
 def aspect(p_arr, q_arr):
-    return np.deg2rad(
-        -90 * (1 - np.sign(q_arr)) * (1 - np.abs(np.sign(p_arr)))
-        + 180 * (1 + np.sign(p_arr))
-        - (180 / np.pi)
-        * np.sign(p_arr)
-        * np.arccos(-q_arr / np.sqrt(p_arr**2 + q_arr**2))
-    )
+    """outputs in radians"""
+    return np.arctan2(p_arr, q_arr) + np.pi
+    # return np.deg2rad(
+    #     -90 * (1 - np.sign(q_arr)) * (1 - np.abs(np.sign(p_arr)))
+    #     + 180 * (1 + np.sign(p_arr))
+    #     - (180 / np.pi)
+    #     * np.sign(p_arr)
+    #     * np.arccos(-q_arr / np.sqrt(p_arr**2 + q_arr**2))
+    # )
 
 
 def hillshade(slope, aspect, altitude=45, azimuth=315, norm=True):
-    hs = np.cos(np.pi * 0.5 - aspect - azimuth) * np.sin(slope) * np.sin(
-        np.pi * 0.5 - altitude
-    ) + np.cos(slope) * np.cos(np.pi * 0.5 - altitude)
+    """accepts degrees"""
+
+    if slope.max() < 2 * np.pi:
+        raise Warning(
+            "Maximum slope value is < 2π (slope.max()), ensure input is degrees"
+        )
+    if aspect.max() < 2 * np.pi:
+        raise Warning(
+            "Maximum slope value is < 2π (aspect.max()), ensure input is degrees"
+        )
+
+    hs = 255.0 * (
+        (np.cos(np.deg2rad(altitude)) * np.cos(np.deg2rad(slope)))
+        + (
+            np.sin(np.deg2rad(altitude))
+            * np.sin(np.deg2rad(slope))
+            * np.cos(np.deg2rad(azimuth) - np.deg2rad(aspect))
+        )
+    ).astype("float32")
+
+    # hs = np.cos(np.pi * 0.5 - np.deg2rad(aspect) - np.deg2rad(azimuth)) * np.sin(
+    #     np.deg2rad(slope)
+    # ) * np.sin(np.pi * 0.5 - np.deg2rad(altitude)) + np.cos(np.deg2rad(slope)) * np.cos(
+    #     np.pi * 0.5 - np.deg2rad(altitude)
+    # )
 
     # return normalised values
     if norm == True: