Merge pull request #211 from ASFHyP3/flood-wat

Some updates to more closely align with HydroSAR notebook

CHANGELOG.md

@@ -10,17 +10,21 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 ## [0.6.0]
 
 ## Added
-* You can choose whether an `fmi`  or `ts` minimization metric is used for the `asf_tools.hydrosar.flood_map.iterative` method:
-  * the `flood_map` console script entrypoint now accepts a `--min-metric` argument
-  * the  `asf_tools.hydrosar.floopd_map.make_flood_map` function now accepts a `min_metric` keyword argument
+* You can choose whether the `ts` (threat score; default) or `fmi` (Fowlkes-Mallows index) minimization metric is used for the flood mapping iterative estimator:
+  * the `flood_map` console script entrypoint now accepts a `--minimization-metric` argument
+  * the  `asf_tools.hydrosar.floopd_map.make_flood_map` function now accepts a `minimization_metric` keyword argument
+* The flood mapping iterative estimator will ignore waterbodies smaller than a minimum number of pixels (default = 0)
+  * the `flood_map` console script entrypoint now accepts a `--iterative-min-size` argument
+  * the  `asf_tools.hydrosar.floopd_map.make_flood_map` function now accepts a `iterative_min_size` keyword argument
 
 ### Changed
 * The HydroSAR code (`flood_map`, `water_map`, and `hand`) in `asf_tools` has been isolated to a `asf_tools.hydrosar` sub-package 
-* The `asf_tools.hydrosar.flood_map.iterative` method now generates an initial guess using the `nmad` method and then runs with a maximum step size of 3 instead of the default 0.5.
-* the known water threshold used to determine perennial water when creating flood maps will be calculated by default, or if `NaN`, using `asf_tools.hydrosar.flood_map.get_pw_threshold`
+* The `asf_tools.hydrosar.flood_map.iterative` estimator now runs with a maximum step size of 3 instead of the default 0.5.
+* The `asf_tools.hydrosar.flood_map.iterative` estimator now uses the mean of the iterative bounds at the initial guess.
+* the known water threshold used to determine perennial water when creating flood maps will be calculated `asf_tools.hydrosar.flood_map.get_pw_threshold` if not provided
 * `get_epsg_code` and `epsg_to_wkt` have been moved from`asf_tools.composite` to `asf_tools.util`
 * `read_as_array` and `write_cog` have been moved from`asf_tools.composite` to `asf_tools.raster`
-* * `get_coordinates` has been moved from`asf_tools.flood_map` to `asf_tools.util`
+* `get_coordinates` has been moved from`asf_tools.flood_map` to `asf_tools.util`
 
 ### Fixed
 * Reverted the special handling of nan values introduced in v0.5.2, now that GDAL v3.7.0 has been released.

prototype/water-extent-map-on-demand.ipynb

@@ -213,7 +213,7 @@
     "    'job_parameters': {\n",
     "        # RTC options\n",
     "        'resolution': 30,               # Desired output pixel spacing in meters\n",
-    "        'speckle_filter': True,         # Apply an Enhanced Lee speckle filter\n",
+    "        'speckle_filter': False,        # Apply an Enhanced Lee speckle filter\n",
     "        # Water extent options\n",
     "        'max_vv_threshold': -15.5,      # Maximum threshold value to use for VV polarized raster in decibels (dB)\n",
     "        'max_vh_threshold': -23.0,      # Maximum threshold value to use for VH polarized raster in decibels (dB)\n",
@@ -223,9 +223,11 @@
     "        # Flood  map options\n",
     "        'flood_depth_estimator': None,  # iterative', 'logstat', 'nmad', 'numpy' or None; flood maps won't be calculated if None\n",
     "        'water_level_sigma': 3.0,       # Standard deviation to estimate max water height for each object. Ignored when flood_depth_estimator is None\n",
-    "        'known_water_threshold': 30.0,  # Threshold for extracting known water area in percent. Ignored when flood_depth_estimator is None\n",
-    "        'iterative_min': 0,             # Minimum bound used for iterative method. Ignored when flood_depth_estimator is None\n",
-    "        'iterative_max': 15,            # Maximum bound used for iterative method. Ignored when flood_depth_estimator is None\n",
+    "        'known_water_threshold': None,  # Threshold for extracting known water area in percent. Computed when None and ignored when flood_depth_estimator is None\n",
+    "        'iterative_min': 0,             # Minimum bound used for iterative estimator. Ignored when flood_depth_estimator is None\n",
+    "        'iterative_max': 15,            # Maximum bound used for iterative estimator. Ignored when flood_depth_estimator is None\n",
+    "        'iterative_min_size': 0,        # Minimum size of a connected waterbody in pixels for calculating flood depths with the iterative estimator\n",
+    "        'minimization_metric': 'ts',    # Minimization method when using the iterative estimator (only); Fowlkes-Mallows index (fmi) or a threat score (ts)\n",
     "    }\n",
     "}"
    ]

src/asf_tools/hydrosar/flood_map.py

@@ -57,21 +57,22 @@ def get_waterbody(input_info: dict, threshold: Optional[float] = None) -> np.arr
     return water_array > threshold
 
 
-def iterative(hand: np.array, extent: np.array, x0: float = 7.5, water_levels: np.array = range(15),
-              minimization_metric: str = 'fmi'):
+def iterative(hand: np.array, extent: np.array, water_levels: np.array = np.arange(15),
+              minimization_metric: str = 'ts'):
     def get_confusion_matrix(w):
         iterative_flood_extent = hand < w  # w=water level
         tp = np.nansum(np.logical_and(iterative_flood_extent == 1, extent == 1))  # true positive
+        tn = np.nansum(np.logical_and(iterative_flood_extent == 0, extent == 0))  # true negative
         fp = np.nansum(np.logical_and(iterative_flood_extent == 1, extent == 0))  # False positive
         fn = np.nansum(np.logical_and(iterative_flood_extent == 0, extent == 1))  # False negative
-        return tp, fp, fn
+        return tp, tn, fp, fn
 
     def _goal_ts(w):
-        tp, fp, fn = get_confusion_matrix(w)
-        return 1 - tp / (tp + fp + fn)  # threat score #we will minimize goal func, hence 1-threat_score.
+        tp, _, fp, fn = get_confusion_matrix(w)
+        return 1 - tp / (tp + fp + fn)  # threat score -- we will minimize goal func, hence `1 - threat_score`.
 
     def _goal_fmi(w):
-        tp, fp, fn = get_confusion_matrix(w)
+        tp, _, fp, fn = get_confusion_matrix(w)
         return 1 - np.sqrt((tp/(tp+fp))*(tp/(tp+fn)))
 
     class MyBounds(object):
@@ -84,10 +85,14 @@ def __call__(self, **kwargs):
             tmax = bool(np.all(x <= self.xmax))
             tmin = bool(np.all(x >= self.xmin))
             return tmax and tmin
+
     bounds = MyBounds()
     MINIMIZATION_FUNCTION = {'fmi': _goal_fmi, 'ts': _goal_ts}
-    opt_res = optimize.basinhopping(MINIMIZATION_FUNCTION[minimization_metric], x0, niter=10000, niter_success=100,
-                                    accept_test=bounds, stepsize=3)
+    with warnings.catch_warnings():
+        warnings.filterwarnings('ignore', category=RuntimeWarning)
+        opt_res = optimize.basinhopping(MINIMIZATION_FUNCTION[minimization_metric], np.mean(water_levels),
+                                        niter=10000, niter_success=100, accept_test=bounds, stepsize=3)
+
     if opt_res.message[0] == 'success condition satisfied' \
             or opt_res.message[0] == 'requested number of basinhopping iterations completed successfully':
         return opt_res.x[0]
@@ -111,19 +116,24 @@ def logstat(data: np.ndarray, func: Callable = np.nanstd) -> Union[np.ndarray, f
     return np.exp(st)
 
 
-def estimate_flood_depth(label, hand, flood_labels, estimator='iterative', water_level_sigma=3.,
-                         iterative_bounds=(0, 15), minimization_metric: str = 'fmi'):
+def estimate_flood_depth(label: int, hand: np.ndarray, flood_labels: np.ndarray, estimator: str = 'iterative',
+                         water_level_sigma: float = 3., iterative_bounds: Tuple[int, int] = (0, 15),
+                         iterative_min_size: int = 0, minimization_metric: str = 'ts') -> float:
     with warnings.catch_warnings():
         warnings.filterwarnings('ignore', r'Mean of empty slice')
 
-        if estimator.lower() == "iterative" or estimator.lower() == "nmad":
+        if estimator.lower() == "iterative":
+            if (flood_labels == label).sum() < iterative_min_size:
+                return np.nan
+
+            water_levels = np.arange(*iterative_bounds)
+            return iterative(hand, flood_labels == label,
+                             water_levels=water_levels, minimization_metric=minimization_metric)
+
+        if estimator.lower() == "nmad":
             hand_mean = np.nanmean(hand[flood_labels == label])
-            hand_std = stats.median_abs_deviation(hand[flood_labels == label], scale='normal',
-                                                  nan_policy='omit')
-            if estimator.lower() == "iterative":
-                return iterative(hand, flood_labels == label,
-                                 x0=hand_mean + water_level_sigma * hand_std,
-                                 water_levels=iterative_bounds, minimization_metric=minimization_metric)
+            hand_std = stats.median_abs_deviation(hand[flood_labels == label], scale='normal', nan_policy='omit')
+
         if estimator.lower() == "numpy":
             hand_mean = np.nanmean(hand[flood_labels == label])
             hand_std = np.nanstd(hand[flood_labels == label])
@@ -144,7 +154,9 @@ def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path],
                    water_level_sigma: float = 3.,
                    known_water_threshold: Optional[float] = None,
                    iterative_bounds: Tuple[int, int] = (0, 15),
