Skip to content

Commit

Permalink
bug-fixes identified in cylc suite (#222)
Browse files Browse the repository at this point in the history 
* make robust to no ocean bathy data within the download area

* fixup: Format Python code with Black

* Add file caching to the roughness step.

* Add optional filtering of estimated waterways by OSM id.

* Update README.md

* make OSM id filter optional

* update version

* Updated test to include a filter for OSM

---------

Co-authored-by: github-actions <[email protected]>
  • Loading branch information
@rosepearson
rosepearson and github-actions authored Dec 1, 2023
1 parent 21b82c3 commit 22bdede
Show file tree
Hide file tree
Showing 8 changed files with 209 additions and 141 deletions.
4 changes: 2 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

The `geofabrics` package includes routines and classes for combining point (i.e. LiDAR), vector (i.e. catchment of interest, infrastructure), and raster (i.e. reference DEM) to generate a hydrologically conditioned raster.

A peer-reviewed journal article of the package and methodogy can be found at: [https://www.sciencedirect.com/science/article/pii/S1364815223002281](https://www.sciencedirect.com/science/article/pii/S1364815223002281)
A peer-reviewed journal article of the package and methodogy can be found at: [https://www.sciencedirect.com/science/article/pii/S1364815223002281](https://www.sciencedirect.com/science/article/pii/S1364815223002281). This is sadly not Open Access. Please email [email protected] if you would like a copy of the article.

## Installation
`geofabrics` is avaliable on [conda-forge](https://anaconda.org/conda-forge/geofabrics) and [PyPI](https://pypi.org/project/geofabrics/). Conda is recommended due to difficulties installing geopandas (a dependency) with pip on Windows. See the [Wiki Install Instructions](https://github.com/rosepearson/GeoFabrics/wiki/Package-Install-Instructions) for more information.
Expand Down Expand Up @@ -74,7 +74,7 @@ Maintainer: Rose Pearson @rosepearson [email protected]
<a href="https://github.com/AliceHarang">
<img src="https://avatars.githubusercontent.com/u/82012970?v=4" width="100;" alt="AliceHarang"/>
<br />
<sub><b>Null</b></sub>
<sub><b>Alice Harang</b></sub>
</a>
</td></tr>
<tr>
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"

[project]
name = "geofabrics"
version = "1.1.3"
version = "1.1.4"
description = "A package for creating geofabrics for flood modelling."
readme = "README.md"
authors = [{ name = "Rose pearson", email = "[email protected]" }]
Expand Down
206 changes: 114 additions & 92 deletions src/geofabrics/dem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1166,6 +1166,104 @@ def _add_lidar_no_chunking(
"_add_lidar_no_chunking must be instantiated in the " "child class"
)

def _load_dem(self, filename: pathlib.Path, chunk_size: int):
"""Load in and replace the DEM with a previously cached version."""
dem = rioxarray.rioxarray.open_rasterio(
filename,
masked=True,
parse_coordinates=True,
chunks={"x": chunk_size, "y": chunk_size},
)
dem = dem.squeeze("band", drop=True)
self._write_netcdf_conventions_in_place(dem, self.catchment_geometry.crs)
self._dem = dem

if "no_values_mask" in self._dem.keys():
self._dem["no_values_mask"] = self._dem.no_values_mask.astype(bool)
if "data_source" in self._dem.keys():
self._dem["data_source"] = self._dem.data_source.astype(
geometry.RASTER_TYPE
)
if "lidar_source" in self._dem.keys():
self._dem["lidar_source"] = self._dem.lidar_source.astype(
geometry.RASTER_TYPE
)
if "z" in self._dem.keys():
self._dem["z"] = self._dem.z.astype(geometry.RASTER_TYPE)

def save_dem(
self,
filename: pathlib.Path,
dem: xarray.Dataset,
):
"""Save the DEM to a netCDF file."""

# Save the file
try:
dem.to_netcdf(
filename,
format="NETCDF4",
engine="netcdf4",
)
# Close the DEM
dem.close()
except (Exception, KeyboardInterrupt) as caught_exception:
pathlib.Path(filename).unlink()
logging.info(
f"Caught error {caught_exception} and deleting"
"partially created netCDF output "
f"{filename} before re-raising error."
)
raise caught_exception

def _save_and_load_dem(
self,
filename: pathlib.Path,
chunk_size: int,
no_values_mask: bool,
buffer_cells: int = None,
):
"""Update the saved file cache for the DEM as a netCDF file. The bool
no_data_layer may optionally be included."""

# Logging update
logging.info(
"In LidarBase._save_and_load_dem saving NetCDF file to "
f"{filename} with no_values_mask={no_values_mask}"
)
# Get the DEM from the property call
dem = self._dem
# Create mask if specified
if no_values_mask:
if self.catchment_geometry.land_and_foreshore.area.sum() > 0:
no_value_mask = (
dem.z.rolling(
dim={"x": buffer_cells * 2 + 1, "y": buffer_cells * 2 + 1},
min_periods=1,
center=True,
)
.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
else:
no_value_mask = xarray.zeros_like(self._dem.z)
dem["no_values_mask"] = no_value_mask
dem.no_values_mask.rio.write_crs(
self.catchment_geometry.crs["horizontal"], inplace=True
)

# Save the DEM with the no_values_layer
self.save_dem(filename=filename, dem=dem)
# Load in the temporarily saved DEM
self._load_dem(filename=filename, chunk_size=chunk_size)


class RawDem(LidarBase):
"""A class to manage the creation of a 'raw' DEM from LiDAR tiles, and/or a
Expand Down Expand Up @@ -1427,10 +1525,11 @@ def add_lidar(

# Save a cached copy of DEM to temporary memory cache
logging.info("In dem.add_lidar - write out temp raw DEM to netCDF")
self._save_and_load_dem_with_no_values_mask(
self._save_and_load_dem(
filename=self.temp_folder / "raw_lidar.nc",
buffer_cells=buffer_cells,
chunk_size=chunk_size,
no_values_mask=True,
)

def _add_tiled_lidar_chunked(
Expand Down Expand Up @@ -1915,11 +2014,12 @@ def dask_interpolation(y, x):
logging.info(
"In dem.add_coarse_dems - write out temp raw DEM to netCDF"
)
self._save_and_load_dem_with_no_values_mask(
self._save_and_load_dem(
filename=self.temp_folder
/ f"raw_dem_{coarse_dem_path.stem}.nc",
buffer_cells=buffer_cells,
chunk_size=chunk_size,
no_values_mask=True,
)
logging.info(
f"In dem.add_coarse_dems - remove previous cached file {previous_cached_file}"
Expand All @@ -1929,96 +2029,6 @@ def dask_interpolation(y, x):
self.temp_folder / f"raw_dem_{coarse_dem_path.stem}.nc"
)

def _load_dem(self, filename: pathlib.Path, chunk_size: int):
"""Load in and replace the DEM with a previously cached version."""
dem = rioxarray.rioxarray.open_rasterio(
filename,
masked=True,
parse_coordinates=True,
chunks={"x": chunk_size, "y": chunk_size},
)
dem = dem.squeeze("band", drop=True)
self._write_netcdf_conventions_in_place(dem, self.catchment_geometry.crs)
self._dem = dem

if "no_values_mask" in self._dem.keys():
self._dem["no_values_mask"] = self._dem.no_values_mask.astype(bool)
if "data_source" in self._dem.keys():
self._dem["data_source"] = self._dem.data_source.astype(
geometry.RASTER_TYPE
)
if "lidar_source" in self._dem.keys():
self._dem["lidar_source"] = self._dem.lidar_source.astype(
geometry.RASTER_TYPE
)
if "z" in self._dem.keys():
self._dem["z"] = self._dem.z.astype(geometry.RASTER_TYPE)

def save_dem(
self,
filename: pathlib.Path,
dem: xarray.Dataset,
):
"""Save the DEM to a netCDF file."""

# Save the file
try:
dem.to_netcdf(
filename,
format="NETCDF4",
engine="netcdf4",
)
# Close the DEM
dem.close()
except (Exception, KeyboardInterrupt) as caught_exception:
pathlib.Path(filename).unlink()
logging.info(
f"Caught error {caught_exception} and deleting"
"partially created netCDF output "
f"{filename} before re-raising error."
)
raise caught_exception

def _save_and_load_dem_with_no_values_mask(
self,
filename: pathlib.Path,
buffer_cells: int,
chunk_size: int,
):
"""Update the saved file cache for the DEM as a netCDF file. The no_data_layer of bol values may
optionally be included."""

# Get the DEM from the property call
dem = self._dem
if self.catchment_geometry.land_and_foreshore.area.sum() > 0:
no_value_mask = (
dem.z.rolling(
dim={"x": buffer_cells * 2 + 1, "y": buffer_cells * 2 + 1},
min_periods=1,
center=True,
)
.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
else:
no_value_mask = xarray.zeros_like(self._dem.z)
dem["no_values_mask"] = no_value_mask
dem.no_values_mask.rio.write_crs(
self.catchment_geometry.crs["horizontal"], inplace=True
)

# Save the DEM with the no_values_layer
self.save_dem(filename=filename, dem=dem)
# Load in the temporarily saved DEM
self._load_dem(filename=filename, chunk_size=chunk_size)


class RoughnessDem(LidarBase):
"""A class to add a roughness (zo) layer to a hydrologically conditioned DEM.
Expand Down Expand Up @@ -2047,6 +2057,7 @@ def __init__(
self,
catchment_geometry: geometry.CatchmentGeometry,
hydrological_dem_path: typing.Union[str, pathlib.Path],
temp_folder: pathlib.Path,
interpolation_method: str,
default_values: dict,
drop_offshore_lidar: dict,
Expand Down Expand Up @@ -2085,6 +2096,7 @@ def __init__(
self.catchment_geometry.catchment.geometry, drop=True
)

self.temp_folder = temp_folder
self.interpolation_method = interpolation_method
self.default_values = default_values
self.drop_offshore_lidar = drop_offshore_lidar
Expand Down Expand Up @@ -2160,6 +2172,11 @@ def add_lidar(
chunk_size=chunk_size,
metadata=metadata,
)
self._save_and_load_dem(
filename=self.temp_folder / "raw_lidar_zo.nc",
chunk_size=chunk_size,
no_values_mask=False,
)
# Set roughness where water
self._dem["zo"] = self._dem.zo.where(
self._dem.data_source != self.SOURCE_CLASSIFICATION["ocean bathymetry"],
Expand Down Expand Up @@ -2414,6 +2431,11 @@ def _add_roughness_to_data_set(
A list of roughnesses over the x, and y coordiantes for each dataset.
"""

logging.info(
"In RoughnessDem._add_roughness_to_data_set creating and "
"adding the Zo layer to the dataset."
)

# Create a DataArray of zo
zos = []
for zo_array in roughnesses:
Expand Down
20 changes: 18 additions & 2 deletions src/geofabrics/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -512,12 +512,14 @@ class EstimatedBathymetryPoints:
DEPTH_LABEL = "elevation"
BANK_HEIGHT_LABEL = "bank_height"
WIDTH_LABEL = "width"
OSM_ID = "OSM_id"

def __init__(
self,
points_files: list,
polygon_files: list,
catchment_geometry: CatchmentGeometry,
filter_osm_ids: list = [],
z_labels: list = None,
):
self.catchment_geometry = catchment_geometry
Expand All @@ -526,9 +528,15 @@ def __init__(
self._points = None
self._polygon = None

self._set_up(points_files, polygon_files, z_labels)
self._set_up(points_files, polygon_files, z_labels, filter_osm_ids)

def _set_up(self, points_files: list, polygon_files: list, z_labels: list):
def _set_up(
self,
points_files: list,
polygon_files: list,
z_labels: list,
filter_osm_ids: list,
):
"""Load point and polygon files and concatentate and clip to the catchment."""

if len(points_files) != len(polygon_files):
Expand Down Expand Up @@ -560,6 +568,8 @@ def _set_up(self, points_files: list, polygon_files: list, z_labels: list):
columns_i.append(self.BANK_HEIGHT_LABEL)
if self.WIDTH_LABEL in points_i.columns:
columns_i.append(self.WIDTH_LABEL)
if self.OSM_ID in points_i.columns:
columns_i.append(self.OSM_ID)
# Only add the points and polygons if their are points
if len(points_i) > 0:
points_i = points_i[columns_i]
Expand Down Expand Up @@ -587,6 +597,12 @@ def _set_up(self, points_files: list, polygon_files: list, z_labels: list):
points = points.clip(self.catchment_geometry.catchment, keep_geom_type=True)
points = points.reset_index(drop=True)

# Filter if there are OSM IDs specified to filter
if len(filter_osm_ids) > 0:
for osm_id in filter_osm_ids:
points = points[points[self.OSM_ID] != osm_id]
polygon = polygon[polygon[self.OSM_ID] != osm_id]

# Set to class members
self._points = points
self._polygon = polygon
Expand Down
Loading

0 comments on commit 22bdede

Please sign in to comment.