diff --git a/src/geofabrics/dem.py b/src/geofabrics/dem.py index bd795816..8420911b 100644 --- a/src/geofabrics/dem.py +++ b/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) diff --git a/src/geofabrics/processor.py b/src/geofabrics/processor.py index 53d805d7..640c0e7b 100644 --- a/src/geofabrics/processor.py +++ b/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" diff --git a/tests/test_stopbanks_waikanae/data/benchmark.nc b/tests/test_stopbanks_waikanae/data/benchmark.nc index a6a577f6..62f6acee 100644 --- a/tests/test_stopbanks_waikanae/data/benchmark.nc +++ b/tests/test_stopbanks_waikanae/data/benchmark.nc @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:b44d1b9f74d529e521fddc4a2d36f9edfe00043ec400f88136ebbf2dbdd6159f -size 38537 +oid sha256:802ba5c1f6b6fef96a2ce4d34d2dd2d386f5cf942b45f14834d4a8a2e0a15ae2 +size 38699 diff --git a/tests/test_stopbanks_waikanae/instruction.json b/tests/test_stopbanks_waikanae/instruction.json index b857030b..f462b2c4 100644 --- a/tests/test_stopbanks_waikanae/instruction.json +++ b/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 } } } diff --git a/tests/test_stopbanks_waikanae/test_case.py b/tests/test_stopbanks_waikanae/test_case.py index 5be76e3a..90585eb2 100644 --- a/tests/test_stopbanks_waikanae/test_case.py +++ b/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