Support for stopbanks down and added in a test

src/geofabrics/dem.py

@@ -356,11 +356,14 @@ def save_dem(
             self._write_netcdf_conventions_in_place(dem, self.catchment_geometry.crs)
             if filename.suffix.lower() == ".nc":
                 if compression is not None:
-                    compression["grid_mapping"] = dem.encoding["grid_mapping"]
+                    encoding_keys = ("_FillValue", "dtype", "scale_factor", "add_offset", "grid_mapping") 
                     encoding = {}
                     for key in dem.data_vars:
-                        compression["dtype"] = dem[key].dtype
-                        encoding[key] = compression
+                        dem[key] = dem[key].astype(geometry.RASTER_TYPE)
+                        encoding[key] = {encoding_key: value for encoding_key, value in dem[key].encoding.items() if encoding_key in encoding_keys}
+                        if "dtype" not in encoding[key]:
+                            encoding[key]["dtype"] = dem[key].dtype
+                        encoding[key] = {**encoding[key], **compression}
                     dem.to_netcdf(
                         filename, format="NETCDF4", engine="netcdf4", encoding=encoding
                     )
@@ -1054,14 +1057,15 @@ def interpolate_ocean_bathymetry(self, bathy_contours, method="linear"):
 
     def clip_within_polygon(self, polygon_paths: list, label: str):
         """Clip existing DEM to remove areas within the polygons"""
+        crs = self.catchment_geometry.crs
         dem_bounds = geopandas.GeoDataFrame(
             geometry=[shapely.geometry.box(*self._dem.rio.bounds())],
-            crs=self.catchment_geometry.crs,
+            crs=crs["horizontal"],
         )
         clip_polygon = []
         for path in polygon_paths:
             clip_polygon.append(
-                geopandas.read_file(path).to_crs(self.catchment_geometry.crs)
+                geopandas.read_file(path).to_crs(crs["horizontal"])
             )
         clip_polygon = pandas.concat(clip_polygon).dissolve()
         clip_polygon = clip_polygon.clip(dem_bounds)

src/geofabrics/processor.py

@@ -1391,7 +1391,7 @@ def add_hydrological_features(
 
             if len(file_names) > 0:
                 hydrologic_dem.clip_within_polygon(
-                    polygons=file_names,
+                    polygon_paths=file_names,
                     label="masked feature",
                 )
                 temp_file = temp_folder / "dem_feature_masking.nc"

tests/test_stopbanks_waikanae/instruction.json

@@ -58,6 +58,7 @@
 			"extents": "catchment.geojson",
 			"raw_dem": "raw_dem.nc",
 			"stopbanks": [{"extents": "stopbanks/stopbank_polygon.geojson", "elevations": "stopbanks/stopbank_elevation.geojson"}],
+            "feature_masking": ["stopbank_masking_polygon.geojson"],
 			"result_dem": "test_dem.nc",
             "benchmark": "benchmark.nc"
 		},
@@ -71,7 +72,8 @@
 		"general": {
 			"z_labels": {"stopbanks": "z"},
 			"lidar_classifications_to_keep": [2, 9],
-			"interpolation": {"no_data": "linear"}
+			"interpolation": {"no_data": "linear"},
+            "compression": 7
 		}
 	}
 }

tests/test_stopbanks_waikanae/test_case.py

@@ -45,21 +45,28 @@ def setUpClass(cls):
             "datasets": {"vector": {"linz": {"key": linz_key}}}
         }
 
-        # create fake catchment boundary
+        crs = cls.instructions["dem"]["output"]["crs"]["horizontal"]
+
+        # create and save fake catchment boundary
         x0 = 1771850
         y0 = 5472250
         x1 = 1772150
         y1 = 5472600
         catchment = shapely.geometry.Polygon([(x0, y0), (x1, y0), (x1, y1), (x0, y1)])
