split out the LiDAR rasterisation from the coarse DEM regridding. Tes…

…t then add a processor step for patching in a raster.

src/geofabrics/dem.py

@@ -1859,6 +1859,59 @@ def _create_data_set(
 
         return dem
 
+class PatchDem(LidarBase):
+    """A class to manage the addition of a DEM to the foreground or background
+    of a preexisting DEM.
+
+    Parameters
+    ----------
+
+    patch_on_top
+        If True only patch the DEM values on top of the initial DEM. If False
+        patch only where values are NaN.
+    drop_patch_offshore
+        If True only keep patch values on land and the foreshore.
+    elevation_range
+        Optitionally specify a range of valid elevations. Any LiDAR points with
+        elevations outside this range will be filtered out.
+    initial_dem_path
+        The DEM to patch the other DEM on top / only where values are NaN.
+    buffer_cells - the number of empty cells to keep around LiDAR cells for
+        interpolation after the coarse DEM added to ensure a smooth boundary.
+    chunk_size
+        The chunk size in pixels for parallel/staged processing
+    """
+
+    def __init__(
+        self,
+        catchment_geometry: geometry.CatchmentGeometry,
+        patch_on_top: bool,
+        drop_patch_offshore: bool,
+        buffer_cells: int,
+        initial_dem_path: pathlib.Path | str,
+        elevation_range: list | None = None,
+        chunk_size: int | None = None
+    ):
+        """Setup base DEM to add future tiles too"""
+
+        super(RawDem, self).__init__(
+            catchment_geometry=catchment_geometry,
+            chunk_size=chunk_size,
+            elevation_range=elevation_range,
+        )
+        self.logger = logging.getLogger(f"{__name__}.{self.__class__.__name__}")
+
+        self.patch_on_top = patch_on_top
+        self.buffer_cells = buffer_cells
+        # Read in the DEM raster
+        initial_dem = rioxarray.rioxarray.open_rasterio(
+            pathlib.Path(initial_dem_path), masked=True, parse_coordinates=True, chunks=True
+        ).squeeze(
+            "band", drop=True
+        )  # remove band coordinate added by rasterio.open()
+        self._write_netcdf_conventions_in_place(initial_dem, catchment_geometry.crs)
+        self._dem = initial_dem
+
     def add_coarse_dem(self, coarse_dem_path: pathlib.Path, area_threshold: float):
         """Check if gaps in DEM on land, if so iterate through coarse DEMs
         adding missing detail.
@@ -2037,7 +2090,6 @@ def no_values_mask(self):
 
         return no_values_mask
 
-
 class RoughnessDem(LidarBase):
     """A class to add a roughness (zo) layer to a hydrologically conditioned DEM.
 

src/geofabrics/processor.py

@@ -984,6 +984,15 @@ def run(self):
                     key="coarse_dems", data_type="raster"
                 )
                 self.logger.info(f"Incorporating coarse DEMs: {coarse_dem_paths}")
+                self.raw_dem = dem.PatchDem(
+                    catchment_geometry=self.catchment_geometry,
+                    patch_on_top=False,
+                    drop_patch_offshore=True,
+                    initial_dem_path=cached_file,
+                    elevation_range=self.get_instruction_general("elevation_range"),
+                    chunk_size=self.get_processing_instructions("chunk_size"),
+                    buffer_cells=self.get_instruction_general("lidar_buffer"),
+                )
 
                 # Add coarse DEMs if there are any and if area
                 for coarse_dem_path in coarse_dem_paths: