Merge pull request #10 from trchudley/flowaware_hillshade

Flowaware hillshade capability

docs/appendix/references.md

@@ -6,6 +6,8 @@ Howat, I., _et al._ (2022a). The Reference Elevation Model of Antarctica – Str
 
 Howat, I., _et al._ (2022b). The Reference Elevation Model of Antarctica – Mosaics, Version 2, _Harvard Dataverse_ https://doi.org/10.7910/DVN/EBW8UC
 
+MacGregror, J. A. _et al._ (2024). Geologic Provinces Beneath the Greenland Ice Sheet Constrained by Geophysical Data Synthesis. _Geophysical Research Letters_, 51, e2023GL107357. https://doi.org/10.1029/2023GL107357
+
 Mark, R. K. (1992). Multidirectional, oblique-weighted, shaded-relief image of the Island of Hawaii. _United States Geological Survey_. Open-File Report 92-422. https://doi.org/10.3133/ofr92422
 
 Morlighem, M. _et al._ (2022a). IceBridge BedMachine Greenland, Version 5 [Data Set]. _NASA National Snow and Ice Data Center Distributed Active Archive Center_. https://doi.org/10.5067/GMEVBWFLWA7X

docs/appendix/version_updates.md

@@ -2,6 +2,7 @@
 
 | Version | Date | Notes |
 | ------- | ---- | ----- |
+| 0.8.1 | June 2024 | Minor modification to hillshade script to allow for MacGregor _et al._ (2024) flow-aware hillshade. Supplementary notebook available in GitHub `notebooks` directory. |
 | 0.8.0 | May 2024 | Iceberg detection now corrects for sea surface; final additions to documentation before public conda/pip release.
 | 0.7.0 | March 2024 | Minor bug fixes and more detailed docstrings to support readthedocs | 
 | 0.6.0 | February 2024 | Added `data.bedrock_mask_from_vector()` function. |

notebooks/mosaic_and_terrain.ipynb

@@ -1052,7 +1052,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python (Geospatial)",
+   "display_name": "Python (geospatial)",
    "language": "python",
    "name": "geospatial"
   },
@@ -1066,7 +1066,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.10.9"
   }
  },
  "nbformat": 4,

src/pdemtools/__init__.py

@@ -6,6 +6,6 @@
 
 from ._accessor import DemAccessor
 
-__version__ = "0.8.0"
+__version__ = "0.8.1"
 
 __all__ = ["search", "DemAccessor"]

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,50 @@ 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: