Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

r.hydro.flatten: add breaklines and filled_elevation output parameter #1281

Merged
merged 2 commits into from
Jan 21, 2025
Merged

r.hydro.flatten: add breaklines and filled_elevation output parameter #1281

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
62d6a8a
r.hydro.flatten: add breaklines and filled_elevation output parameter
petrasovaa Jan 15, 2025
a6cee31
make filled_elevation output optional
petrasovaa Jan 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion src/raster/r.hydro.flatten/r.hydro.flatten.html
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,14 @@ <h2>DESCRIPTION</h2>
output representing the standard deviation of the lidar elevation values
along the water body's edge. Higher deviation suggests problematic areas
that need to be further inspected.
The optional output <b>filled_elevation</b> is a raster map of the input ground surface
filled with the computed <b>water_elevation</b> raster map.

<p>
The <b>breaklines</b> parameter is an optional input that specifies a vector map of lines
that represent e.g., a break between an impoundment and downstream river, allowing
correct elevation computation.
</p>

<p>
To keep the intermediate results for inspection, use flag <b>-k</b>.
Expand Down Expand Up @@ -44,7 +52,7 @@ <h2>EXAMPLE</h2>
# convert elevation from feet to meters
r.mapcalc "ground_m = ground * 0.304800609601219"
# derive elevation of water bodies and standard deviation
r.hydro.flatten input=ground_m water_elevation=water_elevation water_elevation_stddev=water_elevation_stddev percentile=10 misize=4000
r.hydro.flatten input=ground_m water_elevation=water_elevation water_elevation_stddev=water_elevation_stddev filled_elevation=filled percentile=10 misize=4000
</pre></div>

<div align="center" style="margin: 10px">
Expand Down
98 changes: 81 additions & 17 deletions src/raster/r.hydro.flatten/r.hydro.flatten.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,11 @@
# % key: input
# % description: Raster map of binned lidar point elevation
# %end
# %option G_OPT_V_INPUT
# % key: breaklines
# % description: Vector map of breaklines
# % required: no
# %end
# %option G_OPT_R_OUTPUT
# % key: water_elevation
# % description: Represents single elevation value for each water body
Expand All @@ -37,6 +42,11 @@
# % key: water_elevation_stddev
# % description: Raster map of derived water elevation standard deviation
# %end
# %option G_OPT_R_OUTPUT
# % key: filled_elevation
# % required: no
# % description: Raster map representing filled digital elevation model
# %end
# %option
# % key: percentile
# % type: double
Expand Down Expand Up @@ -92,6 +102,7 @@ def get_name(basename):
return name

ground = options["input"]
breaklines = options["breaklines"]
size_threshold = options["min_size"]
if size_threshold:
size_threshold = int(size_threshold)
Expand All @@ -107,6 +118,17 @@ def get_name(basename):
region_m = gs.parse_command("g.region", flags="gm")
resolution_m = (float(region_m["nsres"]) + float(region_m["ewres"])) / 2
tmp_rfillstats = get_name("rfillstats")

if breaklines:
tmp_breaklines = get_name("breaklines")
gs.run_command(
"v.to.rast",
input=breaklines,
output=tmp_breaklines,
use="val",
flags="d",
value=1000,
)
gs.run_command(
"r.fill.stats",
flags="k",
Expand All @@ -125,6 +147,15 @@ def get_name(basename):
distances=[x * resolution_m for x in range(1, buffer_last_strip)],
units="meters",
)
if breaklines:
tmp_buffer_with_breaklines = get_name("buffer_with_breaklines")
gs.run_command(
"r.patch",
input=[tmp_breaklines, tmp_buffer],
output=tmp_buffer_with_breaklines,
)
tmp_buffer = tmp_buffer_with_breaklines