-        catchment = geopandas.GeoSeries([catchment])
-        catchment = catchment.set_crs(
-            cls.instructions["dem"]["output"]["crs"]["horizontal"]
-        )
-
-        # save faked catchment boundary - used as land boundary as well
+        catchment = geopandas.GeoDataFrame(geometry=[catchment], crs=crs)
         catchment_file = cls.results_dir / "catchment.geojson"
         catchment.to_file(catchment_file)
 
+        # create and save stopbank masking polygon
+        x0 = 1771866
+        y0 = 5472488
+        x1 = 1771882
+        y1 = 5472568
+        catchment = shapely.geometry.Polygon([(x0, y0), (x1, y0), (x1, y1), (x0, y1)])
+        catchment = geopandas.GeoDataFrame(geometry=[catchment], crs=crs)
+        catchment_file = cls.results_dir / "stopbank_masking_polygon.geojson"
+        catchment.to_file(catchment_file)
+
         # Run pipeline - download files and generated DEM
         runner.from_instructions_dict(cls.instructions)
 
@@ -77,17 +84,21 @@ def test_result_windows(self):
         )
         with rioxarray.rioxarray.open_rasterio(file_path, masked=True) as test:
             test.load()
-        # Compare the generated and benchmark elevations
-        diff_array = (
-            test.z.data[~numpy.isnan(test.z.data)]
-            - benchmark.z.data[~numpy.isnan(benchmark.z.data)]
-        )
-        logging.info(f"DEM elevation diff is: {diff_array[diff_array != 0]}")
+
+        # Comparisons
         numpy.testing.assert_array_almost_equal(
             test.z.data,
             benchmark.z.data,
-            err_msg="The generated result_geofabric has different elevation data from "
-            "the benchmark_dem",
+            decimal=1,
+            err_msg="The generated test has significantly different elevation from the "
+            f"benchmark.",
+        )
+
+        numpy.testing.assert_array_almost_equal(
+            test.data_source.data,
+            benchmark.data_source.data,
+            err_msg="The generated test data_source layer differs from the "
+            f"benchmark.",
         )
 
         # explicitly free memory as xarray seems to be hanging onto memory
@@ -110,32 +121,20 @@ def test_result_linux(self):
         )
         with rioxarray.rioxarray.open_rasterio(file_path, masked=True) as test:
             test.load()
-        # Get data generated from LiDAR
-        lidar_mask = (test.data_source.data == 1) & (benchmark.data_source.data == 1)
 
-        # Compare the generated and benchmark elevations
-        lidar_diff = test.z.data[lidar_mask] - benchmark.z.data[lidar_mask]
+        # Comparisons
         numpy.testing.assert_array_almost_equal(
-            test.z.data[lidar_mask],
-            benchmark.z.data[lidar_mask],
-            decimal=6,
+            test.z.data,
+            benchmark.z.data,
+            decimal=1,
             err_msg="The generated test has significantly different elevation from the "
-            f"benchmark where there is LiDAR: {lidar_diff}",
-        )
-
-        diff_array = (
-            test.z.data[~numpy.isnan(test.z.data)]
-            - benchmark.z.data[~numpy.isnan(test.z.data)]
+            f"benchmark.",
         )
-        logging.info(f"DEM array diff is: {diff_array[diff_array != 0]}")
-        threshold = 10e-6
-        percent = 2.5
-        number_above_threshold = len(diff_array[numpy.abs(diff_array) > threshold])
-        self.assertTrue(
-            number_above_threshold < len(diff_array) * percent / 100,
-            f"More than {percent}% of DEM values differ by more than {threshold} on Linux test"
-            f" run: {diff_array[numpy.abs(diff_array) > threshold]} or "
-            f"{number_above_threshold / len(diff_array.flatten()) * 100}%",
+        numpy.testing.assert_array_almost_equal(
+            test.data_source.data,
+            benchmark.data_source.data,
+            err_msg="The generated test data_source layer differs from the "
+            f"benchmark.",
         )
 
         # explicitly free memory as xarray seems to be hanging onto memory