From afc2e185d1f638e8d859758f665f200cb01d2418 Mon Sep 17 00:00:00 2001 From: vikineema Date: Mon, 20 Nov 2023 19:07:59 +0300 Subject: [PATCH] Remove hardcoded buffer distance --- .../plugins/ocean_filtering_using_goas.py | 46 +++---------------- .../ocean_filtering_using_hydrosheds.py | 38 ++++----------- deafrica_waterbodies/plugins/utils.py | 29 ++++++++++++ deafrica_waterbodies/tiling.py | 11 ++--- 4 files changed, 48 insertions(+), 76 deletions(-) diff --git a/deafrica_waterbodies/plugins/ocean_filtering_using_goas.py b/deafrica_waterbodies/plugins/ocean_filtering_using_goas.py index 320d52e..4611cfc 100644 --- a/deafrica_waterbodies/plugins/ocean_filtering_using_goas.py +++ b/deafrica_waterbodies/plugins/ocean_filtering_using_goas.py @@ -5,60 +5,26 @@ import geopandas as gpd import numpy as np -import skimage.morphology import xarray as xr from deafrica_tools.spatial import xr_rasterize +from deafrica_waterbodies.plugins.utils import erode_land_sea_mask + # File extensions to recognise as Parquet files. PARQUET_EXTENSIONS = {".pq", ".parquet"} -buffer_pixels = 500 / 30 - - -def transform_hydrosheds_land_mask(hydrosheds_land_mask: xr.DataArray) -> xr.DataArray: - """ - Function to transform the HydroSHEDs Land Mask into a boolean mask where - 0/False are ocean pixels and 1/True are land pixels. - """ - # Indicator values: 1 = land, 2 = ocean sink, 3 = inland sink, 255 is no data. - boolean_mask = (hydrosheds_land_mask != 255) & (hydrosheds_land_mask != 2) - - return boolean_mask - - -def erode_land_sea_mask(boolean_land_sea_mask: xr.DataArray, buffer_pixels: float) -> xr.DataArray: - """ - Shrink the land in the land/sea mask. - - Parameters - ---------- - boolean_land_sea_mask : xr.DataArray - Boolean mask where 0/False are ocean pixels and 1/True are land pixels. - buffer_pixels : float - Number of pixels to erode the land by. - - Returns - ------- - xr.DataArray - Eroded land sea mask where 0/False are ocean pixels and 1/True are land pixels. - """ - eroded_boolean_land_sea_mask = xr.apply_ufunc( - skimage.morphology.binary_erosion, - boolean_land_sea_mask, - skimage.morphology.disk(buffer_pixels), - ) - return eroded_boolean_land_sea_mask - def load_land_sea_mask( land_sea_mask_fp: str, wofs_alltime_summary_ds: xr.DataArray, + buffer_dist_m: float = 500, ) -> xr.DataArray: """ Load the Marine Regions Global Oceans and Seas v01 from the file path provided. Rasterize the vector data to match the loaded datacube WOfS All Time Summary data and transform the raster into a boolean mask where 0/False are ocean pixels and 1/True are land pixels. + Erode the land pixels by the `buffer_dist_m` buffer distance. Parameters ---------- @@ -66,6 +32,8 @@ def load_land_sea_mask( File path to the Marine Regions Global Oceans and Seas v01 vector data. wofs_alltime_summary_ds : xr.DataArray Loaded datacube WOfS All Time Summary data to match to. + buffer_dist_m : float + Distance in meters to erode the land by in the land/sea mask. Returns ------- @@ -87,6 +55,6 @@ def load_land_sea_mask( boolean_land_sea_mask = np.logical_not(land_sea_mask_ds) # Erode the land in the land sea mask - eroded_boolean_land_sea_mask = erode_land_sea_mask(boolean_land_sea_mask, buffer_pixels) + eroded_boolean_land_sea_mask = erode_land_sea_mask(boolean_land_sea_mask, buffer_dist_m) return eroded_boolean_land_sea_mask diff --git a/deafrica_waterbodies/plugins/ocean_filtering_using_hydrosheds.py b/deafrica_waterbodies/plugins/ocean_filtering_using_hydrosheds.py index b5469e0..439df95 100644 --- a/deafrica_waterbodies/plugins/ocean_filtering_using_hydrosheds.py +++ b/deafrica_waterbodies/plugins/ocean_filtering_using_hydrosheds.py @@ -1,11 +1,10 @@ """ Ocean filtering using HydroSHEDS Land Mask """ -import skimage.morphology import xarray as xr from datacube.testutils.io import rio_slurp_xarray -buffer_pixels = 500 / 30 +from deafrica_waterbodies.plugins.utils import erode_land_sea_mask def transform_hydrosheds_land_mask(hydrosheds_land_mask: xr.DataArray) -> xr.DataArray: @@ -19,38 +18,16 @@ def transform_hydrosheds_land_mask(hydrosheds_land_mask: xr.DataArray) -> xr.Dat return boolean_mask -def erode_land_sea_mask(boolean_land_sea_mask: xr.DataArray, buffer_pixels: float) -> xr.DataArray: - """ - Shrink the land in the land/sea mask. - - Parameters - ---------- - boolean_land_sea_mask : xr.DataArray - Boolean mask where 0/False are ocean pixels and 1/True are land pixels. - buffer_pixels : float - Number of pixels to erode the land by. - - Returns - ------- - xr.DataArray - Eroded land sea mask where 0/False are ocean pixels and 1/True are land pixels. - """ - eroded_boolean_land_sea_mask = xr.apply_ufunc( - skimage.morphology.binary_erosion, - boolean_land_sea_mask, - skimage.morphology.disk(buffer_pixels), - ) - return eroded_boolean_land_sea_mask - - def load_land_sea_mask( land_sea_mask_fp: str, wofs_alltime_summary_ds: xr.DataArray, + buffer_dist_m: float = 500, ) -> xr.DataArray: """ Load and reproject the HydroSHEDS Land Mask raster from the file path provided to - match the loaded datacube WOfS All Time Summary data and transform the raster into - a boolean mask where 0/False are ocean pixels and 1/True are land pixels. + match the loaded datacube WOfS All Time Summary data. Transform the loaded raster into + a boolean mask where 0/False are ocean pixels and 1/True are land pixels and erode the land + pixels by the `buffer_dist_m` buffer distance. Parameters ---------- @@ -58,7 +35,8 @@ def load_land_sea_mask( File path to the HydroSHEDS Land Mask raster. wofs_alltime_summary_ds : xr.DataArray Loaded datacube WOfS All Time Summary data to match to - + buffer_dist_m : float + Distance in meters to erode the land by in the land/sea mask. Returns ------- xr.DataArray @@ -74,6 +52,6 @@ def load_land_sea_mask( boolean_land_sea_mask = transform_hydrosheds_land_mask(land_sea_mask_ds) # Erode the land in the land sea mask - eroded_boolean_land_sea_mask = erode_land_sea_mask(boolean_land_sea_mask, buffer_pixels) + eroded_boolean_land_sea_mask = erode_land_sea_mask(boolean_land_sea_mask, buffer_dist_m) return eroded_boolean_land_sea_mask diff --git a/deafrica_waterbodies/plugins/utils.py b/deafrica_waterbodies/plugins/utils.py index af50cac..1c4e8f6 100644 --- a/deafrica_waterbodies/plugins/utils.py +++ b/deafrica_waterbodies/plugins/utils.py @@ -7,6 +7,9 @@ from pathlib import Path from types import ModuleType +import skimage +import xarray as xr + def run_plugin(plugin_path: str | Path) -> ModuleType: """Run a Python plugin from a path. @@ -42,3 +45,29 @@ def validate_plugin(plugin: ModuleType): required_functions = ["load_land_sea_mask"] for name in required_functions: assert hasattr(getattr(plugin, name), "__call__") + + +def erode_land_sea_mask(boolean_land_sea_mask: xr.DataArray, buffer_dist_m: float) -> xr.DataArray: + """ + Shrink the land in the land/sea mask. + + Parameters + ---------- + boolean_land_sea_mask : xr.DataArray + Boolean mask where 0/False are ocean pixels and 1/True are land pixels. + buffer_dist_m : float + Distance in meters to erode the land by in the land/sea mask. + + Returns + ------- + xr.DataArray + Eroded land sea mask where 0/False are ocean pixels and 1/True are land pixels. + """ + buffer_pixels = buffer_dist_m / abs(boolean_land_sea_mask.geobox.resolution[0]) + + eroded_boolean_land_sea_mask = xr.apply_ufunc( + skimage.morphology.binary_erosion, + boolean_land_sea_mask, + skimage.morphology.disk(buffer_pixels), + ) + return eroded_boolean_land_sea_mask diff --git a/deafrica_waterbodies/tiling.py b/deafrica_waterbodies/tiling.py index 0ca412c..152da9c 100644 --- a/deafrica_waterbodies/tiling.py +++ b/deafrica_waterbodies/tiling.py @@ -29,22 +29,20 @@ def check_tile_intersects_polygons( tuple Tile id if the extent of the geobox of a tile intersects with the polygons. """ - tile_id = tile[0] - tile_extent = tile[1].geobox.extent if polygons_gdf is not None: # Reproject the extent of the geobox of a tile to match the polygons. - tile_extent = tile_extent.to_crs(polygons_gdf.crs) + tile_extent = tile[1].geobox.extent.to_crs(polygons_gdf.crs) # Get the shapely geometry of the reprojected extent of the tile's geobox. tile_extent_geom = tile_extent.geom # Check if the extent intersects with any of the polygons. check_intersection = polygons_gdf.geometry.intersects(tile_extent_geom).any() if check_intersection: - return tile_id + return tile[0] else: return () else: - return tile_id + return tile[0] def filter_tiles( @@ -105,8 +103,7 @@ def tile_wofs_ls_summary_alltime( # wofs_ls_summary_alltime grid. # Regular grid. - grid = "africa_30" - grid, gridspec = parse_gridspec_with_name(grid) + grid, gridspec = parse_gridspec_with_name(s="africa_30") # Multiply the tile size. tile_size = tuple(tile_size_factor * elem for elem in gridspec.tile_size)