r.hydro.flatten: new addon to derive elevation of water bodies for hy…

…dro-flattening (#880)

src/raster/r.hydro.flatten/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.hydro.flatten
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

src/raster/r.hydro.flatten/r.hydro.flatten.html

@@ -0,0 +1,82 @@
+<h2>DESCRIPTION</h2>
+The tool derives single elevation value for water bodies
+based on lidar data. These values are used for hydro-flattening a digital elevation model.
+The <b>input</b> raster is expected to represent ground surface created by binning lidar data
+(e.g., using <em><a href="r.in.pdal.html">r.in.pdal</a></em>) with averaged ground elevation.
+Small gaps in the input are expected. Large gaps are interpreted as water bodies.
+The minimum size of a water body can be set with <b>min_size</b> option in map units.
+
+<p>
+The output <b>water_elevation</b> is a raster map of water bodies where each water body
+has a single value representing the water level elevation derived from the lidar data
+at the edge of a water body.
+Since the elevation varies along the edge, option <b>percentile</b> is used to determine
+a single value. The variation along the edge can be examined with the <b>water_elevation_stddev</b>
+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.
+
+<p>
+To keep the intermediate results for inspection, use flag <b>-k</b>.
+
+<h2>NOTES</h2>
+While this tool was designed for water bodies, it can be used for other purposes,
+e.g., for filling a gap in digital elevation models caused by excluding buildings.
+<p>
+This tool does not interpolate gaps in data, rather it derives a single value
+for each gap. The result can be used to fill gaps and the tool can be run on large areas.
+For actual gap interpolation, which is typically more computationally intensive,
+see <em><a href="https://grass.osgeo.org/grass-stable/manuals/r.fillnuls.html">r.fillnulls</a></em>.
+
+<h2>EXAMPLE</h2>
+We will download a lidar tile with <em><a href="r.in.usgs.html">r.in.usgs</a></em> addon,
+use <em><a href="https://grass.osgeo.org/grass-stable/manuals/r.in.pdal.html">r.in.pdal</a></em>
+to bin the elevation points at 1 meter resolution, and derive elevation levels for lakes
+with minimum size of 4000 m^2.
+
+<div class="code"><pre>
+# select study area and resolution
+g.region n=213300 s=211900 w=653900 e=655300 res=1
+# download lidar tile into /tmp
+r.in.usgs product=lidar output_directory=/tmp title_filter=Phase2 -d
+# bin point elevation using ground and road points with reprojection
+r.in.pdal input=/tmp/USGS_LPC_NC_Phase2_2014_LA_37_20164902_.laz output=ground -w class_filter=2,13
+# 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
+</pre></div>
+
+<div align="center" style="margin: 10px">
+<!--
+r.colors map=ground_m,water_elevation color=elevation
+r.to.vect -t input=water_elevation output=water_elevation type=area
+d.rast map=ground_m
+d.rast map=water_elevation
+d.rast map=water_elevation_stddev
+d.vect map=water_elevation color=77:77:77:255 fill_color=none
+d.legend -t -b raster=water_elevation_stddev title="Stddev [m]" digits=2 border_color=none
+-->
+  <img src="r_hydro_flatten_input.png" alt="input for r.hydro.flatten" border="0">
+  <img src="r_hydro_flatten_output_elevation.png" alt="output elevation from r.hydro.flatten" border="0">
+  <img src="r_hydro_flatten_output_stddev.png" alt="output stddev from r.hydro.flatten" border="0">
+  <br>
+  <i>Figure: Input binned elevation representing ground with gaps (left),
+     input overlayed with elevation values estimated for gaps and highlighted with an outline (middle),
+     input overlayed with standard deviation of the elevation along the edge of the gaps (right).</i>
+</div>
+
+<h2>REFERENCE</h2>
+Method based on workflow <a href="https://www.youtube.com/watch?v=p9KCfufNYgE">presented</a>
+at NC GIS Conference 2021 by Doug Newcomb.
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="https://grass.osgeo.org/grass-stable/manuals/r.in.pdal.html">r.in.pdal</a>,
+<a href="r.in.usgs.html">r.in.usgs</a>,
+<em><a href="https://grass.osgeo.org/grass-stable/manuals/r.fillnulls.html">r.fillnulls</a></em>
+<em>
+
+<h2>AUTHOR</h2>
+
+Anna Petrasova, <a href="http://geospatial.ncsu.edu/geoforall/">NCSU GeoForAll Lab</a>

src/raster/r.hydro.flatten/r.hydro.flatten.py

@@ -0,0 +1,203 @@
+#!/usr/bin/env python
+#
+#########################################################################
+#
+# MODULE:     r.hydro.flatten
+#
+# AUTHOR(S):  Anna Petrasova <kratochanna gmail com>
+#
+# PURPOSE:    Derive elevation of water bodies for hydro-flattening
+#
+# COPYRIGHT:  (C) 2023 by Anna Petrasova, and the GRASS Development Team
+#
+#             This program is free software under the GNU General Public
+#             License (>=v2). Read the COPYING file that comes with GRASS
+#             for details.
+#
+#########################################################################
+
+# %module
+# % description: Derive elevation of water bodies for hydro-flattening
+# % keyword: raster
+# % keyword: elevation
+# % keyword: hydrology
+# % keyword: lidar
+# % keyword: LIDAR
+# %end
+# %option G_OPT_R_INPUT
+# % key: input
+# % description: Raster map of binned lidar point elevation
+# %end
+# %option G_OPT_R_OUTPUT
+# % key: water_elevation
+# % description: Represents single elevation value for each water body
+# % label: Raster map of derived water elevation
+# %end
+# %option G_OPT_R_OUTPUT
+# % key: water_elevation_stddev
+# % description: Raster map of derived water elevation standard deviation
+# %end
+# %option
+# % key: percentile
+# % type: double
+# % required: yes
+# % description: Percentile of elevation to determine water level
+# % answer: 5
+# %end
+# %option
+# % key: min_size
+# % type: integer
+# % required: no
+# % description: Minimum size of areas in map units
+# %end
+# %flag
+# % key: k
+# % description: Keep intermediate results
+# %end
+
+import os
+import sys
+import atexit
+from math import sqrt
+
+import grass.script as gs
+
+RAST_REMOVE = []
+
+
+def cleanup():
+    if RAST_REMOVE:
+        gs.run_command("g.remove", flags="f", type="raster", name=RAST_REMOVE)
+
+
+def get_tmp_name(basename):
+    name = gs.append_node_pid(basename)
+    RAST_REMOVE.append(name)
+    return name
+
+
+def main():
+    options, flags = gs.parser()
+    keep = flags["k"]
+    if keep:
+
+        def get_name(basename):
+            return f"intermediate_{basename}"
+
+    else:
+
+        def get_name(basename):
+            name = gs.append_node_pid(basename)
+            RAST_REMOVE.append(name)
+            return name
+
+    ground = options["input"]
+    size_threshold = options["min_size"]
+    if size_threshold:
+        size_threshold = int(size_threshold)
+    else:
+        size_threshold = None
+    # r.fill.stats settings
+    filling_distance = 3
+    filling_cells = 6
+    region = gs.region()
+    radius = sqrt(region["nsres"] * region["ewres"])
+    # distance range to get a 1-cell wide edge in "clump1"
+    max_radius = radius * 4.01
+    min_radius = radius * 3.01
+    tmp_rfillstats = get_name("rfillstats")
+    gs.run_command(
+        "r.fill.stats",
+        flags="k",
+        input=ground,
+        output=tmp_rfillstats,
+        distance=filling_distance,
+        cells=filling_cells,
+    )
+    tmp_water_buffer = get_name("water_buffer")
+    gs.run_command(
+        "r.grow.distance",
+        flags="n",
+        input=tmp_rfillstats,
+        distance=tmp_water_buffer,
+        metric="squared",
+        max_distance=max_radius * max_radius,
+    )
+    tmp_clump1 = get_name("clump1")
+    gs.mapcalc(
+        f"{tmp_clump1} = if ({tmp_water_buffer} < ({min_radius} * {min_radius}), null(), 1)"
+    )
+    tmp_clump2 = get_name("clump2")
+    gs.run_command("r.clump", flags="d", input=tmp_clump1, output=tmp_clump2)
+    tmp_water_elevation = get_name("water_elevation")
+    gs.run_command(
+        "r.stats.quantile",
+        base=tmp_clump2,
+        cover=tmp_rfillstats,
+        percentiles=options["percentile"],
+        output=tmp_water_elevation,
+    )
+    tmp_water_stddev = get_name("water_stddev")
+    gs.run_command(
+        "r.stats.zonal",
+        base=tmp_clump2,
+        cover=tmp_rfillstats,
+        method="stddev",
+        output=tmp_water_stddev,
+    )
+    tmp_water_elevation_dist = get_name("water_elevation_dist")
+    gs.run_command(
+        "r.grow.distance", input=tmp_water_elevation, value=tmp_water_elevation_dist
+    )
+    tmp_water_stddev_dist = get_name("water_stddev_dist")
+    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_water_buffer} < {min_radius} * {min_radius}, "
+        f"{tmp_water_elevation_dist}, null())"
+    )
+    tmp_water_stddev_dist_res = get_name("water_stddev_dist_res")
+    gs.mapcalc(
+        f"{tmp_water_stddev_dist_res} = if ({tmp_water_buffer} < {min_radius} * {min_radius}, "
+        f"{tmp_water_stddev_dist}, 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,
+            output=tmp_reclass,
+            rules="-",
+            stdin="* = 1",
+        )
+        tmp_clump_reclass = get_name("clump_reclass")
+        gs.run_command("r.clump", input=tmp_reclass, output=tmp_clump_reclass)
+        tmp_size = get_name("size")
+        gs.run_command(
+            "r.stats.zonal",
+            base=tmp_clump_reclass,
+            cover=tmp_reclass,
+            method="sum",
+            output=tmp_size,
+        )
+        gs.mapcalc(
+            f"{options['water_elevation']} = if ({tmp_size} > {size_threshold}, {tmp_water_elevation_dist_res}, null())"
+        )
+        gs.mapcalc(
+            f"{options['water_elevation_stddev']} = if ({tmp_size} > {size_threshold}, {tmp_water_stddev_dist_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.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 __name__ == "__main__":
+    atexit.register(cleanup)
+    sys.exit(main())