bug-fixes identified in cylc suite

README.md

@@ -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.
@@ -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>

pyproject.toml

@@ -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]" }]

src/geofabrics/dem.py

@@ -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
@@ -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(
@@ -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}"
@@ -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.
@@ -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,
@@ -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
@@ -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"],
@@ -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:

src/geofabrics/geometry.py

@@ -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
@@ -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):
@@ -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]
@@ -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