From 96184e8bbaecc6b3fb3fde5e74f52b215c7f7baa Mon Sep 17 00:00:00 2001 From: rosepearson Date: Thu, 26 Oct 2023 16:37:23 +1300 Subject: [PATCH] Added support for only offshore data --- src/geofabrics/dem.py | 43 ++++++++++++++++++-------- src/geofabrics/geometry.py | 60 +++++++++++++++++++++++-------------- src/geofabrics/processor.py | 4 +++ 3 files changed, 72 insertions(+), 35 deletions(-) diff --git a/src/geofabrics/dem.py b/src/geofabrics/dem.py index 051e1628..3c81fc6c 100644 --- a/src/geofabrics/dem.py +++ b/src/geofabrics/dem.py @@ -527,6 +527,18 @@ def _sample_offshore_edge(self, resolution) -> numpy.ndarray: offshore_dense_data_edge = self.catchment_geometry.offshore_dense_data_edge( self._raw_extents ) + if offshore_dense_data_edge.area.sum() == 0: + # No offshore edge. Return an empty array. + offshore_edge = numpy.empty( + [0], + dtype=[ + ("X", geometry.RASTER_TYPE), + ("Y", geometry.RASTER_TYPE), + ("Z", geometry.RASTER_TYPE), + ], + ) + return offshore_edge + # Otherwise proceed as normal offshore_edge_dem = self._raw_dem.rio.clip(offshore_dense_data_edge.geometry) # If the sampling resolution is coaser than the catchment_geometry resolution @@ -990,11 +1002,14 @@ def _define_chunk_region( ) # Define region to rasterise inside the chunk area - chunk_region_to_tile = geopandas.GeoDataFrame( - geometry=region_to_rasterise.buffer(radius).clip( - chunk_geometry.buffer(radius), keep_geom_type=True + if region_to_rasterise.area.sum() > 0: + chunk_region_to_tile = geopandas.GeoDataFrame( + geometry=region_to_rasterise.buffer(radius).clip( + chunk_geometry.buffer(radius), keep_geom_type=True + ) ) - ) + else: + chunk_region_to_tile = region_to_rasterise # remove any subpixel polygons chunk_region_to_tile = chunk_region_to_tile[ chunk_region_to_tile.area @@ -1325,7 +1340,8 @@ def add_lidar( # Clip DEM to Catchment and ensure NaN outside region to rasterise dem = dem.rio.clip(self.catchment_geometry.catchment.geometry, drop=True) - dem = dem.rio.clip(clip_region.geometry, drop=False) + if clip_region.area.sum() > 0: # If area of 0 all will be NaN anyway + dem = dem.rio.clip(clip_region.geometry, drop=False) # If drop offshrore LiDAR ensure the foreshore values are 0 or negative if ( @@ -1715,13 +1731,16 @@ def add_coarse_dems( .count() .isnull() ) - no_value_mask &= ( - xarray.ones_like(self._dem.z) - .rio.clip( - self.catchment_geometry.land_and_foreshore.geometry, drop=False - ) - .notnull() - ) # Awkward as clip of a bool xarray doesn't work as expected + if self.catchment_geometry.land_and_foreshore.area.sum() > 0: + no_value_mask &= ( + xarray.ones_like(self._dem.z) + .rio.clip( + self.catchment_geometry.land_and_foreshore.geometry, drop=False + ) + .notnull() + ) # Awkward as clip of a bool xarray doesn't work as expected + else: + no_value_mask = xarray.zeros_like(self._dem.z) if not no_value_mask.any(): logging.info( f"No land areas greater than the cell buffer {buffer_cells}" diff --git a/src/geofabrics/geometry.py b/src/geofabrics/geometry.py index 920b2232..a310801c 100644 --- a/src/geofabrics/geometry.py +++ b/src/geofabrics/geometry.py @@ -70,28 +70,35 @@ def _set_up(self): # Clip and remove any sub pixel regions self._land = self._catchment.clip(self._full_land, keep_geom_type=True) self._land = self._land[self._land.area > self.resolution * self.resolution] + if self._land.area.sum() > 0: + self._foreshore_and_offshore = self.catchment.overlay( + self.land, how="difference" + ) + # Buffer and clip and remove any sub pixel regions + self._land_and_foreshore = geopandas.GeoDataFrame( + {"geometry": self.land.buffer(self.resolution * self.foreshore_buffer)}, + crs=self.crs["horizontal"], + ) + self._land_and_foreshore = self._catchment.clip( + self._land_and_foreshore, keep_geom_type=True + ) + self._land_and_foreshore = self._land_and_foreshore[ + self._land_and_foreshore.area > self.resolution * self.resolution + ] - self._foreshore_and_offshore = self.catchment.overlay( - self.land, how="difference" - ) + self._foreshore = self._land_and_foreshore.overlay(self.land, how="difference") - # Buffer and clip and remove any sub pixel regions - self._land_and_foreshore = geopandas.GeoDataFrame( - {"geometry": self.land.buffer(self.resolution * self.foreshore_buffer)}, - crs=self.crs["horizontal"], - ) - self._land_and_foreshore = self._catchment.clip( - self._land_and_foreshore, keep_geom_type=True - ) - self._land_and_foreshore = self._land_and_foreshore[ - self._land_and_foreshore.area > self.resolution * self.resolution - ] + self._offshore = self.catchment.overlay( + self.land_and_foreshore, how="difference" + ) + else: + # Assume no foreshore + self._foreshore_and_offshore = self.catchment + self._foreshore = self._land + self._land_and_foreshore = self._land + self._offshore = self.catchment - self._foreshore = self._land_and_foreshore.overlay(self.land, how="difference") - - self._offshore = self.catchment.overlay( - self.land_and_foreshore, how="difference" - ) + assert len(self._catchment) == 1, ( "The catchment is made of multiple separate polygons it must be a single " @@ -166,12 +173,15 @@ def land_and_foreshore_without_lidar(self, dense_extents: geopandas.GeoDataFrame self._assert_land_set() - if dense_extents is None: + if dense_extents is None or dense_extents.is_empty.all(): logging.warning( "In CatchmentGeometry dense extents are `None` so " " `land_and_foreshore_without_lidar` is returning `land_and_foreshore`" ) return self.land_and_foreshore + elif self.land_and_foreshore.area.sum() == 0: + # There is no land and foreshore region - return an empty dataframe + return self.land_and_foreshore # Clip to remove any offshore regions before doing a difference overlay. Drop # any sub-pixel polygons. land_and_foreshore_with_lidar = dense_extents.clip( @@ -191,7 +201,7 @@ def offshore_without_lidar(self, dense_extents: geopandas.GeoDataFrame): self._assert_land_set() - if dense_extents is None: + if dense_extents is None or dense_extents.is_empty.all(): logging.warning( "In CatchmentGeometry dense extents are `None` so " "`offshore_without_lidar` is returning `offshore`" @@ -219,12 +229,16 @@ def offshore_dense_data_edge(self, dense_extents: geopandas.GeoDataFrame): self._assert_land_set() - if dense_extents is None: + if dense_extents is None or dense_extents.is_empty.all(): logging.warning( "In CatchmentGeometry dense extents are `None` so " "`offshore_dense_data_edge` is returning `None`" ) - return None + return dense_extents + elif self.foreshore.area.sum(): + # If no foreshore just return the foreshore + return self.foreshore + # the foreshore and whatever lidar extents are offshore offshore_foreshore_dense_data_extents = geopandas.GeoDataFrame( { diff --git a/src/geofabrics/processor.py b/src/geofabrics/processor.py index 1a82b8e2..03dbdd01 100644 --- a/src/geofabrics/processor.py +++ b/src/geofabrics/processor.py @@ -895,6 +895,8 @@ def run(self): with cluster, distributed.Client(cluster) as client: print("Dask client:", client) print("Dask dashboard:", client.dashboard_link) + logging.info("Dask client:", client) + logging.info("Dask dashboard:", client.dashboard_link) # Load in LiDAR tiles self.raw_dem.add_lidar( @@ -1106,6 +1108,8 @@ def run(self): with cluster, distributed.Client(cluster) as client: print("Dask client:", client) print("Dask dashboard:", client.dashboard_link) + logging.info("Dask client:", client) + logging.info("Dask dashboard:", client.dashboard_link) # setup the hydrologically conditioned DEM generator self.hydrologic_dem = dem.HydrologicallyConditionedDem(