From 21b82c3a5804cd4e5b5238fc21fbf7d97e4e4744 Mon Sep 17 00:00:00 2001 From: Rose Pearson Date: Wed, 15 Nov 2023 20:10:56 +1300 Subject: [PATCH] 220 - bugfixes - saving temp files (#221) * ensure robust to no LiDAR dataset * fixup: Format Python code with Black * Fix to no_value_mask error * Update package version --------- Co-authored-by: github-actions --- pyproject.toml | 2 +- src/geofabrics/dem.py | 74 +++++++++++++++++++++++++-------------- src/geofabrics/version.py | 2 +- 3 files changed, 50 insertions(+), 28 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 42bf0a5c..43d70094 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta" [project] name = "geofabrics" -version = "1.1.2" +version = "1.1.3" description = "A package for creating geofabrics for flood modelling." readme = "README.md" authors = [{ name = "Rose pearson", email = "rose.pearson@niwa.co.nz" }] diff --git a/src/geofabrics/dem.py b/src/geofabrics/dem.py index 35fd7702..9b97eb6d 100644 --- a/src/geofabrics/dem.py +++ b/src/geofabrics/dem.py @@ -1335,7 +1335,38 @@ def add_lidar( } # Don't use dask delayed if there is no chunking - if chunk_size is None: + if len(lidar_datasets_info) == 0: + # Create an empty dataset as no LiDAR + logging.warning("No LiDAR dataset. Creating an empty raw DEM dataset.") + bounds = self.catchment_geometry.catchment.geometry.bounds + resolution = self.catchment_geometry.resolution + x = numpy.arange( + numpy.ceil(bounds.minx.min() / resolution) * resolution, + numpy.ceil(bounds.maxx.max() / resolution) * resolution, + resolution, + dtype=raster_options["raster_type"], + ) + y = numpy.arange( + numpy.ceil(bounds.maxy.max() / resolution) * resolution, + numpy.ceil(bounds.miny.min() / resolution) * resolution, + -resolution, + dtype=raster_options["raster_type"], + ) + metadata["instructions"]["dataset_mapping"]["lidar"] = { + "no LiDAR": self.SOURCE_CLASSIFICATION["no data"] + } + elevations = { + "no LiDAR": dask.array.empty( + shape=(len(y), len(x)), dtype=raster_options["raster_type"] + ) + } + dem = self._create_data_set( + x=x, + y=y, + elevations=elevations, + metadata=metadata, + ) + elif chunk_size is None: dem = self._add_lidar_no_chunking( lidar_datasets_info=lidar_datasets_info, options=raster_options, @@ -1396,9 +1427,8 @@ 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._load_and_save_dem( + self._save_and_load_dem_with_no_values_mask( filename=self.temp_folder / "raw_lidar.nc", - no_values_mask=True, buffer_cells=buffer_cells, chunk_size=chunk_size, ) @@ -1750,7 +1780,6 @@ def add_coarse_dems( for coarse_dem_path in coarse_dem_paths: # Check if any areas (on land and foreshore) still without values - exit if none no_value_mask = self._dem.no_values_mask - self._dem = self._dem.drop_vars("no_values_mask") if not no_value_mask.any(): logging.info( f"No land areas greater than the cell buffer {buffer_cells}" @@ -1886,10 +1915,9 @@ def dask_interpolation(y, x): logging.info( "In dem.add_coarse_dems - write out temp raw DEM to netCDF" ) - self._load_and_save_dem( + self._save_and_load_dem_with_no_values_mask( filename=self.temp_folder / f"raw_dem_{coarse_dem_path.stem}.nc", - no_values_mask=True, buffer_cells=buffer_cells, chunk_size=chunk_size, ) @@ -1951,10 +1979,9 @@ def save_dem( ) raise caught_exception - def _load_and_save_dem( + def _save_and_load_dem_with_no_values_mask( self, filename: pathlib.Path, - no_values_mask: bool, buffer_cells: int, chunk_size: int, ): @@ -1963,7 +1990,7 @@ def _load_and_save_dem( # Get the DEM from the property call dem = self._dem - 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}, @@ -1973,24 +2000,19 @@ def _load_and_save_dem( .count() .isnull() ) - if self.catchment_geometry.land_and_foreshore.area.sum() > 0: - 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 - ) - dem.no_values_mask.rio.write_nodata(numpy.nan, encoded=True, inplace=True) + 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: - if "no_values_mask" in dem: - dem = dem.drop("no_values_mask") + 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) diff --git a/src/geofabrics/version.py b/src/geofabrics/version.py index 5352b8d5..46eefa51 100644 --- a/src/geofabrics/version.py +++ b/src/geofabrics/version.py @@ -3,4 +3,4 @@ Contains the package version information """ -__version__ = "1.1.2" +__version__ = "1.1.3"