diff --git a/CHANGELOG.md b/CHANGELOG.md index a191d4d4..4b499a2d 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/prototype/water-extent-map-on-demand.ipynb b/prototype/water-extent-map-on-demand.ipynb index 3ac41675..110f43af 100644 --- a/prototype/water-extent-map-on-demand.ipynb +++ b/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", "}" ] diff --git a/src/asf_tools/hydrosar/flood_map.py b/src/asf_tools/hydrosar/flood_map.py index 9cc0ce43..7bbf189e 100644 --- a/src/asf_tools/hydrosar/flood_map.py +++ b/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. @@ -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}")