Initial code changes to add support for roughness from a landuse map

pyproject.toml

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "geofabrics"
-version = "1.1.29"
+version = "1.1.30"
 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

@@ -2915,20 +2915,27 @@ def dem(self):
                 self._dem.data_source != self.SOURCE_CLASSIFICATION["waterways"],
                 self.default_values["waterways"],
             )
-        print(self.default_values)
         if self.default_values["lakes"] is not None:
             self._dem["zo"] = self._dem.zo.where(
                 self._dem.data_source != self.SOURCE_CLASSIFICATION["lakes"],
                 self.default_values["lakes"],
             )
+        # Set roughness where land and no LiDAR or landuse data
+        if self.default_values["land"] is not None:
+            roi_mask = clip_mask(
+                self._dem.zo,
+                self.catchment_geometry.land_and_foreshore.geometry,
+                self.chunk_size,
+            )
+            notnull_mask = self._dem.zo.isnull()
+            self._dem["zo"] = self._dem.zo.where(
+                ~(roi_mask & notnull_mask),
+                self.default_values["land"],
+            )
 
-        # Set roughness where land and no LiDAR
-        self._dem["zo"] = self._dem.zo.where(
-            self._dem.data_source != self.SOURCE_CLASSIFICATION["coarse DEM"],
-            self.default_values["land"],
-        )  # or LiDAR with no roughness estimate
         # Ensure the defaults are re-added
         self._write_netcdf_conventions_in_place(self._dem, self.catchment_geometry.crs)
