diff --git a/.github/workflows/static-analysis.yml b/.github/workflows/static-analysis.yml index 243c6583..1d08910e 100644 --- a/.github/workflows/static-analysis.yml +++ b/.github/workflows/static-analysis.yml @@ -3,25 +3,10 @@ name: Static analysis on: push jobs: - flake8: - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v4 - - - uses: actions/setup-python@v5 - with: - python-version: "3.10" - - - name: Install dependencies - run: | - python -m pip install --upgrade pip - python -m pip install flake8 flake8-import-order flake8-blind-except flake8-builtins - - - name: Lint with flake8 - run: | - flake8 --max-line-length=120 --import-order-style=pycharm --statistics \ - --application-import-names asf_tools ArcGIS-toolbox/ASF_Tools.pyt src/asf_tools - - call-secrets-analysis-workflow: + # Docs: https://github.com/ASFHyP3/actions uses: ASFHyP3/actions/.github/workflows/reusable-secrets-analysis.yml@v0.12.0 + + call-ruff-workflow: + # Docs: https://github.com/ASFHyP3/actions + uses: ASFHyP3/actions/.github/workflows/reusable-ruff.yml@v0.12.0 diff --git a/CHANGELOG.md b/CHANGELOG.md index 10287df3..1f53da5a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,9 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ### Fixed - The [`release`](.github/workflows/release.yml) Github Actions workflow now uses the `gh` CLI instead of the archived `repo-sync/pull-request` action. +### Changed +- The [`static-analysis`](.github/workflows/static-analysis.yml) Github Actions workflow now uses `ruff` rather than `flake8`. + ## [0.8.0] ### Removed diff --git a/environment.yml b/environment.yml index 1c39b83d..dc66800d 100644 --- a/environment.yml +++ b/environment.yml @@ -9,10 +9,7 @@ dependencies: # For packaging, and testing # - arcpy # windows only - python-build - - flake8 - - flake8-import-order - - flake8-blind-except - - flake8-builtins + - ruff - setuptools>=61 - setuptools_scm>=6.2 - pytest diff --git a/pyproject.toml b/pyproject.toml index da7ccc57..42a13fbf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -52,11 +52,8 @@ flood_map = "asf_tools.hydrosar.flood_map:hyp3" [project.optional-dependencies] develop = [ - "flake8", - "flake8-import-order", - "flake8-blind-except", - "flake8-builtins", "gdal-utils", + "ruff", "pytest", "pytest-cov", "pytest-console-scripts", diff --git a/ruff.toml b/ruff.toml new file mode 100644 index 00000000..2aed2d63 --- /dev/null +++ b/ruff.toml @@ -0,0 +1,6 @@ +exclude = ["prototype"] + +line-length = 120 + +[format] +quote-style = "single" diff --git a/src/asf_tools/__main__.py b/src/asf_tools/__main__.py index c62666ec..83dc9ae1 100644 --- a/src/asf_tools/__main__.py +++ b/src/asf_tools/__main__.py @@ -6,8 +6,10 @@ def main(): parser = argparse.ArgumentParser(prefix_chars='+', formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument( - '++process', choices=['water_map', 'flood_map'], default='water_map', - help='Select the HyP3 entrypoint to use' # HyP3 entrypoints are specified in `pyproject.toml` + '++process', + choices=['water_map', 'flood_map'], + default='water_map', + help='Select the HyP3 entrypoint to use', # HyP3 entrypoints are specified in `pyproject.toml` ) args, unknowns = parser.parse_known_args() @@ -15,9 +17,7 @@ def main(): (process_entry_point,) = set(entry_points(group='hyp3', name=args.process)) sys.argv = [args.process, *unknowns] - sys.exit( - process_entry_point.load()() - ) + sys.exit(process_entry_point.load()()) if __name__ == '__main__': diff --git a/src/asf_tools/aws.py b/src/asf_tools/aws.py index 9e47daab..a72ab525 100644 --- a/src/asf_tools/aws.py +++ b/src/asf_tools/aws.py @@ -10,14 +10,7 @@ def get_tag_set() -> dict: - tag_set = { - 'TagSet': [ - { - 'Key': 'file_type', - 'Value': 'product' - } - ] - } + tag_set = {'TagSet': [{'Key': 'file_type', 'Value': 'product'}]} return tag_set diff --git a/src/asf_tools/composite.py b/src/asf_tools/composite.py index 1cf42599..682cf295 100755 --- a/src/asf_tools/composite.py +++ b/src/asf_tools/composite.py @@ -123,19 +123,27 @@ def reproject_to_target(raster_info: dict, target_epsg_code: int, target_resolut log.info(f'Reprojecting {raster}') reprojected_raster = os.path.join(directory, os.path.basename(raster)) gdal.Warp( - reprojected_raster, raster, dstSRS=f'EPSG:{target_epsg_code}', - xRes=target_resolution, yRes=target_resolution, targetAlignedPixels=True + reprojected_raster, + raster, + dstSRS=f'EPSG:{target_epsg_code}', + xRes=target_resolution, + yRes=target_resolution, + targetAlignedPixels=True, ) area_raster = get_area_raster(raster) log.info(f'Reprojecting {area_raster}') reprojected_area_raster = os.path.join(directory, os.path.basename(area_raster)) gdal.Warp( - reprojected_area_raster, area_raster, dstSRS=f'EPSG:{target_epsg_code}', - xRes=target_resolution, yRes=target_resolution, targetAlignedPixels=True + reprojected_area_raster, + area_raster, + dstSRS=f'EPSG:{target_epsg_code}', + xRes=target_resolution, + yRes=target_resolution, + targetAlignedPixels=True, ) - target_raster_info[reprojected_raster] = gdal.Info(reprojected_raster, format='json') + target_raster_info[reprojected_raster] = gdal.Info(reprojected_raster, format='json') else: log.info(f'No need to reproject {raster}') target_raster_info[raster] = info @@ -173,8 +181,12 @@ def make_composite(out_name: str, rasters: List[str], resolution: float = None): # resample rasters to maximum resolution & common UTM zone with TemporaryDirectory(prefix='reprojected_') as temp_dir: - raster_info = reproject_to_target(raster_info, target_epsg_code=target_epsg_code, target_resolution=resolution, - directory=temp_dir) + raster_info = reproject_to_target( + raster_info, + target_epsg_code=target_epsg_code, + target_resolution=resolution, + directory=temp_dir, + ) # Get extent of union of all images full_ul, full_lr, full_trans = get_full_extent(raster_info) @@ -188,8 +200,10 @@ def make_composite(out_name: str, rasters: List[str], resolution: float = None): for raster, info in raster_info.items(): log.info(f'Processing raster {raster}') - log.debug(f"Raster upper left: {info['cornerCoordinates']['upperLeft']}; " - f"lower right: {info['cornerCoordinates']['lowerRight']}") + log.debug( + f"Raster upper left: {info['cornerCoordinates']['upperLeft']}; " + f"lower right: {info['cornerCoordinates']['lowerRight']}" + ) values = read_as_array(raster) @@ -224,7 +238,13 @@ def make_composite(out_name: str, rasters: List[str], resolution: float = None): out_raster = write_cog(f'{out_name}.tif', outputs, full_trans, target_epsg_code, nodata_value=0) del outputs - out_counts_raster = write_cog(f'{out_name}_counts.tif', counts, full_trans, target_epsg_code, dtype=gdal.GDT_Int16) + out_counts_raster = write_cog( + f'{out_name}_counts.tif', + counts, + full_trans, + target_epsg_code, + dtype=gdal.GDT_Int16, + ) del counts return out_raster, out_counts_raster @@ -237,14 +257,21 @@ def main(): ) parser.add_argument('out_name', help='Base name of output composite GeoTIFF (without extension)') parser.add_argument('rasters', nargs='+', help='Sentinel-1 GeoTIFF rasters to composite') - parser.add_argument('-r', '--resolution', type=float, - help='Desired output resolution in meters ' - '(default is the max resolution of all the input files)') + parser.add_argument( + '-r', + '--resolution', + type=float, + help='Desired output resolution in meters ' '(default is the max resolution of all the input files)', + ) parser.add_argument('-v', '--verbose', action='store_true', help='Turn on verbose logging') args = parser.parse_args() level = logging.DEBUG if args.verbose else logging.INFO - logging.basicConfig(stream=sys.stdout, format='%(asctime)s - %(levelname)s - %(message)s', level=level) + logging.basicConfig( + stream=sys.stdout, + format='%(asctime)s - %(levelname)s - %(message)s', + level=level, + ) log.debug(' '.join(sys.argv)) log.info(f'Creating a composite of {len(args.rasters)} rasters') diff --git a/src/asf_tools/dem.py b/src/asf_tools/dem.py index 8e2a4fbb..4991fa98 100644 --- a/src/asf_tools/dem.py +++ b/src/asf_tools/dem.py @@ -1,4 +1,5 @@ """Prepare a Copernicus GLO-30 DEM virtual raster (VRT) covering a given geometry""" + from pathlib import Path from typing import Union @@ -31,7 +32,7 @@ def prepare_dem_vrt(vrt: Union[str, Path], geometry: Union[ogr.Geometry, BaseGeo geometry = ogr.CreateGeometryFromWkb(geometry.wkb) min_lon, max_lon, _, _ = geometry.GetEnvelope() - if min_lon < -160. and max_lon > 160.: + if min_lon < -160.0 and max_lon > 160.0: raise ValueError(f'asf_tools does not currently support geometries that cross the antimeridian: {geometry}') tile_features = vector.get_features(DEM_GEOJSON) diff --git a/src/asf_tools/hydrosar/__init__.py b/src/asf_tools/hydrosar/__init__.py index eb1b22fc..3caaaa1d 100644 --- a/src/asf_tools/hydrosar/__init__.py +++ b/src/asf_tools/hydrosar/__init__.py @@ -1,7 +1,6 @@ from warnings import warn -HYDROSAR_MOVE_WARNING = \ - """ +HYDROSAR_MOVE_WARNING = """ --------------------------------------------------------------------------- The HydroSAR codes (`flood_map`, `water_map` and `hand` modules) are being moved to the HydroSAR project repository: @@ -12,4 +11,4 @@ ---------------------------------------------------------------------------- """ -warn(HYDROSAR_MOVE_WARNING, category=FutureWarning) +warn(HYDROSAR_MOVE_WARNING, category=FutureWarning) diff --git a/src/asf_tools/hydrosar/flood_map.py b/src/asf_tools/hydrosar/flood_map.py index 7bbf189e..c514ca60 100644 --- a/src/asf_tools/hydrosar/flood_map.py +++ b/src/asf_tools/hydrosar/flood_map.py @@ -46,9 +46,16 @@ def get_waterbody(input_info: dict, threshold: Optional[float] = None) -> np.arr water_extent_vrt = data_dir / 'water_extent.vrt' # All Perennial Flood Data with tempfile.NamedTemporaryFile() as water_extent_file: - gdal.Warp(water_extent_file.name, str(water_extent_vrt), dstSRS=f'EPSG:{epsg}', - outputBounds=[west, south, east, north], - width=width, height=height, resampleAlg='nearest', format='GTiff') + gdal.Warp( + water_extent_file.name, + str(water_extent_vrt), + dstSRS=f'EPSG:{epsg}', + outputBounds=[west, south, east, north], + width=width, + height=height, + resampleAlg='nearest', + format='GTiff', + ) water_array = gdal.Open(water_extent_file.name, gdal.GA_ReadOnly).ReadAsArray() if threshold is None: @@ -57,8 +64,12 @@ def get_waterbody(input_info: dict, threshold: Optional[float] = None) -> np.arr return water_array > threshold -def iterative(hand: np.array, extent: np.array, water_levels: np.array = np.arange(15), - minimization_metric: str = 'ts'): +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 @@ -73,7 +84,7 @@ def _goal_ts(w): def _goal_fmi(w): tp, _, fp, fn = get_confusion_matrix(w) - return 1 - np.sqrt((tp/(tp+fp))*(tp/(tp+fn))) + return 1 - np.sqrt((tp / (tp + fp)) * (tp / (tp + fn))) class MyBounds(object): def __init__(self, xmax=max(water_levels), xmin=min(water_levels)): @@ -81,7 +92,7 @@ def __init__(self, xmax=max(water_levels), xmin=min(water_levels)): self.xmin = np.array(xmin) def __call__(self, **kwargs): - x = kwargs["x_new"] + x = kwargs['x_new'] tmax = bool(np.all(x <= self.xmax)) tmin = bool(np.all(x >= self.xmin)) return tmax and tmin @@ -90,25 +101,33 @@ def __call__(self, **kwargs): MINIMIZATION_FUNCTION = {'fmi': _goal_fmi, 'ts': _goal_ts} 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) + 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': + 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] else: return np.inf # set as inf to mark unstable solution def logstat(data: np.ndarray, func: Callable = np.nanstd) -> Union[np.ndarray, float]: - """ Calculate a function in logarithmic scale and return in linear scale. - INF values inside the data array are set to nan. - - Args: - data: array of data - func: statistical function to calculate in logarithmic scale - Returns: - statistic: statistic of data in linear scale + """Calculate a function in logarithmic scale and return in linear scale. + INF values inside the data array are set to nan. + + Args: + data: array of data + func: statistical function to calculate in logarithmic scale + Returns: + statistic: statistic of data in linear scale """ ld = np.log(data) ld[np.isinf(ld)] = np.nan @@ -116,29 +135,40 @@ def logstat(data: np.ndarray, func: Callable = np.nanstd) -> Union[np.ndarray, f return np.exp(st) -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: +def estimate_flood_depth( + label: int, + hand: np.ndarray, + flood_labels: np.ndarray, + estimator: str = 'iterative', + water_level_sigma: float = 3.0, + 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": + 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": + 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() == "numpy": + if estimator.lower() == 'numpy': hand_mean = np.nanmean(hand[flood_labels == label]) hand_std = np.nanstd(hand[flood_labels == label]) - elif estimator.lower() == "logstat": + elif estimator.lower() == 'logstat': hand_mean = logstat(hand[flood_labels == label], func=np.nanmean) hand_std = logstat(hand[flood_labels == label]) @@ -148,15 +178,18 @@ def estimate_flood_depth(label: int, hand: np.ndarray, flood_labels: np.ndarray, return hand_mean + water_level_sigma * hand_std -def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], - water_raster: Union[str, Path], hand_raster: Union[str, Path], - estimator: str = 'iterative', - water_level_sigma: float = 3., - known_water_threshold: Optional[float] = None, - iterative_bounds: Tuple[int, int] = (0, 15), - iterative_min_size: int = 0, - minimization_metric: str = 'ts', - ): +def make_flood_map( + out_raster: Union[str, Path], + vv_raster: Union[str, Path], + water_raster: Union[str, Path], + hand_raster: Union[str, Path], + estimator: str = 'iterative', + water_level_sigma: float = 3.0, + known_water_threshold: Optional[float] = None, + iterative_bounds: Tuple[int, int] = (0, 15), + 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 @@ -208,8 +241,14 @@ def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], log.info('Fetching perennial flood data.') known_water_mask = get_waterbody(info, threshold=known_water_threshold) - write_cog(str(out_raster).replace('.tif', f'_{estimator}_PW.tif'), known_water_mask, transform=geotransform, - epsg_code=epsg, dtype=gdal.GDT_Byte, nodata_value=False) + write_cog( + str(out_raster).replace('.tif', f'_{estimator}_PW.tif'), + known_water_mask, + transform=geotransform, + epsg_code=epsg, + dtype=gdal.GDT_Byte, + nodata_value=False, + ) water_map = gdal.Open(str(water_raster)).ReadAsArray() flood_mask = np.logical_or(water_map, known_water_mask) @@ -237,8 +276,13 @@ def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], 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, + 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, ) @@ -254,16 +298,34 @@ def make_flood_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], flood_mask_byte = flood_mask.astype(np.uint8) flood_mask_byte[padding_mask] = floodmask_nodata - write_cog(str(out_raster).replace('.tif', f'_{estimator}_WaterDepth.tif'), flood_depth, transform=geotransform, - epsg_code=epsg, dtype=gdal.GDT_Float64, nodata_value=nodata) - write_cog(str(out_raster).replace('.tif', f'_{estimator}_FloodMask.tif'), flood_mask_byte, transform=geotransform, - epsg_code=epsg, dtype=gdal.GDT_Byte, nodata_value=floodmask_nodata) + write_cog( + str(out_raster).replace('.tif', f'_{estimator}_WaterDepth.tif'), + flood_depth, + transform=geotransform, + epsg_code=epsg, + dtype=gdal.GDT_Float64, + nodata_value=nodata, + ) + write_cog( + str(out_raster).replace('.tif', f'_{estimator}_FloodMask.tif'), + flood_mask_byte, + transform=geotransform, + epsg_code=epsg, + dtype=gdal.GDT_Byte, + nodata_value=floodmask_nodata, + ) flood_mask[known_water_mask] = False flood_depth[np.logical_not(flood_mask)] = 0 flood_depth[padding_mask] = nodata - write_cog(str(out_raster).replace('.tif', f'_{estimator}_FloodDepth.tif'), flood_depth, transform=geotransform, - epsg_code=epsg, dtype=gdal.GDT_Float64, nodata_value=nodata) + write_cog( + str(out_raster).replace('.tif', f'_{estimator}_FloodDepth.tif'), + flood_depth, + transform=geotransform, + epsg_code=epsg, + dtype=gdal.GDT_Float64, + nodata_value=nodata, + ) def optional_str(value: str) -> Optional[str]: @@ -279,56 +341,87 @@ def optional_float(value: str) -> Optional[float]: def _get_cli(interface: Literal['hyp3', 'main']) -> argparse.ArgumentParser: - parser = argparse.ArgumentParser( - description=__doc__, - formatter_class=argparse.ArgumentDefaultsHelpFormatter - ) + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) available_estimators = ['iterative', 'logstat', 'nmad', 'numpy'] estimator_help = 'Flood depth estimation approach.' if interface == 'hyp3': parser.add_argument('--bucket') parser.add_argument('--bucket-prefix', default='') - parser.add_argument('--wm-raster', - help='Water map GeoTIFF raster, with suffix `_WM.tif`.') + parser.add_argument('--wm-raster', help='Water map GeoTIFF raster, with suffix `_WM.tif`.') available_estimators.append(None) estimator_help += ' If `None`, flood depth will not be calculated.' elif interface == 'main': - parser.add_argument('out_raster', - help='File to which flood depth map will be saved.') - parser.add_argument('vv_raster', - help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VV polarization') - parser.add_argument('water_extent_map', - help='HyP3-generated water extent raster file.') - parser.add_argument('hand_raster', - help='Height Above Nearest Drainage (HAND) GeoTIFF aligned to the RTC rasters. ' - 'If not specified, HAND data will be extracted from the GLO-30 HAND.') + parser.add_argument('out_raster', help='File to which flood depth map will be saved.') + parser.add_argument( + 'vv_raster', + help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VV polarization', + ) + parser.add_argument('water_extent_map', help='HyP3-generated water extent raster file.') + parser.add_argument( + 'hand_raster', + help='Height Above Nearest Drainage (HAND) GeoTIFF aligned to the RTC rasters. ' + 'If not specified, HAND data will be extracted from the GLO-30 HAND.', + ) else: raise NotImplementedError(f'Unknown interface: {interface}') - parser.add_argument('--estimator', type=optional_str, default='iterative', choices=available_estimators, - help=estimator_help) - parser.add_argument('--water-level-sigma', type=float, default=3., - 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 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.') + parser.add_argument( + '--estimator', + type=optional_str, + default='iterative', + choices=available_estimators, + help=estimator_help, + ) + parser.add_argument( + '--water-level-sigma', + type=float, + default=3.0, + 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 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, - 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.') + 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': - 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.') + 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}') @@ -342,7 +435,11 @@ def hyp3(): args = parser.parse_args() level = logging.DEBUG if args.verbose else logging.INFO - logging.basicConfig(stream=sys.stdout, format='%(asctime)s - %(levelname)s - %(message)s', level=level) + logging.basicConfig( + stream=sys.stdout, + format='%(asctime)s - %(levelname)s - %(message)s', + level=level, + ) log.debug(' '.join(sys.argv)) if args.estimator is None: @@ -370,10 +467,16 @@ def hyp3(): flood_map_raster = product_dir / f'{product_name}.tif' make_flood_map( - 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), - iterative_min_size=args.iterative_min_size, minimization_metric=args.minimization_metric, + 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), + iterative_min_size=args.iterative_min_size, + minimization_metric=args.minimization_metric, ) log.info(f'Flood depth map created successfully: {flood_map_raster}') @@ -390,14 +493,24 @@ def main(): args = parser.parse_args() level = logging.DEBUG if args.verbose else logging.INFO - logging.basicConfig(stream=sys.stdout, format='%(asctime)s - %(levelname)s - %(message)s', level=level) + logging.basicConfig( + stream=sys.stdout, + format='%(asctime)s - %(levelname)s - %(message)s', + level=level, + ) log.debug(' '.join(sys.argv)) 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, + 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}") + log.info(f'Flood depth map created successfully: {args.out_raster}') diff --git a/src/asf_tools/hydrosar/hand/__init__.py b/src/asf_tools/hydrosar/hand/__init__.py index 41dcff0f..558982d2 100644 --- a/src/asf_tools/hydrosar/hand/__init__.py +++ b/src/asf_tools/hydrosar/hand/__init__.py @@ -1,4 +1,7 @@ -from asf_tools.hydrosar.hand.calculate import calculate_hand_for_basins, make_copernicus_hand +from asf_tools.hydrosar.hand.calculate import ( + calculate_hand_for_basins, + make_copernicus_hand, +) from asf_tools.hydrosar.hand.prepare import prepare_hand_vrt __all__ = [ diff --git a/src/asf_tools/hydrosar/hand/calculate.py b/src/asf_tools/hydrosar/hand/calculate.py index 9fe0d7f1..015d95c0 100644 --- a/src/asf_tools/hydrosar/hand/calculate.py +++ b/src/asf_tools/hydrosar/hand/calculate.py @@ -1,4 +1,5 @@ """Calculate Height Above Nearest Drainage (HAND) from the Copernicus GLO-30 DEM""" + import argparse import logging import sys @@ -29,11 +30,9 @@ def fill_nan(array: np.ndarray) -> np.ndarray: """ kernel = astropy.convolution.Gaussian2DKernel(x_stddev=3) # kernel x_size=8*stddev with warnings.catch_warnings(): - warnings.simplefilter("ignore") + warnings.simplefilter('ignore') while np.any(np.isnan(array)): - array = astropy.convolution.interpolate_replace_nans( - array, kernel, convolve=astropy.convolution.convolve - ) + array = astropy.convolution.interpolate_replace_nans(array, kernel, convolve=astropy.convolution.convolve) return array @@ -55,8 +54,13 @@ def fill_hand(hand: np.ndarray, dem: np.ndarray): return hand -def calculate_hand(dem_array, dem_affine: rasterio.Affine, dem_crs: rasterio.crs.CRS, basin_mask, - acc_thresh: Optional[int] = 100): +def calculate_hand( + dem_array, + dem_affine: rasterio.Affine, + dem_crs: rasterio.crs.CRS, + basin_mask, + acc_thresh: Optional[int] = 100, +): """Calculate the Height Above Nearest Drainage (HAND) Calculate the Height Above Nearest Drainage (HAND) using pySHEDS library. Because HAND @@ -88,10 +92,14 @@ def calculate_hand(dem_array, dem_affine: rasterio.Affine, dem_crs: rasterio.crs """ nodata_fill_value = np.finfo(float).eps with NamedTemporaryFile() as temp_file: - write_cog(temp_file.name, dem_array, - transform=dem_affine.to_gdal(), epsg_code=dem_crs.to_epsg(), - # Prevents PySheds from assuming using zero as the nodata value - nodata_value=nodata_fill_value) + write_cog( + temp_file.name, + dem_array, + transform=dem_affine.to_gdal(), + epsg_code=dem_crs.to_epsg(), + # Prevents PySheds from assuming using zero as the nodata value + nodata_value=nodata_fill_value, + ) # From PySheds; see example usage: http://mattbartos.com/pysheds/ grid = sGrid.from_raster(str(temp_file.name)) @@ -134,8 +142,12 @@ def calculate_hand(dem_array, dem_affine: rasterio.Affine, dem_crs: rasterio.crs return hand -def calculate_hand_for_basins(out_raster: Union[str, Path], geometries: GeometryCollection, - dem_file: Union[str, Path], acc_thresh: Optional[int] = 100): +def calculate_hand_for_basins( + out_raster: Union[str, Path], + geometries: GeometryCollection, + dem_file: Union[str, Path], + acc_thresh: Optional[int] = 100, +): """Calculate the Height Above Nearest Drainage (HAND) for watershed boundaries (hydrobasins). For watershed boundaries, see: https://www.hydrosheds.org/page/hydrobasins @@ -156,11 +168,19 @@ def calculate_hand_for_basins(out_raster: Union[str, Path], geometries: Geometr hand = calculate_hand(basin_array, basin_affine_tf, src.crs, basin_mask, acc_thresh=acc_thresh) write_cog( - out_raster, hand, transform=basin_affine_tf.to_gdal(), epsg_code=src.crs.to_epsg(), nodata_value=np.nan, + out_raster, + hand, + transform=basin_affine_tf.to_gdal(), + epsg_code=src.crs.to_epsg(), + nodata_value=np.nan, ) -def make_copernicus_hand(out_raster: Union[str, Path], vector_file: Union[str, Path], acc_thresh: Optional[int] = 100): +def make_copernicus_hand( + out_raster: Union[str, Path], + vector_file: Union[str, Path], + acc_thresh: Optional[int] = 100, +): """Copernicus GLO-30 Height Above Nearest Drainage (HAND) Make a Height Above Nearest Drainage (HAND) GeoTIFF from the Copernicus GLO-30 DEM @@ -195,18 +215,30 @@ def main(): formatter_class=argparse.RawDescriptionHelpFormatter, ) parser.add_argument('out_raster', help='HAND GeoTIFF to create') - parser.add_argument('vector_file', help='Vector file of watershed boundary (hydrobasin) polygons to calculate HAND ' - 'over. Vector file Must be openable by GDAL, see: ' - 'https://gdal.org/drivers/vector/index.html') - parser.add_argument('-a', '--acc-threshold', type=none_or_int, default=100, - help='Accumulation threshold for determining the drainage mask. ' - 'If `None`, the mean accumulation value is used') + parser.add_argument( + 'vector_file', + help='Vector file of watershed boundary (hydrobasin) polygons to calculate HAND ' + 'over. Vector file Must be openable by GDAL, see: ' + 'https://gdal.org/drivers/vector/index.html', + ) + parser.add_argument( + '-a', + '--acc-threshold', + type=none_or_int, + default=100, + help='Accumulation threshold for determining the drainage mask. ' + 'If `None`, the mean accumulation value is used', + ) parser.add_argument('-v', '--verbose', action='store_true', help='Turn on verbose logging') args = parser.parse_args() level = logging.DEBUG if args.verbose else logging.INFO - logging.basicConfig(stream=sys.stdout, format='%(asctime)s - %(levelname)s - %(message)s', level=level) + logging.basicConfig( + stream=sys.stdout, + format='%(asctime)s - %(levelname)s - %(message)s', + level=level, + ) log.debug(' '.join(sys.argv)) log.info(f'Calculating HAND for {args.vector_file}') diff --git a/src/asf_tools/hydrosar/hand/prepare.py b/src/asf_tools/hydrosar/hand/prepare.py index fd7f755b..8b993598 100644 --- a/src/asf_tools/hydrosar/hand/prepare.py +++ b/src/asf_tools/hydrosar/hand/prepare.py @@ -1,4 +1,5 @@ """Prepare a Height Above Nearest Drainage (HAND) virtual raster (VRT) covering a given geometry""" + from pathlib import Path from tempfile import NamedTemporaryFile from typing import Union @@ -36,7 +37,7 @@ def prepare_hand_vrt(vrt: Union[str, Path], geometry: Union[ogr.Geometry, BaseGe geometry = ogr.CreateGeometryFromWkb(geometry.wkb) min_lon, max_lon, _, _ = geometry.GetEnvelope() - if min_lon < -160. and max_lon > 160.: + if min_lon < -160.0 and max_lon > 160.0: raise ValueError(f'asf_tools does not currently support geometries that cross the antimeridian: {geometry}') tile_features = vector.get_features(HAND_GEOJSON) @@ -48,8 +49,11 @@ def prepare_hand_vrt(vrt: Union[str, Path], geometry: Union[ogr.Geometry, BaseGe gdal.BuildVRT(str(vrt), hand_file_paths) -def prepare_hand_for_raster(hand_raster: Union[str, Path], source_raster: Union[str, Path], - resampling_method: str = 'lanczos'): +def prepare_hand_for_raster( + hand_raster: Union[str, Path], + source_raster: Union[str, Path], + resampling_method: str = 'lanczos', +): """Create a HAND raster pixel-aligned to a source raster Args: @@ -61,13 +65,21 @@ def prepare_hand_for_raster(hand_raster: Union[str, Path], source_raster: Union[ info = gdal.Info(str(source_raster), format='json') hand_geometry = shape(info['wgs84Extent']) - hand_bounds = [info['cornerCoordinates']['upperLeft'][0], - info['cornerCoordinates']['lowerRight'][1], - info['cornerCoordinates']['lowerRight'][0], - info['cornerCoordinates']['upperLeft'][1]] + hand_bounds = [ + info['cornerCoordinates']['upperLeft'][0], + info['cornerCoordinates']['lowerRight'][1], + info['cornerCoordinates']['lowerRight'][0], + info['cornerCoordinates']['upperLeft'][1], + ] with NamedTemporaryFile(suffix='.vrt', delete=False) as hand_vrt: prepare_hand_vrt(hand_vrt.name, hand_geometry) - gdal.Warp(str(hand_raster), hand_vrt.name, dstSRS=f'EPSG:{get_epsg_code(info)}', - outputBounds=hand_bounds, width=info['size'][0], height=info['size'][1], - resampleAlg=Resampling[resampling_method].value) + gdal.Warp( + str(hand_raster), + hand_vrt.name, + dstSRS=f'EPSG:{get_epsg_code(info)}', + outputBounds=hand_bounds, + width=info['size'][0], + height=info['size'][1], + resampleAlg=Resampling[resampling_method].value, + ) diff --git a/src/asf_tools/hydrosar/threshold.py b/src/asf_tools/hydrosar/threshold.py index a05662c7..8e5767d4 100644 --- a/src/asf_tools/hydrosar/threshold.py +++ b/src/asf_tools/hydrosar/threshold.py @@ -1,5 +1,3 @@ -# Turned off flake8 because we haven't refactored 3rd party provided functions -# flake8: noqa import numpy as np @@ -21,7 +19,7 @@ def _make_histogram(image): histogram[0, floor_value] = histogram[0, floor_value] + temp1 histogram[0, floor_value - 1] = histogram[0, floor_value - 1] + temp2 histogram = np.convolve(histogram[0], [1, 2, 3, 2, 1]) - histogram = histogram[2:(histogram.size - 3)] + histogram = histogram[2 : (histogram.size - 3)] histogram /= np.sum(histogram) return histogram @@ -68,43 +66,26 @@ def expectation_maximization_threshold(tile: np.ndarray, number_of_classes: int nonzero_indices = np.nonzero(histogram)[0] histogram = histogram[nonzero_indices] histogram = histogram.flatten() - class_means = ( - (np.arange(number_of_classes) + 1) * maximum / - (number_of_classes + 1) - ) + class_means = (np.arange(number_of_classes) + 1) * maximum / (number_of_classes + 1) class_variances = np.ones(number_of_classes) * maximum class_proportions = np.ones(number_of_classes) * 1 / number_of_classes sml = np.mean(np.diff(nonzero_indices)) / 1000 iteration = 0 while True: - class_likelihood = _make_distribution( - class_means, class_variances, class_proportions, nonzero_indices - ) - sum_likelihood = np.sum(class_likelihood, 1) + np.finfo( - class_likelihood[0][0]).eps + class_likelihood = _make_distribution(class_means, class_variances, class_proportions, nonzero_indices) + sum_likelihood = np.sum(class_likelihood, 1) + np.finfo(class_likelihood[0][0]).eps log_likelihood = np.sum(histogram * np.log(sum_likelihood)) for j in range(0, number_of_classes): - class_posterior_probability = ( - histogram * class_likelihood[:, j] / sum_likelihood - ) + class_posterior_probability = histogram * class_likelihood[:, j] / sum_likelihood class_proportions[j] = np.sum(class_posterior_probability) - class_means[j] = ( - np.sum(nonzero_indices * class_posterior_probability) - / class_proportions[j] - ) - vr = (nonzero_indices - class_means[j]) - class_variances[j] = ( - np.sum(vr * vr * class_posterior_probability) - / class_proportions[j] + sml - ) + class_means[j] = np.sum(nonzero_indices * class_posterior_probability) / class_proportions[j] + vr = nonzero_indices - class_means[j] + class_variances[j] = np.sum(vr * vr * class_posterior_probability) / class_proportions[j] + sml del class_posterior_probability, vr class_proportions += 1e-3 class_proportions /= np.sum(class_proportions) - class_likelihood = _make_distribution( - class_means, class_variances, class_proportions, nonzero_indices - ) - sum_likelihood = np.sum(class_likelihood, 1) + np.finfo( - class_likelihood[0, 0]).eps + class_likelihood = _make_distribution(class_means, class_variances, class_proportions, nonzero_indices) + sum_likelihood = np.sum(class_likelihood, 1) + np.finfo(class_likelihood[0, 0]).eps del class_likelihood new_log_likelihood = np.sum(histogram * np.log(sum_likelihood)) del sum_likelihood @@ -126,16 +107,28 @@ def expectation_maximization_threshold(tile: np.ndarray, number_of_classes: int posterior_lookup.update({pixel_val: [0] * number_of_classes}) for n in range(0, number_of_classes): x = _make_distribution( - class_means[n], class_variances[n], class_proportions[n], - image_copy[i, j] + class_means[n], + class_variances[n], + class_proportions[n], + image_copy[i, j], ) posterior[i, j, n] = x * class_proportions[n] posterior_lookup[pixel_val][n] = posterior[i, j, n] sorti = np.argsort(class_means) - xvec = np.arange(class_means[sorti[0]], class_means[sorti[1]], step=.05) - x1 = _make_distribution(class_means[sorti[0]], class_variances[sorti[0]], class_proportions[sorti[0]], xvec) - x2 = _make_distribution(class_means[sorti[1]], class_variances[sorti[1]], class_proportions[sorti[1]], xvec) + xvec = np.arange(class_means[sorti[0]], class_means[sorti[1]], step=0.05) + x1 = _make_distribution( + class_means[sorti[0]], + class_variances[sorti[0]], + class_proportions[sorti[0]], + xvec, + ) + x2 = _make_distribution( + class_means[sorti[1]], + class_variances[sorti[1]], + class_proportions[sorti[1]], + xvec, + ) dx = np.abs(x1 - x2) return xvec[np.argmin(dx)] diff --git a/src/asf_tools/hydrosar/water_map.py b/src/asf_tools/hydrosar/water_map.py index b293830e..ab464978 100644 --- a/src/asf_tools/hydrosar/water_map.py +++ b/src/asf_tools/hydrosar/water_map.py @@ -20,7 +20,9 @@ from asf_tools.aws import get_path_to_s3_file, upload_file_to_s3 from asf_tools.hydrosar.hand.prepare import prepare_hand_for_raster -from asf_tools.hydrosar.threshold import expectation_maximization_threshold as em_threshold +from asf_tools.hydrosar.threshold import ( + expectation_maximization_threshold as em_threshold, +) from asf_tools.raster import read_as_masked_array, write_cog from asf_tools.tile import tile_array, untile_array from asf_tools.util import get_epsg_code @@ -37,11 +39,16 @@ def mean_of_subtiles(tiles: np.ndarray) -> np.ndarray: return sub_tiles_mean -def select_hand_tiles(tiles: Union[np.ndarray, np.ma.MaskedArray], - hand_threshold: float, hand_fraction: float) -> np.ndarray: +def select_hand_tiles( + tiles: Union[np.ndarray, np.ma.MaskedArray], + hand_threshold: float, + hand_fraction: float, +) -> np.ndarray: if np.allclose(tiles, 0.0): - raise ValueError(f'All pixels in scene have a HAND value of {0.0} (all water); ' - f'scene is not a good candidate for water mapping.') + raise ValueError( + f'All pixels in scene have a HAND value of {0.0} (all water); ' + f'scene is not a good candidate for water mapping.' + ) tile_indexes = np.arange(tiles.shape[0]) @@ -84,14 +91,15 @@ def determine_em_threshold(tiles: np.ndarray, scaling: float) -> float: def calculate_slope_magnitude(array: np.ndarray, pixel_size) -> np.ndarray: dx, dy = np.gradient(array) - magnitude = np.sqrt(dx ** 2, dy ** 2) / pixel_size - slope = np.arctan(magnitude) / np.pi * 180. + magnitude = np.sqrt(dx**2, dy**2) / pixel_size + slope = np.arctan(magnitude) / np.pi * 180.0 return slope def determine_membership_limits( - array: np.ndarray, mask_percentile: float = 90., std_range: float = 3.0) -> Tuple[float, float]: - array = np.ma.masked_values(array, 0.) + array: np.ndarray, mask_percentile: float = 90.0, std_range: float = 3.0 +) -> Tuple[float, float]: + array = np.ma.masked_values(array, 0.0) array = np.ma.masked_greater(array, np.nanpercentile(array.filled(np.nan), mask_percentile)) lower_limit = np.ma.median(array) upper_limit = lower_limit + std_range * array.std() + 5.0 @@ -119,7 +127,11 @@ def segment_area_membership(segments: np.ndarray, min_area: int = 3, max_area: i for area in possible_areas: mask = np.isin(segments, (segment_areas == area).nonzero()) - np.putmask(segment_membership, mask, fuzz.interp_membership(possible_areas, activation, area)) + np.putmask( + segment_membership, + mask, + fuzz.interp_membership(possible_areas, activation, area), + ) return segment_membership @@ -146,35 +158,49 @@ def format_raster_data(raster, padding_mask=None, nodata=np.iinfo(np.uint8).max) return raster -def fuzzy_refinement(initial_map: np.ndarray, gaussian_array: np.ndarray, hand_array: np.ndarray, pixel_size: float, - gaussian_thresholds: Tuple[float, float], membership_threshold: float = 0.45) -> np.ndarray: +def fuzzy_refinement( + initial_map: np.ndarray, + gaussian_array: np.ndarray, + hand_array: np.ndarray, + pixel_size: float, + gaussian_thresholds: Tuple[float, float], + membership_threshold: float = 0.45, +) -> np.ndarray: water_map = np.ones_like(initial_map) water_segments = measure.label(initial_map, connectivity=2) water_segment_membership = segment_area_membership(water_segments) - water_map &= ~np.isclose(water_segment_membership, 0.) + water_map &= ~np.isclose(water_segment_membership, 0.0) gaussian_membership = min_max_membership(gaussian_array, gaussian_thresholds[0], gaussian_thresholds[1], 0.005) - water_map &= ~np.isclose(gaussian_membership, 0.) + water_map &= ~np.isclose(gaussian_membership, 0.0) hand_lower_limit, hand_upper_limit = determine_membership_limits(hand_array) hand_membership = min_max_membership(hand_array, hand_lower_limit, hand_upper_limit, 0.1) - water_map &= ~np.isclose(hand_membership, 0.) + water_map &= ~np.isclose(hand_membership, 0.0) hand_slopes = calculate_slope_magnitude(hand_array, pixel_size) - slope_membership = min_max_membership(hand_slopes, 0., 15., 0.1) - water_map &= ~np.isclose(slope_membership, 0.) + slope_membership = min_max_membership(hand_slopes, 0.0, 15.0, 0.1) + water_map &= ~np.isclose(slope_membership, 0.0) - water_map_weights = (gaussian_membership + hand_membership + slope_membership + water_segment_membership) / 4. + water_map_weights = (gaussian_membership + hand_membership + slope_membership + water_segment_membership) / 4.0 water_map &= water_map_weights >= membership_threshold return water_map -def make_water_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], vh_raster: Union[str, Path], - hand_raster: Optional[Union[str, Path]] = None, tile_shape: Tuple[int, int] = (100, 100), - max_vv_threshold: float = -15.5, max_vh_threshold: float = -23.0, - hand_threshold: float = 15., hand_fraction: float = 0.8, membership_threshold: float = 0.45): +def make_water_map( + out_raster: Union[str, Path], + vv_raster: Union[str, Path], + vh_raster: Union[str, Path], + hand_raster: Optional[Union[str, Path]] = None, + tile_shape: Tuple[int, int] = (100, 100), + max_vv_threshold: float = -15.5, + max_vh_threshold: float = -23.0, + hand_threshold: float = 15.0, + hand_fraction: float = 0.8, + membership_threshold: float = 0.45, +): """Creates a surface water extent map from a Sentinel-1 RTC product Create a surface water extent map from a dual-pol Sentinel-1 RTC product and @@ -253,25 +279,28 @@ def make_water_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], vh selected_tiles = None nodata = np.iinfo(np.uint8).max water_extent_maps = [] - for max_db_threshold, raster, pol in ((max_vh_threshold, vh_raster, 'VH'), (max_vv_threshold, vv_raster, 'VV')): + for max_db_threshold, raster, pol in ( + (max_vh_threshold, vh_raster, 'VH'), + (max_vv_threshold, vv_raster, 'VV'), + ): log.info(f'Creating initial {pol} water extent map from {raster}') array = read_as_masked_array(raster) padding_mask = array.mask - tiles = tile_array(array, tile_shape=tile_shape, pad_value=0.) + tiles = tile_array(array, tile_shape=tile_shape, pad_value=0.0) # Masking less than zero only necessary for old HyP3/GAMMA products which sometimes returned negative powers - tiles = np.ma.masked_less_equal(tiles, 0.) + tiles = np.ma.masked_less_equal(tiles, 0.0) if selected_tiles is None: selected_tiles = select_backscatter_tiles(tiles, hand_candidates) log.info(f'Selected tiles {selected_tiles} from {raster}') with np.testing.suppress_warnings() as sup: sup.filter(RuntimeWarning) # invalid value and divide by zero encountered in log10 - tiles = np.log10(tiles) + 30. # linear power scale -> Gaussian scale optimized for thresholding - max_gaussian_threshold = max_db_threshold / 10. + 30. # db -> Gaussian scale optimized for thresholding + tiles = np.log10(tiles) + 30.0 # linear power scale -> Gaussian scale optimized for thresholding + max_gaussian_threshold = max_db_threshold / 10.0 + 30.0 # db -> Gaussian scale optimized for thresholding if selected_tiles.size: scaling = 256 / (np.mean(tiles) + 3 * np.std(tiles)) gaussian_threshold = determine_em_threshold(tiles[selected_tiles, :, :], scaling) - threshold_db = 10. * (gaussian_threshold - 30.) + threshold_db = 10.0 * (gaussian_threshold - 30.0) log.info(f'Threshold determined to be {threshold_db} db') if gaussian_threshold > max_gaussian_threshold: log.warning(f'Threshold too high! Using maximum threshold {max_db_threshold} db') @@ -284,23 +313,37 @@ def make_water_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], vh water_map = np.ma.masked_less_equal(gaussian_array, gaussian_threshold).mask water_map &= ~array.mask - write_cog(str(out_raster).replace('.tif', f'_{pol}_initial.tif'), - format_raster_data(water_map, padding_mask, nodata), - transform=out_transform, epsg_code=out_epsg, dtype=gdal.GDT_Byte, nodata_value=nodata) + write_cog( + str(out_raster).replace('.tif', f'_{pol}_initial.tif'), + format_raster_data(water_map, padding_mask, nodata), + transform=out_transform, + epsg_code=out_epsg, + dtype=gdal.GDT_Byte, + nodata_value=nodata, + ) log.info(f'Refining initial {pol} water extent map using Fuzzy Logic') array = np.ma.masked_where(~water_map, array) - gaussian_lower_limit = np.log10(np.ma.median(array)) + 30. + gaussian_lower_limit = np.log10(np.ma.median(array)) + 30.0 water_map = fuzzy_refinement( - water_map, gaussian_array, hand_array, pixel_size=out_transform[1], - gaussian_thresholds=(gaussian_lower_limit, gaussian_threshold), membership_threshold=membership_threshold + water_map, + gaussian_array, + hand_array, + pixel_size=out_transform[1], + gaussian_thresholds=(gaussian_lower_limit, gaussian_threshold), + membership_threshold=membership_threshold, ) water_map &= ~array.mask - write_cog(str(out_raster).replace('.tif', f'_{pol}_fuzzy.tif'), - format_raster_data(water_map, padding_mask, nodata), - transform=out_transform, epsg_code=out_epsg, dtype=gdal.GDT_Byte, nodata_value=nodata) + write_cog( + str(out_raster).replace('.tif', f'_{pol}_fuzzy.tif'), + format_raster_data(water_map, padding_mask, nodata), + transform=out_transform, + epsg_code=out_epsg, + dtype=gdal.GDT_Byte, + nodata_value=nodata, + ) water_extent_maps.append(water_map) @@ -310,47 +353,83 @@ def make_water_map(out_raster: Union[str, Path], vv_raster: Union[str, Path], vh combined_segments = measure.label(combined_water_map, connectivity=2) combined_water_map = remove_small_segments(combined_segments) - write_cog(out_raster, format_raster_data(combined_water_map, padding_mask, nodata), transform=out_transform, - epsg_code=out_epsg, dtype=gdal.GDT_Byte, nodata_value=nodata) + write_cog( + out_raster, + format_raster_data(combined_water_map, padding_mask, nodata), + transform=out_transform, + epsg_code=out_epsg, + dtype=gdal.GDT_Byte, + nodata_value=nodata, + ) def _get_cli(interface: Literal['hyp3', 'main']) -> argparse.ArgumentParser: - parser = argparse.ArgumentParser( - description=__doc__, - formatter_class=argparse.ArgumentDefaultsHelpFormatter - ) + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) if interface == 'hyp3': parser.add_argument('--bucket') parser.add_argument('--bucket-prefix', default='') - parser.add_argument('--vv-raster', - help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VV polarization.') + parser.add_argument( + '--vv-raster', + help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VV polarization.', + ) elif interface == 'main': parser.add_argument('out_raster', help='Water map GeoTIFF to create') # FIXME: Decibel RTCs would be real nice. - parser.add_argument('vv_raster', - help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VV polarization') - parser.add_argument('vh_raster', - help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VH polarization') - - parser.add_argument('--hand-raster', - help='Height Above Nearest Drainage (HAND) GeoTIFF aligned to the RTC rasters. ' - 'If not specified, HAND data will be extracted from the GLO-30 HAND.') - parser.add_argument('--tile-shape', type=int, nargs=2, default=(100, 100), - help='image tiles will have this shape (height, width) in pixels') + parser.add_argument( + 'vv_raster', + help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VV polarization', + ) + parser.add_argument( + 'vh_raster', + help='Sentinel-1 RTC GeoTIFF raster, in power scale, with VH polarization', + ) + + parser.add_argument( + '--hand-raster', + help='Height Above Nearest Drainage (HAND) GeoTIFF aligned to the RTC rasters. ' + 'If not specified, HAND data will be extracted from the GLO-30 HAND.', + ) + parser.add_argument( + '--tile-shape', + type=int, + nargs=2, + default=(100, 100), + help='image tiles will have this shape (height, width) in pixels', + ) else: raise NotImplementedError(f'Unknown interface: {interface}') - parser.add_argument('--max-vv-threshold', type=float, default=-15.5, - help='Maximum threshold value to use for `vv_raster` in decibels (db)') - parser.add_argument('--max-vh-threshold', type=float, default=-23.0, - help='Maximum threshold value to use for `vh_raster` in decibels (db)') - parser.add_argument('--hand-threshold', type=float, default=15., - help='The maximum height above nearest drainage in meters to consider a pixel valid') - parser.add_argument('--hand-fraction', type=float, default=0.8, - help='The minimum fraction of valid HAND pixels required in a tile for thresholding') - parser.add_argument('--membership-threshold', type=float, default=0.45, - help='The average membership to the fuzzy indicators required for a water pixel') + parser.add_argument( + '--max-vv-threshold', + type=float, + default=-15.5, + help='Maximum threshold value to use for `vv_raster` in decibels (db)', + ) + parser.add_argument( + '--max-vh-threshold', + type=float, + default=-23.0, + help='Maximum threshold value to use for `vh_raster` in decibels (db)', + ) + parser.add_argument( + '--hand-threshold', + type=float, + default=15.0, + help='The maximum height above nearest drainage in meters to consider a pixel valid', + ) + parser.add_argument( + '--hand-fraction', + type=float, + default=0.8, + help='The minimum fraction of valid HAND pixels required in a tile for thresholding', + ) + parser.add_argument( + '--membership-threshold', + type=float, + default=0.45, + help='The average membership to the fuzzy indicators required for a water pixel', + ) parser.add_argument('-v', '--verbose', action='store_true', help='Turn on verbose logging') @@ -362,7 +441,11 @@ def hyp3(): args = parser.parse_args() level = logging.DEBUG if args.verbose else logging.INFO - logging.basicConfig(stream=sys.stdout, format='%(asctime)s - %(levelname)s - %(message)s', level=level) + logging.basicConfig( + stream=sys.stdout, + format='%(asctime)s - %(levelname)s - %(message)s', + level=level, + ) log.debug(' '.join(sys.argv)) if args.vv_raster: @@ -382,10 +465,14 @@ def hyp3(): water_map_raster = product_dir / f'{product_name}.tif' make_water_map( - out_raster=water_map_raster, vv_raster=vv_raster, vh_raster=vh_raster, - max_vv_threshold=args.max_vv_threshold, max_vh_threshold=args.max_vh_threshold, - hand_threshold=args.hand_threshold, hand_fraction=args.hand_fraction, - membership_threshold=args.membership_threshold + out_raster=water_map_raster, + vv_raster=vv_raster, + vh_raster=vh_raster, + max_vv_threshold=args.max_vv_threshold, + max_vh_threshold=args.max_vh_threshold, + hand_threshold=args.hand_threshold, + hand_fraction=args.hand_fraction, + membership_threshold=args.membership_threshold, ) log.info(f'Water map created successfully: {water_map_raster}') @@ -402,11 +489,24 @@ def main(): args = parser.parse_args() level = logging.DEBUG if args.verbose else logging.INFO - logging.basicConfig(stream=sys.stdout, format='%(asctime)s - %(levelname)s - %(message)s', level=level) + logging.basicConfig( + stream=sys.stdout, + format='%(asctime)s - %(levelname)s - %(message)s', + level=level, + ) log.debug(' '.join(sys.argv)) - make_water_map(args.out_raster, args.vv_raster, args.vh_raster, args.hand_raster, args.tile_shape, - args.max_vv_threshold, args.max_vh_threshold, args.hand_threshold, args.hand_fraction, - args.membership_threshold) + make_water_map( + args.out_raster, + args.vv_raster, + args.vh_raster, + args.hand_raster, + args.tile_shape, + args.max_vv_threshold, + args.max_vh_threshold, + args.hand_threshold, + args.hand_fraction, + args.membership_threshold, + ) log.info(f'Water map created successfully: {args.out_raster}') diff --git a/src/asf_tools/raster.py b/src/asf_tools/raster.py index 096ae1cf..f030857f 100644 --- a/src/asf_tools/raster.py +++ b/src/asf_tools/raster.py @@ -13,8 +13,11 @@ log = logging.getLogger(__name__) -def convert_scale(array: Union[np.ndarray, np.ma.MaskedArray], in_scale: Literal['db', 'amplitude', 'power'], - out_scale: Literal['db', 'amplitude', 'power']) -> Union[np.ndarray, np.ma.MaskedArray]: +def convert_scale( + array: Union[np.ndarray, np.ma.MaskedArray], + in_scale: Literal['db', 'amplitude', 'power'], + out_scale: Literal['db', 'amplitude', 'power'], +) -> Union[np.ndarray, np.ma.MaskedArray]: """Convert calibrated raster scale between db, amplitude and power""" if in_scale == out_scale: warnings.warn(f'Nothing to do! {in_scale} is same as {out_scale}.') @@ -30,9 +33,9 @@ def convert_scale(array: Union[np.ndarray, np.ma.MaskedArray], in_scale: Literal if in_scale == 'amplitude': if out_scale == 'power': - return array ** 2 + return array**2 if out_scale == 'db': - return 10 * log10(array ** 2) + return 10 * log10(array**2) if in_scale == 'power': if out_scale == 'amplitude': @@ -81,8 +84,14 @@ def read_as_array(raster: str, band: int = 1) -> np.array: return data -def write_cog(file_name: Union[str, Path], data: np.ndarray, transform: List[float], epsg_code: int, - dtype=gdal.GDT_Float32, nodata_value=None): +def write_cog( + file_name: Union[str, Path], + data: np.ndarray, + transform: List[float], + epsg_code: int, + dtype=gdal.GDT_Float32, + nodata_value=None, +): """Creates a Cloud Optimized GeoTIFF Args: @@ -108,7 +117,12 @@ def write_cog(file_name: Union[str, Path], data: np.ndarray, transform: List[flo temp_geotiff.SetProjection(epsg_to_wkt(epsg_code)) driver = gdal.GetDriverByName('COG') - options = ['COMPRESS=LZW', 'OVERVIEW_RESAMPLING=AVERAGE', 'NUM_THREADS=ALL_CPUS', 'BIGTIFF=YES'] + options = [ + 'COMPRESS=LZW', + 'OVERVIEW_RESAMPLING=AVERAGE', + 'NUM_THREADS=ALL_CPUS', + 'BIGTIFF=YES', + ] driver.CreateCopy(str(file_name), temp_geotiff, options=options) del temp_geotiff # How to close w/ gdal diff --git a/src/asf_tools/tile.py b/src/asf_tools/tile.py index 1cf2f310..c3e12fba 100644 --- a/src/asf_tools/tile.py +++ b/src/asf_tools/tile.py @@ -3,8 +3,11 @@ import numpy as np -def tile_array(array: Union[np.ndarray, np.ma.MaskedArray], tile_shape: Tuple[int, int] = (200, 200), - pad_value: float = None) -> Union[np.ndarray, np.ma.MaskedArray]: +def tile_array( + array: Union[np.ndarray, np.ma.MaskedArray], + tile_shape: Tuple[int, int] = (200, 200), + pad_value: float = None, +) -> Union[np.ndarray, np.ma.MaskedArray]: """Tile a 2D numpy array Turn a 2D numpy array like: @@ -65,8 +68,9 @@ def tile_array(array: Union[np.ndarray, np.ma.MaskedArray], tile_shape: Tuple[in return tiled -def untile_array(tiled_array: Union[np.ndarray, np.ma.MaskedArray], array_shape: Tuple[int, int]) \ - -> Union[np.ndarray, np.ma.MaskedArray]: +def untile_array( + tiled_array: Union[np.ndarray, np.ma.MaskedArray], array_shape: Tuple[int, int] +) -> Union[np.ndarray, np.ma.MaskedArray]: """Untile a tiled array into a 2D numpy array This is the reverse of `tile_array` and will turn a tiled array like: @@ -106,7 +110,10 @@ def untile_array(tiled_array: Union[np.ndarray, np.ma.MaskedArray], array_shape: untiled_rows = int(np.ceil(array_rows / tile_rows)) untiled_columns = int(np.ceil(array_columns / tile_columns)) - untiled = np.zeros((untiled_rows*tile_rows, untiled_columns*tile_columns), dtype=tiled_array.dtype) + untiled = np.zeros( + (untiled_rows * tile_rows, untiled_columns * tile_columns), + dtype=tiled_array.dtype, + ) if (array_size := array_rows * array_columns) > tiled_array.size: raise ValueError( @@ -116,8 +123,10 @@ def untile_array(tiled_array: Union[np.ndarray, np.ma.MaskedArray], array_shape: for ii in range(untiled_rows): for jj in range(untiled_columns): - untiled[ii*tile_rows:(ii+1)*tile_rows, jj*tile_columns:(jj+1)*tile_columns] = \ - tiled_array[ii * untiled_columns + jj, :, :] + untiled[ + ii * tile_rows : (ii + 1) * tile_rows, + jj * tile_columns : (jj + 1) * tile_columns, + ] = tiled_array[ii * untiled_columns + jj, :, :] if isinstance(tiled_array, np.ma.MaskedArray): untiled_mask = untile_array(tiled_array.mask, untiled.shape) diff --git a/src/asf_tools/util.py b/src/asf_tools/util.py index ed25bbbb..8af4f485 100644 --- a/src/asf_tools/util.py +++ b/src/asf_tools/util.py @@ -7,6 +7,7 @@ class GDALConfigManager: """Context manager for setting GDAL config options temporarily""" + def __init__(self, **options): """ Args: diff --git a/src/asf_tools/vector.py b/src/asf_tools/vector.py index a6707499..3f278f94 100644 --- a/src/asf_tools/vector.py +++ b/src/asf_tools/vector.py @@ -18,8 +18,7 @@ def get_property_values_for_intersecting_features(geometry: ogr.Geometry, featur return True -def intersecting_feature_properties(geometry: ogr.Geometry, features: Iterator, - feature_property: str) -> List[str]: +def intersecting_feature_properties(geometry: ogr.Geometry, features: Iterator, feature_property: str) -> List[str]: property_values = [] for feature in features: if feature.GetGeometryRef().Intersects(geometry): diff --git a/src/asf_tools/watermasking/fill_missing_tiles.py b/src/asf_tools/watermasking/fill_missing_tiles.py index ea9b64a2..bddaa6d9 100644 --- a/src/asf_tools/watermasking/fill_missing_tiles.py +++ b/src/asf_tools/watermasking/fill_missing_tiles.py @@ -12,17 +12,32 @@ def main(): - parser = argparse.ArgumentParser( prog='fill_missing_tiles.py', - description='Script for creating filled tifs in areas with missing tiles.' + description='Script for creating filled tifs in areas with missing tiles.', ) parser.add_argument('--fill-value', help='The value to fill the data array with.', default=0) - parser.add_argument('--lat-begin', help='The minimum latitude of the dataset in EPSG:4326.', default=-85) - parser.add_argument('--lat-end', help='The maximum latitude of the dataset in EPSG:4326.', default=85) - parser.add_argument('--lon-begin', help='The minimum longitude of the dataset in EPSG:4326.', default=-180) - parser.add_argument('--lon-end', help='The maximum longitude of the dataset in EPSG:4326.', default=180) + parser.add_argument( + '--lat-begin', + help='The minimum latitude of the dataset in EPSG:4326.', + default=-85, + ) + parser.add_argument( + '--lat-end', + help='The maximum latitude of the dataset in EPSG:4326.', + default=85, + ) + parser.add_argument( + '--lon-begin', + help='The minimum longitude of the dataset in EPSG:4326.', + default=-180, + ) + parser.add_argument( + '--lon-end', + help='The maximum longitude of the dataset in EPSG:4326.', + default=180, + ) parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) @@ -41,7 +56,6 @@ def main(): for lat in lat_range: for lon in lon_range: - tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') tile_tif = 'tiles/' + tile + '.tif' tile_cog = 'tiles/cogs/' + tile + '.tif' @@ -57,7 +71,13 @@ def main(): data.fill(fill_value) driver = gdal.GetDriverByName('GTiff') - dst_ds = driver.Create(tile_tif, xsize=data.shape[0], ysize=data.shape[1], bands=1, eType=gdal.GDT_Byte) + dst_ds = driver.Create( + tile_tif, + xsize=data.shape[0], + ysize=data.shape[1], + bands=1, + eType=gdal.GDT_Byte, + ) dst_ds.SetGeoTransform([xmin, pixel_size_x, 0, ymin, 0, pixel_size_y]) srs = osr.SpatialReference() srs.ImportFromEPSG(4326) diff --git a/src/asf_tools/watermasking/generate_osm_tiles.py b/src/asf_tools/watermasking/generate_osm_tiles.py index 5233c64b..d42b6816 100644 --- a/src/asf_tools/watermasking/generate_osm_tiles.py +++ b/src/asf_tools/watermasking/generate_osm_tiles.py @@ -6,7 +6,11 @@ import geopandas as gpd from osgeo import gdal -from asf_tools.watermasking.utils import lat_lon_to_tile_string, remove_temp_files, setup_directories +from asf_tools.watermasking.utils import ( + lat_lon_to_tile_string, + remove_temp_files, + setup_directories, +) gdal.UseExceptions() @@ -54,7 +58,7 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig tile = lat_lon_to_tile_string(lat, lon, is_worldcover=False, postfix='') tile_tif = output_dir + tile + '.tif' - xmin, xmax, ymin, ymax = lon, lon+tile_width_deg, lat, lat+tile_height_deg + xmin, xmax, ymin, ymax = lon, lon + tile_width_deg, lat, lat + tile_height_deg pixel_size_x = 0.00009009009 pixel_size_y = 0.00009009009 @@ -70,10 +74,16 @@ def process_ocean_tiles(ocean_polygons_path, lat, lon, tile_width_deg, tile_heig burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, - creationOptions=GDAL_OPTIONS + creationOptions=GDAL_OPTIONS, ) - temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx', tile + '.shp'] + temp_files = [ + tile + '.dbf', + tile + '.cpg', + tile + '.prj', + tile + '.shx', + tile + '.shp', + ] remove_temp_files(temp_files) @@ -112,7 +122,7 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio water_gdf.to_file(tile_shp, mode='w', engine='pyogrio') water_gdf = None - xmin, xmax, ymin, ymax = lon, lon+tile_width_deg, lat, lat+tile_height_deg + xmin, xmax, ymin, ymax = lon, lon + tile_width_deg, lat, lat + tile_height_deg pixel_size_x = 0.00009009009 pixel_size_y = 0.00009009009 @@ -124,10 +134,18 @@ def extract_water(water_file, lat, lon, tile_width_deg, tile_height_deg, interio burnValues=1, outputBounds=[xmin, ymin, xmax, ymax], outputType=gdal.GDT_Byte, - creationOptions=GDAL_OPTIONS + creationOptions=GDAL_OPTIONS, ) - temp_files = [tile + '.dbf', tile + '.cpg', tile + '.prj', tile + '.shx', tile_shp, tile_pbf, tile_geojson] + temp_files = [ + tile + '.dbf', + tile + '.cpg', + tile + '.prj', + tile + '.shx', + tile_shp, + tile_pbf, + tile_geojson, + ] remove_temp_files(temp_files) @@ -159,7 +177,7 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di '--outfile', output_tile, '--calc', - '"logical_or(A, B)"' + '"logical_or(A, B)"', ] subprocess.run(command) @@ -185,18 +203,33 @@ def merge_interior_and_ocean(internal_tile_dir, ocean_tile_dir, finished_tile_di def main(): - parser = argparse.ArgumentParser( prog='generate_osm_tiles.py', - description='Main script for creating a tiled watermask dataset from OSM data.' + description='Main script for creating a tiled watermask dataset from OSM data.', ) parser.add_argument('--planet-file-path', help='The path to the global planet.pbf file.') parser.add_argument('--ocean-polygons-path', help='The path to the global OSM ocean polygons.') - parser.add_argument('--lat-begin', help='The minimum latitude of the dataset in EPSG:4326.', default=-85) - parser.add_argument('--lat-end', help='The maximum latitude of the dataset in EPSG:4326.', default=85) - parser.add_argument('--lon-begin', help='The minimum longitude of the dataset in EPSG:4326.', default=-180) - parser.add_argument('--lon-end', help='The maximum longitude of the dataset in EPSG:4326.', default=180) + parser.add_argument( + '--lat-begin', + help='The minimum latitude of the dataset in EPSG:4326.', + default=-85, + ) + parser.add_argument( + '--lat-end', + help='The maximum latitude of the dataset in EPSG:4326.', + default=85, + ) + parser.add_argument( + '--lon-begin', + help='The minimum longitude of the dataset in EPSG:4326.', + default=-180, + ) + parser.add_argument( + '--lon-end', + help='The maximum longitude of the dataset in EPSG:4326.', + default=180, + ) parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) @@ -235,7 +268,7 @@ def main(): lon, tile_width, tile_height, - interior_tile_dir=INTERIOR_TILE_DIR + interior_tile_dir=INTERIOR_TILE_DIR, ) process_ocean_tiles( args.ocean_polygons_path, @@ -243,7 +276,7 @@ def main(): lon, tile_width, tile_height, - output_dir=OCEAN_TILE_DIR + output_dir=OCEAN_TILE_DIR, ) end_time = time.time() total_time = end_time - start_time diff --git a/src/asf_tools/watermasking/generate_worldcover_tiles.py b/src/asf_tools/watermasking/generate_worldcover_tiles.py index 508f8883..92fd90d5 100644 --- a/src/asf_tools/watermasking/generate_worldcover_tiles.py +++ b/src/asf_tools/watermasking/generate_worldcover_tiles.py @@ -6,7 +6,12 @@ import numpy as np from osgeo import gdal -from asf_tools.watermasking.utils import lat_lon_to_tile_string, merge_tiles, remove_temp_files, setup_directories +from asf_tools.watermasking.utils import ( + lat_lon_to_tile_string, + merge_tiles, + remove_temp_files, + setup_directories, +) gdal.UseExceptions() @@ -41,12 +46,12 @@ def filename_filter(filename): in_lat_range = (latitude >= mnlat) and (latitude <= mxlat) in_lon_range = (longitude >= mnlon) and (longitude <= mxlon) return in_lat_range and in_lon_range + filenames_filtered = [f for f in filenames if filename_filter(f)] index = 0 num_tiles = len(filenames_filtered) for filename in filenames_filtered: - start_time = time.time() tile_name = filename.split('_')[5] @@ -70,7 +75,7 @@ def filename_filter(filename): water_arr.shape[1], 1, gdal.GDT_Byte, - options=GDAL_OPTIONS + options=GDAL_OPTIONS, ) dst_ds.SetGeoTransform(src_ds.GetGeoTransform()) dst_ds.SetProjection(src_ds.GetProjection()) @@ -120,7 +125,7 @@ def create_missing_tiles(tile_dir, lat_range, lon_range): ysize=y_size, bands=1, eType=gdal.GDT_Byte, - options=GDAL_OPTIONS + options=GDAL_OPTIONS, ) ds.SetProjection('EPSG:4326') ds.SetGeoTransform(geotransform) @@ -198,18 +203,18 @@ def crop_tile(tile, lat, lon, tile_width, tile_height): gdal.Translate( out_filename, src_ds, - projWin=[lon, lat+tile_height, lon+tile_width, lat], + projWin=[lon, lat + tile_height, lon + tile_width, lat], xRes=pixel_size_x, yRes=pixel_size_y, outputSRS='EPSG:4326', format='COG', - creationOptions=['NUM_THREADS=all_cpus'] + creationOptions=['NUM_THREADS=all_cpus'], ) remove_temp_files(['tmp_px_size.tif', 'tmp.shp']) def build_dataset(worldcover_tile_dir, lat_range, lon_range, tile_width, tile_height): - """ Main function for generating a dataset with worldcover tiles. + """Main function for generating a dataset with worldcover tiles. Args: worldcover_tile_dir: The directory containing the unprocessed worldcover tiles. @@ -232,22 +237,36 @@ def build_dataset(worldcover_tile_dir, lat_range, lon_range, tile_width, tile_he def main(): - parser = argparse.ArgumentParser( prog='generate_worldcover_tiles.py', - description='Main script for creating a tiled watermask dataset from the ESA WorldCover dataset.' + description='Main script for creating a tiled watermask dataset from the ESA WorldCover dataset.', ) - parser.add_argument('--worldcover-tiles-dir', help='The path to the directory containing the worldcover tifs.') + parser.add_argument( + '--worldcover-tiles-dir', + help='The path to the directory containing the worldcover tifs.', + ) parser.add_argument( '--lat-begin', help='The minimum latitude of the dataset in EPSG:4326.', default=-85, - required=True + required=True, + ) + parser.add_argument( + '--lat-end', + help='The maximum latitude of the dataset in EPSG:4326.', + default=85, + ) + parser.add_argument( + '--lon-begin', + help='The minimum longitude of the dataset in EPSG:4326.', + default=-180, + ) + parser.add_argument( + '--lon-end', + help='The maximum longitude of the dataset in EPSG:4326.', + default=180, ) - parser.add_argument('--lat-end', help='The maximum latitude of the dataset in EPSG:4326.', default=85) - parser.add_argument('--lon-begin', help='The minimum longitude of the dataset in EPSG:4326.', default=-180) - parser.add_argument('--lon-end', help='The maximum longitude of the dataset in EPSG:4326.', default=180) parser.add_argument('--tile-width', help='The desired width of the tile in degrees.', default=5) parser.add_argument('--tile-height', help='The desired height of the tile in degrees.', default=5) @@ -270,12 +289,12 @@ def main(): wc_lat_range = range( lat_begin - (lat_begin % WORLDCOVER_TILE_SIZE), lat_end + (lat_end % WORLDCOVER_TILE_SIZE), - WORLDCOVER_TILE_SIZE + WORLDCOVER_TILE_SIZE, ) wc_lon_range = range( lon_begin - (lon_begin % WORLDCOVER_TILE_SIZE), lon_end + (lon_end % WORLDCOVER_TILE_SIZE), - WORLDCOVER_TILE_SIZE + WORLDCOVER_TILE_SIZE, ) # Ocean only tiles are missing from WorldCover, so we need to create blank (water-only) ones. @@ -286,7 +305,7 @@ def main(): lat_range, lon_range, tile_width=tile_width, - tile_height=tile_height + tile_height=tile_height, ) diff --git a/src/asf_tools/watermasking/utils.py b/src/asf_tools/watermasking/utils.py index 7806be0d..7f7449b0 100644 --- a/src/asf_tools/watermasking/utils.py +++ b/src/asf_tools/watermasking/utils.py @@ -43,11 +43,14 @@ def merge_tiles(tiles, out_filename, out_format, compress=False): else: translate_command = [ 'gdal_translate', - '-of', out_format, - '-co', 'COMPRESS=LZW', - '-co', 'NUM_THREADS=all_cpus', + '-of', + out_format, + '-co', + 'COMPRESS=LZW', + '-co', + 'NUM_THREADS=all_cpus', vrt, - out_filename + out_filename, ] subprocess.run(build_vrt_command) subprocess.run(translate_command) diff --git a/tests/conftest.py b/tests/conftest.py index 96364eac..7fa46136 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,14 +2,17 @@ def pytest_configure(config): - config.addinivalue_line( - "markers", "integration: marks tests as integration" - ) + config.addinivalue_line('markers', 'integration: marks tests as integration') def pytest_addoption(parser): - parser.addoption('--integration', action='store_true', default=False, dest='integration', - help='enable integration tests') + parser.addoption( + '--integration', + action='store_true', + default=False, + dest='integration', + help='enable integration tests', + ) def pytest_collection_modifyitems(config, items): diff --git a/tests/hydrosar/test_flood_map.py b/tests/hydrosar/test_flood_map.py index 072dcc2e..f92062e8 100644 --- a/tests/hydrosar/test_flood_map.py +++ b/tests/hydrosar/test_flood_map.py @@ -9,14 +9,18 @@ @pytest.mark.integration def test_get_waterbody(): - water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif' + water_raster = ( + '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif' + ) info = gdal.Info(water_raster, format='json') known_water_mask = flood_map.get_waterbody(info, threshold=30) - test_mask = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' \ - 'flood_map/known_water_mask.tif' + test_mask = ( + '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' + 'flood_map/known_water_mask.tif' + ) test_mask_array = gdal.Open(test_mask, gdal.GA_ReadOnly).ReadAsArray() assert np.all(known_water_mask == test_mask_array) @@ -42,35 +46,58 @@ def test_estimate_flood_depths_iterative(flood_window, hand_window): @pytest.mark.integration def test_estimate_flood_depths_logstat(flood_window, hand_window): - water_height = flood_map.estimate_flood_depth(1, hand_window, flood_window, estimator='logstat', - water_level_sigma=3, - iterative_bounds=(0, 15)) + water_height = flood_map.estimate_flood_depth( + 1, + hand_window, + flood_window, + estimator='logstat', + water_level_sigma=3, + iterative_bounds=(0, 15), + ) assert water_height == 21.02364492416382 @pytest.mark.integration def test_estimate_flood_depths_nmad(flood_window, hand_window): - water_height = flood_map.estimate_flood_depth(1, hand_window, flood_window, estimator='nmad', water_level_sigma=3, - iterative_bounds=(0, 15)) + water_height = flood_map.estimate_flood_depth( + 1, + hand_window, + flood_window, + estimator='nmad', + water_level_sigma=3, + iterative_bounds=(0, 15), + ) assert water_height == 7.887911175434299 @pytest.mark.integration def test_estimate_flood_depths_numpy(flood_window, hand_window): - water_height = flood_map.estimate_flood_depth(1, hand_window, flood_window, estimator='numpy', water_level_sigma=3, - iterative_bounds=(0, 15)) + water_height = flood_map.estimate_flood_depth( + 1, + hand_window, + flood_window, + estimator='numpy', + water_level_sigma=3, + iterative_bounds=(0, 15), + ) assert water_height == 16.353520154953003 @pytest.mark.integration def test_make_flood_map(tmp_path): - water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif' - vv_raster = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/RTC_VV.tif' - hand_raster = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap_HAND.tif' + water_raster = ( + '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif' + ) + vv_raster = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/RTC_VV.tif' + ) + hand_raster = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap_HAND.tif' + ) out_flood_map = tmp_path / 'flood_map.tif' flood_map.make_flood_map(out_flood_map, vv_raster, water_raster, hand_raster, estimator='nmad') @@ -78,9 +105,11 @@ def test_make_flood_map(tmp_path): assert out_flood_map.exists() - golden_flood_map = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' \ - 'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' \ - 'flood_map/flood_map_nmad.tif' + golden_flood_map = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' + 'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/' + 'flood_map/flood_map_nmad.tif' + ) diffs = find_diff(golden_flood_map, str(out_flood_map)) assert diffs == 0 diff --git a/tests/hydrosar/test_hand.py b/tests/hydrosar/test_hand.py index 035b3714..59291268 100644 --- a/tests/hydrosar/test_hand.py +++ b/tests/hydrosar/test_hand.py @@ -8,10 +8,14 @@ from asf_tools import vector from asf_tools.hydrosar import hand -HAND_BASINS = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' \ - 'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/hand/hybas_af_lev12_v1c_firstpoly.geojson' -GOLDEN_HAND = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' \ - 'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/hand/hybas_af_lev12_v1c_firstpoly.tif' +HAND_BASINS = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' + 'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/hand/hybas_af_lev12_v1c_firstpoly.geojson' +) +GOLDEN_HAND = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/' + 'asf-tools/S1A_IW_20230228T120437_DVR_RTC30/hand/hybas_af_lev12_v1c_firstpoly.tif' +) gdal.UseExceptions() @@ -29,7 +33,6 @@ def nodata_equal_nan(golden_hand, out_hand): @pytest.mark.integration def test_make_copernicus_hand(tmp_path): - out_hand = tmp_path / 'hand.tif' hand.make_copernicus_hand(out_hand, HAND_BASINS) @@ -53,13 +56,15 @@ def test_prepare_hand_vrt(tmp_path): hand_vrt = tmp_path / 'hand.tif' geojson = { 'type': 'Polygon', - 'coordinates': [[ - [0.4, 10.16], - [0.4, 10.86], - [0.6, 10.86], - [0.6, 10.16], - [0.4, 10.16], - ]], + 'coordinates': [ + [ + [0.4, 10.16], + [0.4, 10.86], + [0.6, 10.86], + [0.6, 10.16], + [0.4, 10.16], + ] + ], } geometry = ogr.CreateGeometryFromJson(json.dumps(geojson)) @@ -67,8 +72,14 @@ def test_prepare_hand_vrt(tmp_path): assert hand_vrt.exists() info = gdal.Info(str(hand_vrt), format='json') - assert info['geoTransform'] == \ - [-0.0001388888888889, 0.0002777777777778, 0.0, 11.00013888888889, 0.0, -0.0002777777777778] + assert info['geoTransform'] == [ + -0.0001388888888889, + 0.0002777777777778, + 0.0, + 11.00013888888889, + 0.0, + -0.0002777777777778, + ] assert info['size'] == [3600, 3600] @@ -76,20 +87,24 @@ def test_prepare_hand_vrt_antimeridian(): geojson = { 'type': 'MultiPolygon', 'coordinates': [ - [[ - [179.5, 51.4], - [179.5, 51.6], - [180.0, 51.6], - [180.0, 51.4], - [179.5, 51.4], - ]], - [[ - [-180.0, 51.4], - [-180.0, 51.6], - [-179.5, 51.6], - [-179.5, 51.4], - [-180.0, 51.4], - ]], + [ + [ + [179.5, 51.4], + [179.5, 51.6], + [180.0, 51.6], + [180.0, 51.4], + [179.5, 51.4], + ] + ], + [ + [ + [-180.0, 51.4], + [-180.0, 51.6], + [-179.5, 51.6], + [-179.5, 51.4], + [-180.0, 51.4], + ] + ], ], } geometry = ogr.CreateGeometryFromJson(json.dumps(geojson)) diff --git a/tests/hydrosar/test_water_map.py b/tests/hydrosar/test_water_map.py index d844e7cf..443786b8 100644 --- a/tests/hydrosar/test_water_map.py +++ b/tests/hydrosar/test_water_map.py @@ -19,18 +19,20 @@ def test_select_hand_tiles(hand_candidates): hand_array = read_as_array(str(hand_geotif)) hand_tiles = np.ma.masked_invalid(tile_array(hand_array, tile_shape=(100, 100), pad_value=np.nan)) - selected_tiles = water_map.select_hand_tiles(hand_tiles, 15., 0.8) + selected_tiles = water_map.select_hand_tiles(hand_tiles, 15.0, 0.8) assert np.all(selected_tiles == hand_candidates) with pytest.raises(ValueError): - _ = water_map.select_hand_tiles(np.zeros(shape=(10, 10, 10), dtype=float), 15., 0.8) + _ = water_map.select_hand_tiles(np.zeros(shape=(10, 10, 10), dtype=float), 15.0, 0.8) @pytest.mark.integration def test_select_backscatter_tiles(hand_candidates): backscatter_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/water-map/20200603_VH.tif' backscatter_array = np.ma.masked_invalid(read_as_array(backscatter_geotif)) - backscatter_tiles = np.ma.masked_less_equal(tile_array(backscatter_array, tile_shape=(100, 100), pad_value=0.), 0.) + backscatter_tiles = np.ma.masked_less_equal( + tile_array(backscatter_array, tile_shape=(100, 100), pad_value=0.0), 0.0 + ) tile_indexes = water_map.select_backscatter_tiles(backscatter_tiles, hand_candidates) assert np.all(tile_indexes == np.array([771, 1974, 2397, 1205, 2577])) @@ -38,19 +40,27 @@ def test_select_backscatter_tiles(hand_candidates): @pytest.mark.integration def test_make_water_map(tmp_path): - vv_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/water_map/RTC_VV.tif' - vh_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/water_map/RTC_VH.tif' - hand_geotif = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/water_map/HAND.tif' + vv_geotif = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/water_map/RTC_VV.tif' + ) + vh_geotif = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/water_map/RTC_VH.tif' + ) + hand_geotif = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/water_map/HAND.tif' + ) out_water_map = tmp_path / 'water_map.tif' water_map.make_water_map(out_water_map, vv_geotif, vh_geotif, hand_geotif) assert out_water_map.exists() - golden_water_map = '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/water_map/fuzzy_water_map.tif' + golden_water_map = ( + '/vsicurl/https://hyp3-testing.s3-us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/water_map/fuzzy_water_map.tif' + ) diffs = find_diff(golden_water_map, str(out_water_map)) assert diffs == 0 diff --git a/tests/test_aws.py b/tests/test_aws.py index 5a71ed3a..ffe767d4 100644 --- a/tests/test_aws.py +++ b/tests/test_aws.py @@ -12,14 +12,7 @@ def s3_stubber(): def test_get_tag_set(): - assert aws.get_tag_set() == { - 'TagSet': [ - { - 'Key': 'file_type', - 'Value': 'product' - } - ] - } + assert aws.get_tag_set() == {'TagSet': [{'Key': 'file_type', 'Value': 'product'}]} def test_get_content_type(): @@ -40,11 +33,7 @@ def test_upload_file_to_s3(tmp_path, s3_stubber): tag_params = { 'Bucket': 'myBucket', 'Key': 'myFile.zip', - 'Tagging': { - 'TagSet': [ - {'Key': 'file_type', 'Value': 'product'} - ] - } + 'Tagging': {'TagSet': [{'Key': 'file_type', 'Value': 'product'}]}, } s3_stubber.add_response(method='put_object', expected_params=expected_params, service_response={}) s3_stubber.add_response(method='put_object_tagging', expected_params=tag_params, service_response={}) @@ -64,11 +53,7 @@ def test_upload_file_to_s3_with_prefix(tmp_path, s3_stubber): tag_params = { 'Bucket': 'myBucket', 'Key': 'myPrefix/myFile.txt', - 'Tagging': { - 'TagSet': [ - {'Key': 'file_type', 'Value': 'product'} - ] - } + 'Tagging': {'TagSet': [{'Key': 'file_type', 'Value': 'product'}]}, } s3_stubber.add_response(method='put_object', expected_params=expected_params, service_response={}) s3_stubber.add_response(method='put_object_tagging', expected_params=tag_params, service_response={}) @@ -91,14 +76,23 @@ def test_get_path_to_s3_file(s3_stubber): ], } - s3_stubber.add_response(method='list_objects_v2', expected_params=expected_params, - service_response=service_response) + s3_stubber.add_response( + method='list_objects_v2', + expected_params=expected_params, + service_response=service_response, + ) assert aws.get_path_to_s3_file('myBucket', 'myPrefix', file_type='.nc') == '/vsis3/myBucket/myPrefix/foo.nc' - s3_stubber.add_response(method='list_objects_v2', expected_params=expected_params, - service_response=service_response) + s3_stubber.add_response( + method='list_objects_v2', + expected_params=expected_params, + service_response=service_response, + ) assert aws.get_path_to_s3_file('myBucket', 'myPrefix', file_type='.txt') == '/vsis3/myBucket/myPrefix/foo.txt' - s3_stubber.add_response(method='list_objects_v2', expected_params=expected_params, - service_response=service_response) + s3_stubber.add_response( + method='list_objects_v2', + expected_params=expected_params, + service_response=service_response, + ) assert aws.get_path_to_s3_file('myBucket', 'myPrefix', file_type='.csv') is None diff --git a/tests/test_composite.py b/tests/test_composite.py index 9c58a83b..d1e4ee0f 100644 --- a/tests/test_composite.py +++ b/tests/test_composite.py @@ -25,9 +25,7 @@ def test_get_target_epsg_code(): assert composite.get_target_epsg_code([32701, 32760, 32701]) == 32701 assert composite.get_target_epsg_code([32701, 32760, 32760]) == 32760 - assert composite.get_target_epsg_code( - [32731, 32631, 32731, 32631, 32732, 32633, 32733, 32633, 32733] - ) == 32732 + assert composite.get_target_epsg_code([32731, 32631, 32731, 32631, 32732, 32633, 32733, 32633, 32733]) == 32732 # bounds with pytest.raises(ValueError): @@ -67,7 +65,11 @@ def test_get_full_extents(): expected_upper_left = (10.0, 130.0) expected_lower_right = (110.0, 30.0) expected_geotransform = [10.0, 2.0, 0.0, 130.0, 0.0, -2.0] - assert composite.get_full_extent(data) == (expected_upper_left, expected_lower_right, expected_geotransform) + assert composite.get_full_extent(data) == ( + expected_upper_left, + expected_lower_right, + expected_geotransform, + ) data['b'] = { 'cornerCoordinates': { @@ -80,7 +82,11 @@ def test_get_full_extents(): expected_upper_left = (10.0, 140.0) expected_lower_right = (120.0, 30.0) expected_geotransform = [10.0, 2.0, 0.0, 140.0, 0.0, -2.0] - assert composite.get_full_extent(data) == (expected_upper_left, expected_lower_right, expected_geotransform) + assert composite.get_full_extent(data) == ( + expected_upper_left, + expected_lower_right, + expected_geotransform, + ) def test_make_composite(tmp_path): @@ -88,26 +94,34 @@ def test_make_composite(tmp_path): epsg_code = 32601 transform = [0.0, 30.0, 0.0, 60.0, 0.0, -30.0] - data = np.array([ - [1, 1, 1, 1], - [1, 1, 1, 1], - ]) - area = np.array([ - [1, 1, 1, 1], - [1, 1, 1, 1], - ]) + data = np.array( + [ + [1, 1, 1, 1], + [1, 1, 1, 1], + ] + ) + area = np.array( + [ + [1, 1, 1, 1], + [1, 1, 1, 1], + ] + ) asf_tools.raster.write_cog('first_data.tif', data, transform, epsg_code, nodata_value=0) asf_tools.raster.write_cog('first_area.tif', area, transform, epsg_code) transform = [30.0, 30.0, 0.0, 30.0, 0.0, -30.0] - data = np.array([ - [3, 0, 3, 3], - [3, 0, 3, 3], - ]) - area = np.array([ - [1, 1, 3, 1], - [1, 1, 2, 1], - ]) + data = np.array( + [ + [3, 0, 3, 3], + [3, 0, 3, 3], + ] + ) + area = np.array( + [ + [1, 1, 3, 1], + [1, 1, 2, 1], + ] + ) asf_tools.raster.write_cog('second_data.tif', data, transform, epsg_code) asf_tools.raster.write_cog('second_area.tif', area, transform, epsg_code) @@ -119,17 +133,21 @@ def test_make_composite(tmp_path): assert os.path.exists(count_file) data = np.nan_to_num(asf_tools.raster.read_as_array(out_file)) - expected = np.array([ - [1, 1, 1, 1, 0], - [1, 2, 1, 1.5, 3], - [0, 3, 0, 3, 3], - ]) + expected = np.array( + [ + [1, 1, 1, 1, 0], + [1, 2, 1, 1.5, 3], + [0, 3, 0, 3, 3], + ] + ) assert np.allclose(data, expected) counts = asf_tools.raster.read_as_array(count_file) - expected = np.array([ - [1, 1, 1, 1, 0], - [1, 2, 1, 2, 1], - [0, 1, 0, 1, 1], - ]) + expected = np.array( + [ + [1, 1, 1, 1, 0], + [1, 2, 1, 2, 1], + [0, 1, 0, 1, 1], + ] + ) assert np.allclose(counts, expected) diff --git a/tests/test_dem.py b/tests/test_dem.py index fe742783..8d1f0796 100644 --- a/tests/test_dem.py +++ b/tests/test_dem.py @@ -22,13 +22,15 @@ def test_prepare_dem_vrt(tmp_path): dem_vrt = tmp_path / 'dem.tif' geojson = { 'type': 'Polygon', - 'coordinates': [[ - [0.4, 10.16], - [0.4, 10.86], - [0.6, 10.86], - [0.6, 10.16], - [0.4, 10.16], - ]], + 'coordinates': [ + [ + [0.4, 10.16], + [0.4, 10.86], + [0.6, 10.86], + [0.6, 10.16], + [0.4, 10.16], + ] + ], } geometry = ogr.CreateGeometryFromJson(json.dumps(geojson)) @@ -36,8 +38,14 @@ def test_prepare_dem_vrt(tmp_path): assert dem_vrt.exists() info = gdal.Info(str(dem_vrt), format='json') - assert info['geoTransform'] == \ - [-0.0001388888888889, 0.0002777777777778, 0.0, 11.00013888888889, 0.0, -0.0002777777777778] + assert info['geoTransform'] == [ + -0.0001388888888889, + 0.0002777777777778, + 0.0, + 11.00013888888889, + 0.0, + -0.0002777777777778, + ] assert info['size'] == [3600, 3600] @@ -46,20 +54,24 @@ def test_prepare_dem_geotiff_antimeridian(tmp_path): geojson = { 'type': 'MultiPolygon', 'coordinates': [ - [[ - [179.5, 51.4], - [179.5, 51.6], - [180.0, 51.6], - [180.0, 51.4], - [179.5, 51.4], - ]], - [[ - [-180.0, 51.4], - [-180.0, 51.6], - [-179.5, 51.6], - [-179.5, 51.4], - [-180.0, 51.4], - ]], + [ + [ + [179.5, 51.4], + [179.5, 51.6], + [180.0, 51.6], + [180.0, 51.4], + [179.5, 51.4], + ] + ], + [ + [ + [-180.0, 51.4], + [-180.0, 51.6], + [-179.5, 51.6], + [-179.5, 51.4], + [-180.0, 51.4], + ] + ], ], } geometry = ogr.CreateGeometryFromJson(json.dumps(geojson)) diff --git a/tests/test_raster.py b/tests/test_raster.py index 6be64526..143bc123 100644 --- a/tests/test_raster.py +++ b/tests/test_raster.py @@ -11,7 +11,7 @@ def test_convert_scale(): assert np.allclose(c, np.array([100, 25, 0, 25, 100])) c = raster.convert_scale(np.array([-10, -5, 0, 5, 10]), 'amplitude', 'db') - assert np.allclose(c, np.array([20., 13.97940009, -np.inf, 13.97940009, 20.])) + assert np.allclose(c, np.array([20.0, 13.97940009, -np.inf, 13.97940009, 20.0])) c = raster.convert_scale(np.array([-1, 0, 1, 4, 9]), 'power', 'amplitude') assert np.isnan(c[0]) @@ -19,14 +19,17 @@ def test_convert_scale(): c = raster.convert_scale(np.array([-1, 0, 1, 4, 9]), 'power', 'db') assert np.isnan(c[0]) - assert np.allclose(c[1:], np.array([-np.inf, 0., 6.02059991, 9.54242509]), ) + assert np.allclose( + c[1:], + np.array([-np.inf, 0.0, 6.02059991, 9.54242509]), + ) - c = raster.convert_scale(np.array([np.nan, -np.inf, 0., 6.02059991, 9.54242509]), 'db', 'power') + c = raster.convert_scale(np.array([np.nan, -np.inf, 0.0, 6.02059991, 9.54242509]), 'db', 'power') assert np.isnan(c[0]) assert np.allclose(c[1:], np.array([0, 1, 4, 9])) - c = raster.convert_scale(np.array([-np.inf, -20., 0., 13.97940009, 20.]), 'db', 'amplitude') - assert np.allclose(c, np.array([0., 0.1, 1., 5., 10.])) + c = raster.convert_scale(np.array([-np.inf, -20.0, 0.0, 13.97940009, 20.0]), 'db', 'amplitude') + assert np.allclose(c, np.array([0.0, 0.1, 1.0, 5.0, 10.0])) a = np.array([-10, -5, 0, 5, 10]) with pytest.raises(ValueError): @@ -35,7 +38,10 @@ def test_convert_scale(): _ = raster.convert_scale(a, 'bar', 'amplitude') with pytest.warns(UserWarning): - assert np.allclose(raster.convert_scale(a, 'amplitude', 'amplitude'), np.array([-10, -5, 0, 5, 10])) + assert np.allclose( + raster.convert_scale(a, 'amplitude', 'amplitude'), + np.array([-10, -5, 0, 5, 10]), + ) with pytest.warns(UserWarning): assert np.allclose(raster.convert_scale(a, 'power', 'power'), np.array([-10, -5, 0, 5, 10])) with pytest.warns(UserWarning): @@ -47,7 +53,11 @@ def test_convert_scale_masked_arrays(): c = raster.convert_scale(masked_array, 'power', 'db') assert np.allclose(c.mask, [True, True, False, False, False]) assert np.allclose( - c, np.ma.MaskedArray([np.nan, -np.inf, 0., 6.02059991, 9.54242509], mask=[True, True, False, False, False]) + c, + np.ma.MaskedArray( + [np.nan, -np.inf, 0.0, 6.02059991, 9.54242509], + mask=[True, True, False, False, False], + ), ) a = raster.convert_scale(c, 'db', 'power') diff --git a/tests/test_tile.py b/tests/test_tile.py index ea87ff66..c888bda6 100644 --- a/tests/test_tile.py +++ b/tests/test_tile.py @@ -5,10 +5,7 @@ def test_tile_array(): - a = np.array([[0, 0, 1, 1], - [0, 0, 1, 1], - [2, 2, 3, 3], - [2, 2, 3, 3]]) + a = np.array([[0, 0, 1, 1], [0, 0, 1, 1], [2, 2, 3, 3], [2, 2, 3, 3]]) tiled = tile.tile_array(a, tile_shape=(2, 2)) assert tiled.shape == (4, 2, 2) @@ -35,18 +32,19 @@ def test_tile_array(): def test_tile_masked_array(): - a = np.array([[0, 0, 1, 1], - [0, 0, 1, 1], - [2, 2, 3, 3], - [2, 2, 3, 3]]) + a = np.array([[0, 0, 1, 1], [0, 0, 1, 1], [2, 2, 3, 3], [2, 2, 3, 3]]) with pytest.raises(AttributeError): _ = tile.tile_array(a, tile_shape=(2, 2)).mask - m = np.array([[False, False, False, True], - [False, False, False, False], - [False, False, False, False], - [False, False, False, True]]) + m = np.array( + [ + [False, False, False, True], + [False, False, False, False], + [False, False, False, False], + [False, False, False, True], + ] + ) ma = np.ma.MaskedArray(a, mask=m) tiled = tile.tile_array(ma, tile_shape=(2, 2)) @@ -54,14 +52,15 @@ def test_tile_masked_array(): assert tiled.shape == (4, 2, 2) assert isinstance(tiled, np.ma.MaskedArray) assert np.all( - tiled.mask == np.array([[[False, False], - [False, False]], - [[False, True], - [False, False]], - [[False, False], - [False, False]], - [[False, False], - [False, True]]]) + tiled.mask + == np.array( + [ + [[False, False], [False, False]], + [[False, True], [False, False]], + [[False, False], [False, False]], + [[False, False], [False, True]], + ] + ) ) tiled = tile.tile_array(ma, tile_shape=(3, 3), pad_value=4) @@ -72,19 +71,20 @@ def test_tile_masked_array(): tiled[0, :, :].mask == np.array([[False, False, False], [False, False, False], [False, False, False]]) ) assert np.all(np.ma.getdata(tiled[-1, :, :]) == np.array([[3, 4, 4], [4, 4, 4], [4, 4, 4]])) - assert np.all( - tiled[-1, :, :].mask == np.array([[True, True, True], [True, True, True], [True, True, True]]) - ) + assert np.all(tiled[-1, :, :].mask == np.array([[True, True, True], [True, True, True], [True, True, True]])) def test_untile_array(): - a = np.array([[0, 0, 1, 1, 2, 2], - [0, 0, 1, 1, 2, 2], - [3, 3, 4, 4, 5, 5], - [3, 3, 4, 4, 5, 5], - [6, 6, 7, 7, 8, 8], - [6, 6, 7, 7, 8, 8], - ]) + a = np.array( + [ + [0, 0, 1, 1, 2, 2], + [0, 0, 1, 1, 2, 2], + [3, 3, 4, 4, 5, 5], + [3, 3, 4, 4, 5, 5], + [6, 6, 7, 7, 8, 8], + [6, 6, 7, 7, 8, 8], + ] + ) assert np.all(a == tile.untile_array(tile.tile_array(a, tile_shape=(2, 2)), array_shape=a.shape)) assert np.all(a == tile.untile_array(tile.tile_array(a, tile_shape=(4, 4), pad_value=9), array_shape=a.shape)) @@ -99,24 +99,25 @@ def test_untile_array(): # array shape will subset some of the padding that was required to tile `a` with `tile_shape` assert np.all( - np.pad(a, ((0, 0), (0, 2)), constant_values=9) - == tile.untile_array(tile.tile_array(a, tile_shape=(2, 4), pad_value=9), array_shape=(6, 8)) + np.pad(a, ((0, 0), (0, 2)), constant_values=9) + == tile.untile_array(tile.tile_array(a, tile_shape=(2, 4), pad_value=9), array_shape=(6, 8)) ) def test_untile_masked_array(): - a = np.array([[0, 0, 1, 1], - [0, 0, 1, 1], - [2, 2, 3, 3], - [2, 2, 3, 3]]) + a = np.array([[0, 0, 1, 1], [0, 0, 1, 1], [2, 2, 3, 3], [2, 2, 3, 3]]) with pytest.raises(AttributeError): _ = tile.untile_array(tile.tile_array(a, tile_shape=(2, 2)), array_shape=a.shape).mask - m = np.array([[False, False, False, True], - [False, False, False, False], - [False, False, False, False], - [False, False, False, True]]) + m = np.array( + [ + [False, False, False, True], + [False, False, False, False], + [False, False, False, False], + [False, False, False, True], + ] + ) ma = np.ma.MaskedArray(a, mask=m) untiled = tile.untile_array(tile.tile_array(ma.copy(), tile_shape=(2, 2)), array_shape=a.shape) diff --git a/tests/test_util.py b/tests/test_util.py index 4366b6c5..59f3c169 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -22,8 +22,10 @@ def test_get_epsg_code(): def test_get_coordinates(): - water_raster = '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' \ - 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif' + water_raster = ( + '/vsicurl/https://hyp3-testing.s3.us-west-2.amazonaws.com/asf-tools/' + 'S1A_IW_20230228T120437_DVR_RTC30/flood_map/watermap.tif' + ) info = gdal.Info(water_raster, format='json') west, south, east, north = util.get_coordinates(info) diff --git a/tests/test_vector.py b/tests/test_vector.py index d8987b5f..24ddbf99 100644 --- a/tests/test_vector.py +++ b/tests/test_vector.py @@ -45,7 +45,7 @@ def test_get_intersecting_feature_properties(): geojson = { 'type': 'MultiPoint', - 'coordinates': [[0, 0], [169, -45], [-121.5, 73.5]] + 'coordinates': [[0, 0], [169, -45], [-121.5, 73.5]], } geometry = ogr.CreateGeometryFromJson(json.dumps(geojson)) assert vector.intersecting_feature_properties(geometry, dem_tile_features, 'file_path') == [