Skip to content

Commit

Permalink
Fixed lake behaviour
Browse files Browse the repository at this point in the history
  • Loading branch information
@rosepearson
rosepearson committed Nov 20, 2024
1 parent dcc694f commit 8d756f8
Show file tree
Hide file tree
Showing 2 changed files with 116 additions and 91 deletions.
173 changes: 89 additions & 84 deletions src/geofabrics/dem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1253,7 +1253,8 @@ def add_points_within_polygon_nearest_chunked(
method: str,
cache_path: pathlib.Path,
label: str,
k_nearest_neighbours: int = 100,
k_nearest_neighbours: int,
include_edges: bool = True,
) -> xarray.Dataset:
"""Performs interpolation from estimated bathymetry points within a polygon
using the specified interpolation approach after filtering the points based
Expand All @@ -1267,7 +1268,7 @@ def add_points_within_polygon_nearest_chunked(
"method": method,
"crs": crs,
"k_nearest_neighbours": k_nearest_neighbours,
"use_edge": True,
"use_edge": include_edges,
"strict": False,
}
if method == "rbf":
Expand All @@ -1290,92 +1291,94 @@ def add_points_within_polygon_nearest_chunked(
pdal_pipeline.execute()

# Tempoarily save the adjacent points from the DEM
edge_dem = self._dem.rio.clip(
region_to_rasterise.dissolve().buffer(self.catchment_geometry.resolution),
drop=True,
)
self._write_netcdf_conventions_in_place(edge_dem, self.catchment_geometry.crs)
edge_dem["z"] = edge_dem.z.rio.interpolate_na(method="nearest")
edge_dem = self._dem.rio.clip(
region_to_rasterise.dissolve().buffer(self.catchment_geometry.resolution),
drop=True,
)
edge_dem = edge_dem.rio.clip(
region_to_rasterise.dissolve().geometry,
invert=True,
drop=True,
)

# Save provided points
grid_x, grid_y = numpy.meshgrid(edge_dem.x, edge_dem.y)
flat_x = grid_x.flatten()
flat_y = grid_y.flatten()
flat_z = edge_dem.z.values.flatten()
mask_z = ~numpy.isnan(flat_z)

# Interpolate the estimated bank heights around the polygon if they exist
if elevations.bank_heights_exist():
# Get the estimated bank heights and define a mask where nan
bank_points = elevations.bank_height_points()
bank_nan_mask = numpy.logical_not(numpy.isnan(bank_points["Z"]))
# Interpolate from the estimated bank heights
xy_out = numpy.concatenate(
[[flat_x[mask_z]], [flat_y[mask_z]]], axis=0
).transpose()
options = {
"radius": elevations.points["width"].max(),
"raster_type": geometry.RASTER_TYPE,
"method": "linear",
"strict": False,
}
estimated_edge_z = elevation_from_points(
point_cloud=bank_points[bank_nan_mask],
xy_out=xy_out,
options=options,
if include_edges:
edge_dem = self._dem.rio.clip(
region_to_rasterise.dissolve().buffer(self.catchment_geometry.resolution),
drop=True,
)
self._write_netcdf_conventions_in_place(edge_dem, self.catchment_geometry.crs)
edge_dem["z"] = edge_dem.z.rio.interpolate_na(method="nearest")
edge_dem = self._dem.rio.clip(
region_to_rasterise.dissolve().buffer(self.catchment_geometry.resolution),
drop=True,
)
edge_dem = edge_dem.rio.clip(
region_to_rasterise.dissolve().geometry,
invert=True,
drop=True,
)

# Use the estimated bank heights where lower than the DEM edge values
mask_z_edge = mask_z.copy()
mask_z_edge[:] = False
mask_z_edge[mask_z] = flat_z[mask_z] > estimated_edge_z
flat_z[mask_z_edge] = estimated_edge_z[flat_z[mask_z] > estimated_edge_z]
# Save provided points
grid_x, grid_y = numpy.meshgrid(edge_dem.x, edge_dem.y)
flat_x = grid_x.flatten()
flat_y = grid_y.flatten()
flat_z = edge_dem.z.values.flatten()
mask_z = ~numpy.isnan(flat_z)

# Interpolate the estimated bank heights around the polygon if they exist
if elevations.bank_heights_exist():
# Get the estimated bank heights and define a mask where nan
bank_points = elevations.bank_height_points()
bank_nan_mask = numpy.logical_not(numpy.isnan(bank_points["Z"]))
# Interpolate from the estimated bank heights
xy_out = numpy.concatenate(
[[flat_x[mask_z]], [flat_y[mask_z]]], axis=0
).transpose()
options = {
"radius": elevations.points["width"].max(),
"raster_type": geometry.RASTER_TYPE,
"method": "linear",
"strict": False,
}
estimated_edge_z = elevation_from_points(
point_cloud=bank_points[bank_nan_mask],
xy_out=xy_out,
options=options,
)

# Use the flat_x/y/z to define edge points and heights
edge_points = numpy.empty(
[mask_z.sum().sum()],
dtype=[
("X", geometry.RASTER_TYPE),
("Y", geometry.RASTER_TYPE),
("Z", geometry.RASTER_TYPE),
],
)
edge_points["X"] = flat_x[mask_z]
edge_points["Y"] = flat_y[mask_z]
edge_points["Z"] = flat_z[mask_z]
# Use the estimated bank heights where lower than the DEM edge values
mask_z_edge = mask_z.copy()
mask_z_edge[:] = False
mask_z_edge[mask_z] = flat_z[mask_z] > estimated_edge_z
flat_z[mask_z_edge] = estimated_edge_z[flat_z[mask_z] > estimated_edge_z]

edge_file = cache_path / f"{label}_edge_points.laz"
pdal_pipeline_instructions = [
{
"type": "writers.las",
"a_srs": f"EPSG:" f"{crs['horizontal']}+" f"{crs['vertical']}",
"filename": str(edge_file),
"compression": "laszip",
}
]
pdal_pipeline = pdal.Pipeline(
json.dumps(pdal_pipeline_instructions), [edge_points]
)
pdal_pipeline.execute()
# Use the flat_x/y/z to define edge points and heights
edge_points = numpy.empty(
[mask_z.sum().sum()],
dtype=[
("X", geometry.RASTER_TYPE),
("Y", geometry.RASTER_TYPE),
("Z", geometry.RASTER_TYPE),
],
)
edge_points["X"] = flat_x[mask_z]
edge_points["Y"] = flat_y[mask_z]
edge_points["Z"] = flat_z[mask_z]

edge_file = cache_path / f"{label}_edge_points.laz"
pdal_pipeline_instructions = [
{
"type": "writers.las",
"a_srs": f"EPSG:" f"{crs['horizontal']}+" f"{crs['vertical']}",
"filename": str(edge_file),
"compression": "laszip",
}
]
pdal_pipeline = pdal.Pipeline(
json.dumps(pdal_pipeline_instructions), [edge_points]
)
pdal_pipeline.execute()

if (
len(points) < raster_options["k_nearest_neighbours"]
or len(edge_points) < raster_options["k_nearest_neighbours"]
or (include_edges and len(edge_points) < raster_options["k_nearest_neighbours"])
):
k_nearest_neighbours = min(len(points), len(edge_points)) if include_edges else len(points)
logging.info(
f"Fewer points or edge points than the default expected {raster_options['k_nearest_neighbours']}. "
f"Updating k_nearest_neighbours to {min(len(points), len(edge_points))}."
f"Updating k_nearest_neighbours to {k_nearest_neighbours}."
)
raster_options["k_nearest_neighbours"] = min(len(points), len(edge_points))
raster_options["k_nearest_neighbours"] = k_nearest_neighbours
if raster_options["k_nearest_neighbours"] < 3:
logging.warning(
f"Not enough points or edge points to meaningfully include {raster_options['k_nearest_neighbours']}. "
Expand Down Expand Up @@ -1409,13 +1412,15 @@ def add_points_within_polygon_nearest_chunked(
chunk_region_to_tile=None,
crs=raster_options["crs"],
)

edge_points = delayed_load_tiles_in_chunk(
lidar_files=[edge_file],
source_crs=raster_options["crs"],
chunk_region_to_tile=None,
crs=raster_options["crs"],
)
if include_edges:
edge_points = delayed_load_tiles_in_chunk(
lidar_files=[edge_file],
source_crs=raster_options["crs"],
chunk_region_to_tile=None,
crs=raster_options["crs"],
)
else:
edge_points = None

# Rasterise tiles
delayed_chunked_x.append(
Expand Down
34 changes: 27 additions & 7 deletions src/geofabrics/processor.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,10 +304,18 @@ def get_instruction_general(self, key: str, subkey: str = None):
"ocean": None,
},
"ignore_clipping": False,
"ocean": {
"nearest_k_for_interpolation": 40,
"use_edge": False,
"is_depth": False,
"nearest_k_for_interpolation": {
"ocean": 40,
"lake": 500,
"rivers": 100,
},
"use_edge": {
"ocean": False,
"lake": False,
},
"is_depth": {
"ocean": False,
"lake": False,
},
"filter_waterways_by_osm_ids": [],
"compression": 1,
Expand Down Expand Up @@ -1178,18 +1186,18 @@ def add_hydrological_features(
key="z_labels", subkey="ocean"
),
is_depth=self.get_instruction_general(
key="ocean", subkey="is_depth"
key="is_depth", subkey="ocean"
),
)
# Interpolate
hydrologic_dem.interpolate_ocean_chunked(
ocean_points=ocean_points,
cache_path=temp_folder,
use_edge=self.get_instruction_general(
key="ocean", subkey="use_edge"
key="use_edge", subkey="ocean"
),
k_nearest_neighbours=self.get_instruction_general(
key="ocean", subkey="nearest_k_for_interpolation"
key="nearest_k_for_interpolation", subkey="ocean"
),
buffer=self.get_instruction_general(key="lidar_buffer"),
method=self.get_instruction_general(
Expand Down Expand Up @@ -1295,6 +1303,9 @@ def add_hydrological_features(
polygon_files=[polygon],
catchment_geometry=self.catchment_geometry,
z_labels=z_labels[index],
is_depth=self.get_instruction_general(
key="is_depth", subkey="lakes"
),
)

if (
Expand All @@ -1316,6 +1327,12 @@ def add_hydrological_features(
),
cache_path=temp_folder,
label="lakes",
include_edges=self.get_instruction_general(
key="use_edge", subkey="lakes"
),
k_nearest_neighbours=self.get_instruction_general(
key="nearest_k_for_interpolation", subkey="lakes"
),
)
temp_file = temp_folder / f"dem_added_{index + 1}_lake.nc"
self.logger.info(
Expand Down Expand Up @@ -1377,6 +1394,9 @@ def add_hydrological_features(
),
cache_path=temp_folder,
label="rivers and fans",
k_nearest_neighbours=self.get_instruction_general(
key="nearest_k_for_interpolation", subkey="rivers"
),
)
temp_file = temp_folder / f"dem_added_{index + 1}_rivers.nc"
self.logger.info(
Expand Down

0 comments on commit 8d756f8

Please sign in to comment.