add c-hand method and example

doc/index.rst

@@ -40,6 +40,7 @@ Examples
 .. nbgallery::
    examples/fim_workflow_combined
    examples/fim_workflow_indiv_steps
+   examples/c_hand_example.ipynb
 
 .. Getting Started
 .. ###############

src/pygeoflood/pygeoflood.py

@@ -11,8 +11,8 @@
 from os import PathLike
 from pathlib import Path
 from rasterio.features import rasterize, shapes
-from rasterio.transform import from_bounds
-from rasterio.warp import transform_bounds,  reproject, calculate_default_transform, Resampling
+from rasterio.transform import from_bounds, rowcol
+from rasterio.warp import transform_bounds, reproject, calculate_default_transform, Resampling
 
 
 class PyGeoFlood(object):
@@ -81,6 +81,7 @@ def dem_path(self, value: str | PathLike):
     fsm_dephier_path = t.path_property("fsm_dephier_path")
     fsm_labels_path = t.path_property("fsm_labels_path")
     fsm_flowdir_path = t.path_property("fsm_flowdir_path")
+    coastal_inundation_path = t.path_property("coastal_inundation_path")
 
 
 
@@ -176,6 +177,7 @@ def __init__(
         self.fsm_dephier_path = f"{default_prefix}_fsm_dephier.json"
         self.fsm_labels_path = f"{default_prefix}_fsm_labels.npy"
         self.fsm_flowdir_path = f"{default_prefix}_fsm_flowdir.npy"
+        self.coastal_inundation_path = f"{default_prefix}_coastal_inundation.tif"
 
 
 
@@ -2614,9 +2616,72 @@ def dict_to_instance(data):
             dst.write(water_depth, 1)
 
         print(f"FSM inundation saved to {fsm_inundation}")
-
 
+    @t.time_it
+    @t.use_config_defaults
+    def c_hand(
+        self,
+        ocean_coords: tuple[float, float] = None,
+        xy: bool = True,
+        gage_el: float = None,
+        custom_dem: str | PathLike = None,
+        custom_path: str | PathLike = None,
+    ):
+        """
+        Calculate coastal inundation with the c-HAND method.
 
+        Parameters
+        ----------
+        ocean_coords : `tuple`
+            Tuple of coordinates (x, y) of a DEM cell in the ocean. Coordinates should be in the same CRS as the DEM.
+        xy : `bool`, optional
+            If False, ocean_coords are (row, col) of an array rather than (x, y) in a CRS. Default is True.
+        gage_el : `float`, optional
+            Constant water surface elevation to apply to DEM. Units and vertical datum must match DEM.
+        custom_dem : `str`, `os.PathLike`, optional
+            Custom file path to input dem. If not provided default 
+            DEM is used. A custom_path is required when a custom_dem is provided. Intermediary files
+            will be saved to match the name of the input custom_dem and be within the default working folder.
+        custom_path :  `str`, `os.PathLike`, optional
+            Custom path to save the coastal inundation raster. If not provided, coastal inundation raster 
+            will be saved in the project directory with the default name.
+        """
+        if custom_dem is None:
+            dem=self.dem_path
+        else:
+            dem=custom_dem
+            if custom_path is None:
+                raise ValueError("A custom path is required when a custom DEM is provided.")
+
+        t.check_attributes([("DEM", dem)], "c_hand")
+
+        if ocean_coords is None:
+            raise ValueError("ocean_coords must be provided")
+        if gage_el is None:
+            raise ValueError("gage_el must be provided")
+
+        # read original DEM
+        dem, dem_profile = t.read_raster(dem)
+
+        if xy:
+            ocean_coords = rowcol(dem_profile["transform"], *ocean_coords)
+
+        coastal_inun = t.get_c_hand(dem, gage_el, ocean_coords)
+
+        # set file path
+        if custom_path is None:
+            output_coastal_inundation_path = self.coastal_inundation_path    
+        else:
+            output_coastal_inundation_path = f"{custom_path}.tif"
+
+        t.write_raster(
+            raster=coastal_inun,
+            profile=dem_profile,
+            file_path=output_coastal_inundation_path,
+        )
+
+        print(f"Coastal inundation raster written to {output_coastal_inundation_path}")
+
 
     @t.time_it
     def run_fim_workflow(self):

src/pygeoflood/tools.py

@@ -24,6 +24,7 @@
 from shapely.geometry import LineString, Point, shape
 from shapely.ops import linemerge, snap, split
 from skimage.graph import route_through_array
+from skimage.measure import label
 from whitebox import WhiteboxTools
 
 warnings.filterwarnings(
@@ -1785,3 +1786,26 @@ def get_inun(hand, seg_catch, df):
     )
     inun = np.where(catch_h_mapped > hand, catch_h_mapped - hand, np.nan)
     return inun
+
+def get_c_hand(dem, gage_el, opix):
+
+    # initialize array with nan values
+    inun = np.full(dem.shape, np.nan, dtype=np.float32)
+
+    # initial inun array: 0 if (DEM ≥ gage_el) else (gage_el - DEM)
+    inun = np.where(dem >= (gage_el), 0, (gage_el) - dem)
+
+    # masked inun array: 255 if inun > 0 else 0
+    inun_mask = np.where(inun == 0, 0, 255)
+
+    # label connected regions of inundation
+    regions = label(inun_mask)
+
+    # only keep region containing the ocean pixel
+    inun = np.where(regions == regions[opix], inun, 0)
+
+    # return masked array if fed one
+    if isinstance(dem, np.ma.MaskedArray):
+        inun = np.ma.masked_array(inun, dem.mask)
+
+    return inun