271 bugfix min value roughness where interpolated (#272)

* correct order of min/max enforced Zo so that only applied after all values set

* Fix missing clip in ElevationPoints function

* Updated stopbanks to actually use width as width and not as radius

* minor update to stopbanks benchmark

* updated benchmark given minor numerical changes

* updated version information

pyproject.toml

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "geofabrics"
-version = "1.1.25"
+version = "1.1.26"
 description = "A package for creating geofabrics for flood modelling."
 readme = "README.md"
 authors = [{ name = "Rose pearson", email = "rose.pearson@niwa.co.nz" }]

src/geofabrics/dem.py

@@ -2863,6 +2863,18 @@ def add_lidar(
             if numpy.isnan(self._dem.zo.data).any():
                 self._dem["zo"] = self._dem.zo.rio.interpolate_na(method="nearest")
 
+        # Ensure roughness values are bounded by the defaults
+        if self.default_values["minimum"] is not None:
+            self._dem["zo"] = self._dem.zo.where(
+                self._dem.zo > self.default_values["minimum"],
+                self.default_values["minimum"],
+            )
+        if self.default_values["maximum"] is not None:
+            self._dem["zo"] = self._dem.zo.where(
+                self._dem.zo < self.default_values["maximum"],
+                self.default_values["maximum"],
+            )
+
         mask = clip_mask(
             self._dem.z, self.catchment_geometry.catchment.geometry, self.chunk_size
         )
@@ -3156,17 +3168,6 @@ def _add_roughness_to_data_set(
             zo = rioxarray.merge.merge_arrays(zos, method="first")
         # Resize zo to share the same dimensions at the DEM
         self._dem["zo"] = zo.sel(x=self._dem.x, y=self._dem.y, method="nearest")
-        # Ensure roughness values are bounded by the defaults
-        if self.default_values["minimum"] is not None:
-            self._dem["zo"] = self._dem.zo.where(
-                self._dem.zo > self.default_values["minimum"],
-                self.default_values["minimum"],
-            )
-        if self.default_values["maximum"] is not None:
-            self._dem["zo"] = self._dem.zo.where(
-                self._dem.zo < self.default_values["maximum"],
-                self.default_values["maximum"],
-            )
 
         # update metadata
         history = self._dem.attrs["history"]

src/geofabrics/geometry.py

@@ -644,6 +644,7 @@ def _set_up(
             crs=polygon_list[0].crs,
         )
         polygon = polygon.to_crs(self.catchment_geometry.crs["horizontal"])
+        polygon = polygon.clip(self.catchment_geometry.catchment, keep_geom_type=True)
         points = points.clip(polygon.buffer(0), keep_geom_type=True)
         points = points.clip(self.catchment_geometry.catchment, keep_geom_type=True)
         points = points.reset_index(drop=True)
@@ -858,7 +859,8 @@ class RiverMouthFan:
         cross_section_spacing  The spacing in (m) of the sampled cross sections.
         aligned_channel_file  Thefile name for the aligned river channel file.
         river_bathymetry_file  The file name for the river bathymetry values.
-        ocean_contour_file  The file name for the ocean contours.
+        ocean_contour_file  The file name for the ocean contours. Depths are positive
+        ocean_points_file The file name for the ocean points. Elevations are negative.
         ocean_contour_depth_label  The column label for the depth values.
     """
 
@@ -1150,12 +1152,13 @@ def polygon_and_bathymetry(self):
         )
 
         if self.ocean_points_file is not None:
-            length, end_depth = self._get_distance_to_nearest_ocean_point_in_fan(
+            length, end_elevation = self._get_distance_to_nearest_ocean_point_in_fan(
                 river_mouth_z=max(river_mouth_elevations),
                 fan=fan_polygon,
                 mouth_point=mouth_point,
                 depth_multiplier=2
             )
+            end_depth = -1 * end_elevation
             fan_polygon = self._fixed_length_fan_polygon(
                 river_mouth_width=river_mouth_width,
                 mouth_point=mouth_point,

src/geofabrics/processor.py

@@ -3943,7 +3943,7 @@ def sample(geometry):
         points = points.sort_index(ascending=True).explode(
             ignore_index=False, index_parts=True, column="geometry"
         )
-        points["polygons"] = points.buffer(points["width"].to_numpy())
+        points["polygons"] = points.buffer(points["width"].to_numpy() / 2)
         # Sample maximum elevation in polygon around each point
         points["z"] = numpy.nan
         for index, rows in points.groupby(level=0):
@@ -3965,7 +3965,7 @@ def sample(geometry):
                 "estimate their creast elevations."
             )
         # Create, filter to remove NaN areas and save overall polygon
-        stopbanks["polygon"] = stopbanks.buffer(stopbanks["width"].to_numpy())
+        stopbanks["polygon"] = stopbanks.buffer(stopbanks["width"].to_numpy() / 2)
         stopbanks = stopbanks.sort_index(ascending=True)
         stopbanks = stopbanks[nan_filter]
         stopbanks.set_geometry("polygon", drop=True)[["geometry"]].to_file(polygon_file)
@@ -3985,7 +3985,7 @@ def create_dem(self, stopbanks: geopandas.GeoDataFrame) -> xarray.Dataset:
                     key="stopbank_polygon", index=index
                 )
                 stopbank_polygon = geopandas.GeoDataFrame(
-                    geometry=[row.geometry.buffer(row.width)], crs=stopbanks.crs
+                    geometry=[row.geometry.buffer(row.width / 2)], crs=stopbanks.crs
                 )
                 stopbank_polygon.to_file(stopbank_polygon_file)
 

src/geofabrics/version.py

@@ -3,4 +3,4 @@
 Contains the package version information
 """
 
-__version__ = "1.1.25"
+__version__ = "1.1.26"

tests/test_stopbanks_waikanae/instruction.json

@@ -32,7 +32,7 @@
 			"lidar_classifications_to_keep": [2, 9]
 		},
 		"stopbanks": {
-			"width": 6,
+			"width": 12,
 		"source": "file"
 		}
 	},