tmp_reclass_for_clump = get_name("reclass_for_clump")
gs.write_command(
"r.reclass",
Expand Down Expand Up @@ -155,30 +186,53 @@ def get_name(basename):
method="stddev",
output=tmp_water_stddev,
)
tmp_water_elevation_dist = get_name("water_elevation_dist")
tmp_water_elevation_zonal = get_name("water_elevation_zonal")
gs.run_command(
"r.grow.distance", input=tmp_water_elevation, value=tmp_water_elevation_dist
"r.stats.zonal",
base=tmp_clump,
cover=tmp_water_elevation,
method="average",
output=tmp_water_elevation_zonal,
)
tmp_water_stddev_dist = get_name("water_stddev_dist")
tmp_water_elevation_stddev_zonal = get_name("water_elevation_stddev_zonal")
gs.run_command(
"r.grow.distance", input=tmp_water_stddev, value=tmp_water_stddev_dist
)
tmp_water_elevation_dist_res = get_name("water_elevation_dist_res")
gs.mapcalc(
f"{tmp_water_elevation_dist_res} = if ({tmp_buffer} < {buffer_last_strip}, "
f"{tmp_water_elevation_dist}, null())"
"r.stats.zonal",
base=tmp_clump,
cover=tmp_water_stddev,
method="average",
output=tmp_water_elevation_stddev_zonal,
)
tmp_water_stddev_dist_res = get_name("water_stddev_dist_res")
tmp_water_elevation_zonal_res = get_name("water_elevation_zonal_res")
if breaklines:
water_elevation_zonal_res_breaklines = get_name(
"water_elevation_zonal_res_breaklines"
)
gs.mapcalc(
f"{water_elevation_zonal_res_breaklines} = if (isnull({tmp_strip}), {tmp_water_elevation_zonal}, null())"
)
# heal the breakline holes
gs.run_command(
"r.neighbors",
input=water_elevation_zonal_res_breaklines,
selection=tmp_breaklines,
output=tmp_water_elevation_zonal_res,
size=5,
)
else:
gs.mapcalc(
f"{tmp_water_elevation_zonal_res} = if (isnull({tmp_strip}), {tmp_water_elevation_zonal}, null())"
)

tmp_water_elevation_stddev_zonal_res = get_name("water_elevation_stddev_zonal_res")
gs.mapcalc(
f"{tmp_water_stddev_dist_res} = if ({tmp_buffer} < {buffer_last_strip}, "
f"{tmp_water_stddev_dist}, null())"
f"{tmp_water_elevation_stddev_zonal_res} = if (isnull({tmp_strip}), {tmp_water_elevation_stddev_zonal}, null())"
)
if size_threshold:
size_threshold /= region["nsres"] * region["ewres"]
tmp_reclass = get_name("reclass")
gs.write_command(
"r.reclass",
input=tmp_water_elevation_dist_res,
input=tmp_water_elevation_zonal_res,
output=tmp_reclass,
rules="-",
stdin="* = 1",
Expand All @@ -194,18 +248,28 @@ def get_name(basename):
output=tmp_size,
)
gs.mapcalc(
f"{options['water_elevation']} = if ({tmp_size} > {size_threshold}, {tmp_water_elevation_dist_res}, null())"
f"{options['water_elevation']} = if ({tmp_size} > {size_threshold}, {tmp_water_elevation_zonal_res}, null())"
)
gs.mapcalc(
f"{options['water_elevation_stddev']} = if ({tmp_size} > {size_threshold}, {tmp_water_stddev_dist_res}, null())"
f"{options['water_elevation_stddev']} = if ({tmp_size} > {size_threshold}, {tmp_water_elevation_stddev_zonal_res}, null())"
)
else:
gs.mapcalc(f"{options['water_elevation']} = {tmp_water_elevation_dist_res}")
gs.mapcalc(f"{options['water_elevation_stddev']} = {tmp_water_stddev_dist_res}")
gs.mapcalc(f"{options['water_elevation']} = {tmp_water_elevation_zonal_res}")
gs.mapcalc(
f"{options['water_elevation_stddev']} = {tmp_water_elevation_stddev_zonal_res}"
)
gs.run_command("r.colors", map=options["water_elevation"], raster=ground)
gs.run_command("r.colors", map=options["water_elevation_stddev"], color="reds")
gs.raster_history(options["water_elevation"])
gs.raster_history(options["water_elevation_stddev"])
if options["filled_elevation"]:
gs.run_command(
"r.patch",
input=[tmp_rfillstats, options["water_elevation"]],
output=options["filled_elevation"],
)
gs.run_command("r.colors", map=options["filled_elevation"], raster=ground)
gs.raster_history(options["filled_elevation"])


if __name__ == "__main__":
Expand Down
Loading