Skip to content

Commit

Permalink
Added support for only offshore data
Browse files Browse the repository at this point in the history
  • Loading branch information
rosepearson committed Oct 26, 2023
1 parent 7cac6f0 commit 96184e8
Show file tree
Hide file tree
Showing 3 changed files with 72 additions and 35 deletions.
43 changes: 31 additions & 12 deletions src/geofabrics/dem.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 (
Expand Down Expand Up @@ -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}"
Expand Down
60 changes: 37 additions & 23 deletions src/geofabrics/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 "
Expand Down Expand Up @@ -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(
Expand All @@ -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`"
Expand Down Expand Up @@ -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(
{
Expand Down
4 changes: 4 additions & 0 deletions src/geofabrics/processor.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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(
Expand Down

0 comments on commit 96184e8

Please sign in to comment.