-                   minimization_metric: str = 'fmi'):
+                   iterative_min_size: int = 0,
+                   minimization_metric: str = 'ts',
+                   ):
     """Create a flood depth map from a surface water extent map.
 
     WARNING: This functionality is still under active development and the products
@@ -178,8 +190,12 @@ def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path],
         water_level_sigma: Max water height used in logstat, nmad, and numpy estimations
         known_water_threshold: Threshold for extracting the known water area in percent.
             If `None`, the threshold is calculated.
-        iterative_bounds: Bounds on basin-hopping algorithm used in iterative estimation
-        min_metric : Evaluation method to minimize in iterative estimation
+        iterative_bounds: Minimum and maximum bound on the flood depths calculated by the basin-hopping algorithm
+            used in the iterative estimator
+        iterative_min_size: Minimum size of a connected waterbody in pixels for calculating flood depths with the
+            iterative estimator. Waterbodies smaller than this wll be skipped.
+        minimization_metric: Evaluation method to minimize when using the iterative estimator.
+            Options include a Fowlkes-Mallows index (fmi) or a threat score (ts).
 
     References:
         Jean-Francios Pekel, Andrew Cottam, Noel Gorelik, Alan S. Belward. 2016. <https://doi:10.1038/nature20584>
@@ -206,7 +222,9 @@ def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path],
 
     labeled_flood_mask, num_labels = ndimage.label(flood_mask)
     object_slices = ndimage.find_objects(labeled_flood_mask)
-    log.info(f'Detected {num_labels} water bodies...')
+    log.info(f'Detected {num_labels} waterbodies...')
+    if estimator.lower() == 'iterative':
+        log.info(f'Skipping waterbodies less than {iterative_min_size} pixels.')
 
     flood_depth = np.zeros(flood_mask.shape)
 
@@ -218,9 +236,11 @@ def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path],
         flood_window = labeled_flood_mask[min0:max0, min1:max1]
         hand_window = hand_array[min0:max0, min1:max1]
 
-        water_height = estimate_flood_depth(ll, hand_window, flood_window, estimator=estimator,
-                                            water_level_sigma=water_level_sigma, iterative_bounds=iterative_bounds,
-                                            minimization_metric=minimization_metric)
+        water_height = estimate_flood_depth(
+            ll, hand_window, flood_window, estimator=estimator, water_level_sigma=water_level_sigma,
+            iterative_bounds=iterative_bounds, minimization_metric=minimization_metric,
+            iterative_min_size=iterative_min_size,
+        )
 
         flood_depth_window = flood_depth[min0:max0, min1:max1]
         flood_depth_window[flood_window == ll] = water_height - hand_window[flood_window == ll]
@@ -292,16 +312,23 @@ def _get_cli(interface: Literal['hyp3', 'main']) -> argparse.ArgumentParser:
                         help='Estimate max water height for each object.')
     parser.add_argument('--known-water-threshold', type=optional_float, default=None,
                         help='Threshold for extracting known water area in percent.'
-                             ' If `None`, threshold will be caluclated.')
-    parser.add_argument('--minimization-metric', type=str, default='fmi', choices=['fmi', 'ts'],
-                        help='Evaluation method to minimize in iterative estimation')
+                             ' If `None`, threshold will be calculated.')
+    parser.add_argument('--minimization-metric', type=str, default='ts', choices=['fmi', 'ts'],
+                        help='Evaluation method to minimize when using the iterative estimator. '
+                             'Options include a Fowlkes-Mallows index (fmi) or a threat score (ts).')
+    parser.add_argument('--iterative-min-size', type=int, default=0,
+                        help='Minimum size of a connected waterbody in pixels for calculating flood depths with the '
+                             'iterative estimator. Waterbodies smaller than this wll be skipped.')
 
     if interface == 'hyp3':
-        parser.add_argument('--iterative-min', type=int, default=0)
-        parser.add_argument('--iterative-max', type=int, default=15)
+        parser.add_argument('--iterative-min', type=int, default=0,
+                            help='Minimum bound on the flood depths calculated using the iterative estimator.')
+        parser.add_argument('--iterative-max', type=int, default=15,
+                            help='Maximum bound on the flood depths calculated using the iterative estimator.')
     elif interface == 'main':
-        # FIXME: add a help string
-        parser.add_argument('--iterative-bounds', type=int, nargs=2, default=[0, 15])
+        parser.add_argument('--iterative-bounds', type=int, nargs=2, default=[0, 15],
+                            help='Minimum and maximum bound on the flood depths calculated using the iterative '
+                                 'estimator.')
     else:
         raise NotImplementedError(f'Unknown interface: {interface}')
 
@@ -346,10 +373,10 @@ def hyp3():
         out_raster=flood_map_raster, vv_raster=vv_raster, water_raster=water_map_raster, hand_raster=hand_raster,
         estimator=args.estimator, water_level_sigma=args.water_level_sigma,
         known_water_threshold=args.known_water_threshold, iterative_bounds=(args.iterative_min, args.iterative_max),
-        minimization_metric=args.minimization_metric
+        iterative_min_size=args.iterative_min_size, minimization_metric=args.minimization_metric,
     )
 
-    log.info(f"Flood depth map created successfully: {flood_map_raster}")
+    log.info(f'Flood depth map created successfully: {flood_map_raster}')
 
     if args.bucket:
         output_zip = make_archive(base_name=product_name, format='zip', base_dir=product_name)
@@ -366,8 +393,11 @@ def main():
     logging.basicConfig(stream=sys.stdout, format='%(asctime)s - %(levelname)s - %(message)s', level=level)
     log.debug(' '.join(sys.argv))
 
-    make_flood_map(args.out_raster, args.vv_raster, args.water_extent_map, args.hand_raster,
-                   args.estimator, args.water_level_sigma, args.known_water_threshold, tuple(args.iterative_bounds),
-                   args.minimization_metric)
+    make_flood_map(
+        out_raster=args.out_raster, vv_raster=args.vv_raster, water_raster=args.water_extent_map,
+        hand_raster=args.hand_raster, estimator=args.estimator, water_level_sigma=args.water_level_sigma,
+        known_water_threshold=args.known_water_threshold, iterative_bounds=tuple(args.iterative_bounds),
+        iterative_min_size=args.iterative_min_size, minimization_metric=args.minimization_metric,
+    )
 
     log.info(f"Flood depth map created successfully: {args.out_raster}")