ENH: Add support for aspect function (only in rasters). #29

CHANGES.md

@@ -3,9 +3,10 @@
 ## 1.44.0 (2024-mm-dd)
 
 - **BREAKING CHANGE**: Renaming of `rasters(_rio).write` argument `path` to `output_path` to avoid shadowing `sertit.path` module. Older argument is deprecated.
-- ENH: Use `odc.geo.xr.mask` instead of `rasterio.mask` to be dask-compatible in `rasters` (true for `paint` and `mask`)
-- ENH: Use `odc.geo.xr.xr_reproject` instead of `rioxarray.reproject_match` to be dask-compatible in `rasters.collocate`
-- ENH: Use `xarray-spatial` to be dask-compatible in `rasters.slope` and `rasters.hillshade`
+- **ENH: Use `odc.geo.xr.mask` instead of `rasterio.mask` to be dask-compatible in `rasters` (true for `paint` and `mask`)** ([#27](https://github.com/sertit/sertit-utils/issues/27))
+- **ENH: Use `odc.geo.xr.xr_reproject` instead of `rioxarray.reproject_match` to be dask-compatible in `rasters.collocate`** ([#27](https://github.com/sertit/sertit-utils/issues/27))
+- **ENH: Use `xarray-spatial` to be dask-compatible in `rasters.slope` and `rasters.hillshade`.** ([#27](https://github.com/sertit/sertit-utils/issues/27))
+- **ENH: Add support for `aspect` function (only in `rasters`).** ([#29](https://github.com/sertit/sertit-utils/issues/29))
 - FIX: Fix the ability to save COGs with any dtype with Dask, with the workaround described [here](https://github.com/opendatacube/odc-geo/issues/189#issuecomment-2513450481) (don't compute statistics for problematic dtypes)
 - FIX: Better separability of `dask` (it has its own module now): don't create a client if the user doesn't specify it (as it is not required anymore in `Lock`). This should remove the force-use of `dask`.
 - OPTIM: For arrays with same shape and CRS, replace only the coordinates of `other` by `ref`'s in `rasters.collocate`

CI/SCRIPTS/test_rasters.py

@@ -542,6 +542,7 @@ def test_vrt(tmp_path, raster_path):
 )
 def test_merge_different_crs(tmp_path):
     """Test merge_vrt (with different CRS) function"""
+    tmp_path = "/home/data/CI"
     # DIFFERENT CRS
     true_vrt_path = rasters_path().joinpath("merge_32-31.vrt")
     true_tif_path = rasters_path().joinpath("merge_32-31.tif")
@@ -760,18 +761,27 @@ def test_where():
 @dask_env
 def test_dem_fct(tmp_path):
     """Test DEM fct, i.e. slope and hillshade"""
+    tmp_path = "/home/data/CI"
+
     # Paths IN
     dem_path = rasters_path().joinpath("dem.tif")
+    aspect_path = rasters_path().joinpath("aspect.tif")
     hlsd_path = rasters_path().joinpath("hillshade.tif")
     slope_path = rasters_path().joinpath("slope.tif")
     slope_r_path = rasters_path().joinpath("slope_r.tif")
     slope_p_path = rasters_path().joinpath("slope_p.tif")
 
     # Path OUT
+    aspect_path_out = os.path.join(tmp_path, "aspect_out.tif")
     hlsd_path_out = os.path.join(tmp_path, "hillshade_out.tif")
-    slope_path_out = os.path.join(tmp_path, "slope.tif")
-    slope_r_path_out = os.path.join(tmp_path, "slope_r.tif")
-    slope_p_path_out = os.path.join(tmp_path, "slope_p.tif")
+    slope_path_out = os.path.join(tmp_path, "slope_out.tif")
+    slope_r_path_out = os.path.join(tmp_path, "slope_r_out.tif")
+    slope_p_path_out = os.path.join(tmp_path, "slope_p_out.tif")
+
+    # Aspect
+    aspect = rasters.aspect(dem_path)
+    rasters.write(aspect, aspect_path_out, dtype="float32")
+    ci.assert_raster_almost_equal(aspect_path, aspect_path_out, decimal=4)
 
     # Hillshade
     hlsd = rasters.hillshade(dem_path, 34.0, 45.2)

sertit/rasters.py

@@ -1879,7 +1879,7 @@ def hillshade(
     - `2 <http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=How%20Hillshade%20works>`_
 
     Args:
-        xds (AnyRasterType): Path to the raster, its dataset, its :code:`xarray` or a tuple containing its array and metadata
+        xds (AnyRasterType): Path to the DEM, its dataset, its :code:`xarray` or a tuple containing its array and metadata
         azimuth (float): Azimuth of the light, in degrees. 0 if it comes from the top of the raster, 90 from the east, ...
         zenith (float): Zenith angle in degrees
 
@@ -1951,7 +1951,7 @@ def slope(
     Goal: replace `gdaldem CLI <https://gdal.org/programs/gdaldem.html>`_
 
     Args:
-        xds (AnyRasterType): Path to the raster, its dataset, its :code:`xarray` or a tuple containing its array and metadata
+        xds (AnyRasterType): Path to the DEM, its dataset, its :code:`xarray` or a tuple containing its array and metadata
         in_pct (bool): Outputs slope in percents
         in_rad (bool): Outputs slope in radians. Not taken into account if :code:`in_pct == True`
 
@@ -1981,4 +1981,26 @@ def slope(
     return xds
 
 
-# TODO: add other dem-related functions like 'aspect'. Create a dedicated module?
+@any_raster_to_xr_ds
+@_3d_to_2d
+def aspect(xds: AnyRasterType, **kwargs) -> AnyXrDataStructure:
+    """
+    Compute the aspect of a DEM.
+
+    Args:
+        xds (AnyRasterType): Path to the DEM, its dataset, its :code:`xarray` or a tuple containing its array and metadata
+
+    Returns:
+        AnyXrDataStructure: Aspect
+    """
+    try:
+        from xrspatial import aspect
+
+        return aspect(xds, name=kwargs.get("name", "aspect"))
+    except ImportError:
+        raise NotImplementedError(
+            "'Aspect' cannot be computed when 'xarray-spatial' is not installed."
+        )
+
+
+# TODO: add other DEM-related functions like 'curvature', etc if needed. Create a dedicated module?