Skip to content

Commit

Permalink
Tidy-up to remove the now unused coarse DEM chunking code
Browse files Browse the repository at this point in the history
  • Loading branch information
rosepearson committed Nov 3, 2023
1 parent a9c791d commit 4037956
Showing 1 changed file with 0 additions and 157 deletions.
157 changes: 0 additions & 157 deletions src/geofabrics/dem.py
Original file line number Diff line number Diff line change
Expand Up @@ -1836,163 +1836,6 @@ def add_coarse_dems(
0,
)

"""if chunk_size is None:
self._add_coarse_dem_no_chunking(
coarse_dem_path=coarse_dem_path,
mask=no_value_mask,
)
else:
self._add_coarse_dem_chunked(
coarse_dem_path=coarse_dem_path,
chunk_size=chunk_size,
mask=no_value_mask,
radius=coarse_dem_resolution * numpy.sqrt(2),
)"""

def _add_coarse_dem_no_chunking(
self,
coarse_dem_path: pathlib.Path,
mask: numpy.ndarray,
):
"""Fill gaps in dense DEM from areas with no LiDAR with the coarse DEM.
Perform linear interpolation.
Parameters
----------
coarse_dem - The coarse DEM to use for rasterising
mask - the pixel mask of where to set these values
"""

# Load in the coarse DEM
extents = {
"total": self.catchment_geometry.catchment,
"land": self.catchment_geometry.full_land,
"foreshore": self.catchment_geometry.foreshore,
}
coarse_dem = CoarseDem(
dem_file=coarse_dem_path,
extents=extents, # by reference, so "total" trimmed to coarse DEM bounds
set_foreshore=self.drop_offshore_lidar,
)

# create dictionary defining raster options
# Set search radius to the diagonal cell length to ensure corners covered
raster_options = {
"raster_type": geometry.RASTER_TYPE,
"radius": coarse_dem.resolution * numpy.sqrt(2),
"method": "linear",
}
# Get the grid locations overwhich to perform averaging
grid_x, grid_y = numpy.meshgrid(self._dem.x, self._dem.y)

# Mask to only rasterise where there aren't already LiDAR derived values
xy_out = numpy.empty((mask.values.sum(), 2))
xy_out[:, 0] = grid_x[mask]
xy_out[:, 1] = grid_y[mask]

# Perform specified averaging from the coarse DEM where there is no data
z_flat = elevation_from_points(
point_cloud=coarse_dem.points,
xy_out=xy_out,
options=raster_options,
)
# Update the DEM
self._dem.z.data[mask] = z_flat
# Update the data source layer
self._dem["data_source"] = self._dem.data_source.where(
~(mask & self._dem.z.notnull().values),
self.SOURCE_CLASSIFICATION["coarse DEM"],
)

def _add_coarse_dem_chunked(
self,
coarse_dem_path: pathlib.Path,
chunk_size: int,
mask: numpy.ndarray,
radius: float,
):
"""Fill gaps in dense DEM from areas with no LiDAR with the coarse DEM.
Perform linear interpolation.
Parameters
----------
coarse_dem - The coarse DEM to use for rasterising
extents - a dictionary of area to rasterise using the coarse DEM
chunk_size - the size of each chunk
"""

# Define raster options
raster_options = {
"raster_type": geometry.RASTER_TYPE,
"elevation_range": self.elevation_range,
"radius": radius,
"method": "linear",
}
extents = {
"total": self.catchment_geometry.catchment,
"land": self.catchment_geometry.full_land,
"foreshore": self.catchment_geometry.foreshore,
}

# get chunking information
chunked_dim_x, chunked_dim_y = self._chunks_from_dem(chunk_size, self._dem)
elevations = {}

# cycle through index chunks - and collect in a delayed array
delayed_chunked_matrix = []
for i, dim_y in enumerate(chunked_dim_y):
delayed_chunked_x = []
for j, dim_x in enumerate(chunked_dim_x):
logging.info(f"\tCoarse chunk {[i, j]}")

# Define the region of the chunk to rasterise
chunk_region_to_tile = self._define_chunk_region(
region_to_rasterise=extents["total"],
dim_x=dim_x,
dim_y=dim_y,
radius=raster_options["radius"],
)
chunk_extents = { # reset to unbuffered land and foreshore
"total": chunk_region_to_tile,
"land": self.catchment_geometry.full_land,
"foreshore": self.catchment_geometry.foreshore,
}

# Send through the clipped Coarse DEM points
coarse_points = delayed_chunk_coarse_dem(
dem_file=coarse_dem_path,
extents=chunk_extents,
set_foreshore=self.drop_offshore_lidar,
)

# Use the coase DEM to sample any missing values
delayed_chunked_x.append(
dask.array.from_delayed(
delayed_elevation_over_chunk(
dim_x=dim_x,
dim_y=dim_y,
tile_points=coarse_points,
options=raster_options,
),
shape=(len(dim_y), len(dim_x)),
dtype=geometry.RASTER_TYPE,
)
)
delayed_chunked_matrix.append(delayed_chunked_x)

# Combine chunks into a array and replace missing values in dataset
elevations = dask.array.block(delayed_chunked_matrix)
self._dem["z"] = self._dem.z.where(~mask, elevations)
# Update the data source layer
self._dem["data_source"] = self._dem.data_source.where(
~(mask & self._dem.z.notnull()),
self.SOURCE_CLASSIFICATION["coarse DEM"],
)


class RoughnessDem(LidarBase):
"""A class to add a roughness (zo) layer to a hydrologically conditioned DEM.
Expand Down

0 comments on commit 4037956

Please sign in to comment.