+
         # Interpolate any missing roughness values
         if self.interpolation_method is not None:
             self._dem["zo"] = self._dem.zo.rio.interpolate_na(
@@ -3011,7 +3018,33 @@ def add_lidar(
             )
         return status
 
-    def add_roads(self, roads_polygon: dict):
+    def add_landuse(self, landuse_polygon: geopandas.GeoDataFrame):
+        """Set landuse to roughness values.
+
+        Parameters
+        ----------
+
+        landuse_polygon
+            Dataframe with polygon and associated roughness values
+        """
+
+        # Set unpaved roads
+        for land_category in landuse_polygon["landcover"].unique():
+            mask = clip_mask(
+                self._dem.z,
+                landuse_polygon[landuse_polygon["landcover"] == land_category].geometry,
+                self.chunk_size,
+            )
+            mask &= self._dem.zo.isnull()  # Only set roughness where not set by LiDAR
+            roughness_value = landuse_polygon[
+                landuse_polygon["landcover"] == land_category
+            ]["roughness"].values[0]
+            self._dem["zo"] = self._dem.zo.where(~mask, roughness_value)
+            self._write_netcdf_conventions_in_place(
+                self._dem, self.catchment_geometry.crs
+            )
+
+    def add_roads(self, roads_polygon: geopandas.GeoDataFrame):
         """Set roads to paved and unpaved roughness values.
 
         Parameters

src/geofabrics/processor.py

@@ -696,11 +696,15 @@ def get_lidar_datasets_info(self) -> dict:
 
             # download the specified datasets from the data service - then get the
             # local file path
-            search_polygon = (
-                self.catchment_geometry.catchment
-                if self.catchment_geometry is not None
-                else None
-            )
+            if self.catchment_geometry is None:
+                search_polygon = None
+            else:
+                search_polygon = geopandas.GeoDataFrame(
+                    self.catchment_geometry.catchment.buffer(
+                        self.catchment_geometry.resolution / numpy.sqrt(2)
+                    ),
+                    columns=["geometry"],
+                )
             self.lidar_fetcher = geoapis.lidar.OpenTopography(
                 cache_path=self.get_instruction_path("downloads") / "lidar",
                 search_polygon=search_polygon,
@@ -1773,12 +1777,13 @@ def get_roughness_instruction(self, key: str):
                     "motorway": 12,
                 },
             },
+            "landuse": {"drop_offshore": True},
         }
 
         if "roughness" in self.instructions and key in self.instructions["roughness"]:
             roughness_instruction = self.instructions["roughness"][key]
             # ensure all default keys included if a dictionary
-            if key in ["default_values", "parameters"]:
+            if key in defaults.keys() and isinstance(roughness_instruction, dict):
                 for sub_key in defaults[key]:
                     if sub_key not in roughness_instruction:
                         roughness_instruction[sub_key] = defaults[key][sub_key]
@@ -1794,6 +1799,87 @@ def get_roughness_instruction(self, key: str):
                 " does not have a default value."
             )
 
+    def load_landuse_lris(self) -> bool:
+        """Download a Landuse map from LRISwithin the catchment BBox."""
+
+        defaults = {
+            "landuse_to_roughness": "landuse_to_roughness.geojson",
+        }
+        landuse_to_roughness_path = self.get_instruction_path(
+            "landuse_to_roughness", defaults=defaults
+        )
+
+        if landuse_to_roughness_path.is_file():
+            landuse_polygon = geopandas.read_file(landuse_to_roughness_path)
+            if landuse_polygon.area.sum() == 0:
+                message = (
+                    "Warning zero area landuse polygon provided. Will ignore. "
+                    f"Please check {landuse_to_roughness_path} if unexpected."
+                )
+                self.logger.warning(message)
+                return landuse_polygon
+            if "roughness" not in landuse_polygon.columns:
+                message = (
+                    "No roughnesses defined in the landuse polygon file. This is "
+                    f"required. Please check {landuse_to_roughness_path} and add."
+                )
+                self.logger.error(message)
+                raise ValueError(message)
+            return landuse_polygon
+
+        else:  # Download from LRIS
+            if not self.check_vector_or_raster("landuse", "vector"):
+                message = (
+                    "No roughnesses landuse dataset specified. This is required"
+                    f" if using `landuse`. Please as in the instruction file."
+                )
+                self.logger.error(message)
+                raise ValueError(message)
+            landuse_paths = self.get_vector_or_raster_paths(
+                "landuse", "vector", required=True
+            )
+            if len(landuse_paths) > 1:
+                self.logger.warning(
+                    f"{len(landuse_paths)} landuse datasets provided. "
+                    f"Specficially {landuse_paths}. Only consider the "
+                    "first if multiple."
+                )
+            landuse = geopandas.read_file(landuse_paths[0])
+
+            # Standardise columns and add rougness values
+            landuse_instructions = self.get_roughness_instruction("landuse")
+            if landuse_instructions["landcover"] not in landuse.columns:
+                message = (
+                    f"Column name {landuse_instructions['landcover']} not in "
+                    f"landuse dataset. Please check the instruction file and the "
+                    f"dataset. Dataset columns names are: {landuse.columns}."
+                )
+                self.logger.error(message)
+                raise ValueError(message)
+            # Remove any landuse classes without a specified roughness value
+            landuse = landuse[
+                landuse[landuse_instructions["landcover"]].isin(
+                    landuse_instructions["classes_to_roughness"].keys()
+                )
+            ]
+            landuse["roughness"] = landuse[landuse_instructions["landcover"]].map(
+                landuse_instructions["classes_to_roughness"]
+            )
+
+            if landuse_instructions["drop_offshore"]:  # Clip to land
+                landuse = landuse.clip(self.catchment_geometry.land).sort_index(
+                    ascending=True
+                )
+
+            landuse.rename(
+                columns={landuse_instructions["landcover"]: "landcover"}, inplace=True
+            )
+            landuse = landuse[["geometry", "roughness", "landcover"]]
+
+            # Save files
+            landuse.to_file(landuse_to_roughness_path)
+        return landuse
+
     def load_roads_osm(self) -> bool:
         """Download OpenStreetMap roads within the catchment BBox."""
 
@@ -1956,6 +2042,12 @@ def run(self):
         else:
             roads = None
 
+        landuse = self.get_roughness_instruction("landuse")
+        if "source" in landuse and "lris" in landuse["source"]:
+            landuse = self.load_landuse_lris()
+        else:
+            landuse = None
+
         # Create folder for caching raw DEM files during DEM generation
         temp_folder = self.setup_temp_folder()
         cached_file = temp_folder / "not_yet_created_file"
@@ -2012,6 +2104,19 @@ def run(self):
                 self.logger.info(f"Save temp raw DEM to netCDF: {cached_file}")
                 roughness_dem.save_and_load_dem(cached_file)
 
+            # If landuse save temp then add in the landuse
+            if landuse is not None and landuse.area.sum() > 0:
+
+                # Add landuse to roughness
+                roughness_dem.add_landuse(landuse_polygon=landuse)
+
+                # cache the results
+                temp_file = temp_folder / "geofabric_added_landuse.nc"
+                self.logger.info(f"Save geofabric with landuse to netCDF: {temp_file}")
+                roughness_dem.save_and_load_dem(temp_file)
+                self.clean_cached_file(cached_file)
+                cached_file = temp_file
+
             # If roads save temp then add in the roads
             if roads is not None and roads.area.sum() > 0:
 

src/geofabrics/version.py

@@ -3,4 +3,4 @@
 Contains the package version information
 """
 
-__version__ = "1.1.29"
+__version__ = "1.1.30"

tests/test_add_patches_ngaruroro/instruction.json

@@ -20,7 +20,7 @@
 		"patches": ["../patch_1.nc", "../patch_2.tif"],
         "raw_dem": "../initial_dem.nc",
 		"result_dem": "test_dem.nc",
-		"benchmark_dem": "benchmark_dem.nc"
+		"benchmark_dem": "benchmark.nc"
 	},
 	"patch": {
 		"layers": "z"

tests/test_dem_generation_local_1/instruction.json

@@ -17,7 +17,7 @@
 		"coarse_dems": ["coarse_dem.nc"],
 		"ocean_contours": ["bathymetry.geojson"],
 		"result_dem": "test_dem.tif",
-		"benchmark_dem": {"z": "benchmark_dem_z.tif", "data_source": "benchmark_dem_data_source.tif", "lidar_source": "benchmark_dem_lidar_source.tif"}
+		"benchmark_dem": {"z": "benchmark_z.tif", "data_source": "benchmark_data_source.tif", "lidar_source": "benchmark_lidar_source.tif"}
 	},
 	"general": {
 		"drop_offshore_lidar": true,

tests/test_many_stages_waikanae/instruction.json

@@ -194,32 +194,34 @@
 			"benchmark": "benchmark.nc"
 		},
 		"datasets": {
-			"lidar": {
-				"open_topography": {
-					"Wellington_2013": {
-						"crs": {
-							"horizontal": 2193,
-							"vertical": 7839
-						}
-					}
-				}
-			},
 			"vector": {
 				"linz": {
 					"land": {
 						"layers": [51153]
 					}
+				},
+				"lris": {
+					"landuse": {
+						"layers": [104400]
+					}
 				}
 			}
 		},
 		"general": {
-			"drop_offshore_lidar": true,
-			"lidar_classifications_to_keep": [1, 2, 4, 9],
+
+			"lidar_classifications_to_keep": [1, 2, 4, 6, 9],
 			"interpolation": {"no_data": "linear"}
 		},
-		"roughness": {
+		"roughness": { "default_values": {"rivers": 0.006, "lakes": null, "land": 0.1, "ocean": 0.004},
 			"roads": {"source": "osm", "ignore": ["pedestrian", "footway", "footpath", "track", "path", "cycleway"],
-                      "widths": {"default": 10, "residential": 20, "tertiary": 20, "secondary": 20, "motorway": 20}}
+                      "widths": {"default": 10, "residential": 20, "tertiary": 20, "secondary": 20, "motorway": 20}},
+			"landuse": {"source": "lris", "landcover": "Name_2018",
+						"classes_to_roughness": {"Transport Infrastructure": 0.001, "Exotic Forest": 0.6, "Low Producing Grassland": 0.09, "River": 0.006, "Built-up Area (settlement)": 0.05,
+									"Herbaceous Freshwater Vegetation": 0.082, "Indigenous Forest": 0.85, "Broadleaved Indigenous Hardwoods": 0.85, "Lake or Pond": 0.1,
+									"Manuka and/or Kanuka": 0.6, "Gorse and/or Broom": 0.4, "High Producing Exotic Grassland": 0.09, "Deciduous Hardwoods": 0.85, "Sand or Gravel": 0.008,
+									"Mixed Exotic Shrubland": 0.6, "Surface Mine or Dump": 0.004, "Orchard, Vineyard or Other Perennial Crop": 0.09, "Forest - Harvested": 0.004,
+									"Gravel or Rock": 0.0016, "Fernland": 0.6, "Matagouri or Grey Scrub": 0.6, "Urban Parkland/Open Space": 0.002}
+            }
 		}
 	}
 }