Split HAG filter into three (#2927)

* Updates.

* Add filters.hag_dem

* make filters.hag_dem streamable, improve docs

* update copyright info

* Make raster a positional argument in HAGDEMFilter

* Break HAG into pieces.

* Fix validation error message for invalid band parameter

* Formatting changes.  Remove "static" in processOne()..

* Review changes.

* Update documentation with new graphics and proper refs

* Add back HAG filter for now.

* Fix spelling.

Co-authored-by: Julian Fell <hi@jtfell.com>
Co-authored-by: Bradley J Chambers <brad.chambers@gmail.com>

doc/stages/filters.hag.rst

@@ -3,18 +3,22 @@
 filters.hag
 ===============================================================================
 
-The **Height Above Ground (HAG) filter** takes as input a point cloud with
-``Classification`` set to 2 for ground points.  It creates a new dimension,
-``HeightAboveGround``, that contains the normalized height values.
-
 .. note::
 
-   We expect ground returns to have the classification value of 2 in keeping
-   with the `ASPRS Standard LIDAR Point Classes <http://www.asprs.org/a/society/committees/standards/LAS_1_4_r13.pdf>`_.
+    The HAG filter is deprecated and has been replaced by :ref:`filters.hag_nn`,
+    :ref:`filters.hag_dem` and :ref:`filters.hag_delaunay`.  Please use
+    the new filter
+    that matches your needs as this filter will be removed in a future release.
 
-Ground points may be generated by :ref:`filters.pmf` or
-:ref:`filters.smrf`, but you can use any method
-you choose, as long as the ground returns are marked.
+The **Height Above Ground (HAG) filter** takes as input a point cloud with
+``Classification`` set to 2 for ground points (see
+`ASPRS Standard LIDAR Point Classes
+<http://www.asprs.org/a/society/committees/standards/LAS_1_4_r13.pdf>`_).
+It creates a new dimension,
+``HeightAboveGround``, that contains the normalized height values.  Ground
+points may be generated by :ref:`filters.pmf` or :ref:`filters.smrf`,
+but you can use any method you choose, as long as the ground returns are
+marked with the classification value of 2.
 
 Normalized heights are a commonly used attribute of point cloud data. This can
 also be referred to as *height above ground* (HAG) or *above ground level* (AGL)

doc/stages/filters.hag_delaunay.rst

@@ -0,0 +1,90 @@
+.. _filters.hag_delaunay:
+
+filters.hag_delaunay
+===============================================================================
+
+The **Height Above Ground Delaunay filter** takes as input a point cloud with
+``Classification`` set to 2 for ground points.  It creates a new dimension,
+``HeightAboveGround``, that contains the normalized height values.
+
+.. note::
+
+   We expect ground returns to have the classification value of 2 in keeping
+   with the `ASPRS Standard LIDAR Point Classes
+   <http://www.asprs.org/a/society/committees/standards/LAS_1_4_r13.pdf>`_.
+
+Ground points may be generated by :ref:`filters.pmf` or :ref:`filters.smrf`,
+but you can use any method you choose, as long as the ground returns are
+marked.
+
+Normalized heights are a commonly used attribute of point cloud data. This can
+also be referred to as *height above ground* (HAG) or *above ground level*
+(AGL) heights. In the end, it is simply a measure of a point's relative height
+as opposed to its raw elevation value.
+
+The filter creates a delaunay triangulation of the `count`_ ground points
+closest to the non-ground point in question.  If the non-ground point is within
+the trianulated area, the assigned ``HeightAboveGround`` is the difference
+between its ``Z`` value and a ground height interpolated from the three
+vertices of the containing triangle.  If the non-ground point is outside of the
+triangulated area, its ``HeightAboveGround`` is calculated as the difference
+between its ``Z`` value and the ``Z`` value of the nearest ground point.
+
+Choosing a value for `count`_ is difficult, as placing the non-ground point in
+the triangulated area depends on the layout of the nearby points.  If, for
+example, all the ground points near a non-ground point lay on one side of that
+non-ground point, finding a containing triangle will fail.
+
+.. embed::
+
+Example #1
+----------
+
+Using the autzen dataset (here shown colored by elevation), which already has
+points classified as ground
+
+.. image:: ./images/autzen-elevation.png
+   :height: 400px
+
+we execute the following pipeline
+
+.. code-block:: json
+
+  [
+      "autzen.laz",
+      {
+          "type":"filters.hag_delaunay"
+      },
+      {
+          "type":"writers.laz",
+          "filename":"autzen_hag_delaunay.laz",
+          "extra_dims":"HeightAboveGround=float32"
+      }
+  ]
+
+which is equivalent to the ``pdal translate`` command
+
+::
+
+    $ pdal translate autzen.laz autzen_hag_delaunay.laz hag_delaunay \
+        --writers.las.extra_dims="HeightAboveGround=float32"
+
+In either case, the result, when colored by the normalized height instead of
+elevation is
+
+.. image:: ./images/autzen-hag-delaunay.png
+   :height: 400px
+
+Options
+-------------------------------------------------------------------------------
+
+_`count`
+    The number of ground neighbors to consider when determining the height
+    above ground for a non-ground point.  [Default: 10]
+
+allow_extrapolation
+    If false and a non-ground point lies outside of the bounding box of
+    all ground points, its ``HeightAboveGround`` is set to 0.  If true
+    and ``delaunay`` is set, the ``HeightAboveGround`` is set to the
+    difference between the heights of the non-ground point and nearest
+    ground point.  [Default: false]

doc/stages/filters.hag_dem.rst

@@ -0,0 +1,84 @@
+.. _filters.hag_dem:
+
+filters.hag_dem
+===============================================================================
+
+The **Height Above Ground (HAG) Digital Elevation Model (DEM) filter** loads
+a GDAL-readable raster image specifying the DEM. The ``Z`` value of each point
+in the input is compared against the value at the corresponding X,Y location
+in the DEM raster. It creates a new dimension, ``HeightAboveGround``, that
+contains the normalized height values.
+
+Normalized heights are a commonly used attribute of point cloud data. This can
+also be referred to as *height above ground* (HAG) or *above ground level* (AGL)
+heights. In the end, it is simply a measure of a point's relative height as
+opposed to its raw elevation value.
+
+.. embed::
+
+.. streamable::
+
+Example #1
+----------
+
+Using the autzen dataset (here shown colored by elevation)
+
+.. image:: ./images/autzen-elevation.png
+   :height: 400px
+
+we generate a DEM based on the points already classified as ground
+
+::
+  
+    $ pdal translate autzen.laz autzen_dem.tif range \
+        --filters.range.limits="Classification[2:2]" \
+        --writers.gdal.output_type="idw" \
+        --writers.gdal.resolution=6 \
+        --writers.gdal.window_size=24
+
+and execute the following pipeline
+
+.. code-block:: json
+
+  [
+      "autzen.laz",
+      {
+          "type":"filters.hag_dem",
+          "raster": "autzen_dem.tif"
+      },
+      {
+          "type":"writers.las",
+          "filename":"autzen_hag_dem.laz",
+          "extra_dims":"HeightAboveGround=float32"
+      }
+  ]
+
+which is equivalent to the ``pdal translate`` command
+
+::
+
+    $ pdal translate autzen.laz autzen_hag_dem.laz hag_dem \
+        --filters.hag_dem.raster=autzen_dem.tif \
+        --writers.las.extra_dims="HeightAboveGround=float32"
+
+In either case, the result, when colored by the normalized height instead of
+elevation is
+
+.. image:: ./images/autzen-hag-dem.png
+   :height: 400px
+
+Options
+-------------------------------------------------------------------------------
+
+_`raster`
+    GDAL-readable raster to use for DEM.
+
+band
+    GDAL Band number to read (count from 1).
+    [Default: 1]
+
+zero_ground
+    If true, set HAG of ground-classified points to 0 rather than comparing
+    ``Z`` value to raster DEM.
+    [Default: true]
+

doc/stages/filters.hag_nn.rst

@@ -0,0 +1,131 @@
+.. _filters.hag_nn:
+
+filters.hag_nn
+===============================================================================
+
+The **Height Above Ground Nearest Neighbor filter** takes as input a point
+cloud with ``Classification`` set to 2 for ground points.  It creates a new
+dimension, ``HeightAboveGround``, that contains the normalized height values.
+
+.. note::
+
+   We expect ground returns to have the classification value of 2 in keeping
+   with the `ASPRS Standard LIDAR Point Classes
+   <http://www.asprs.org/a/society/committees/standards/LAS_1_4_r13.pdf>`_.
+
+Ground points may be generated by :ref:`filters.pmf` or :ref:`filters.smrf`,
+but you can use any method you choose, as long as the ground returns are
+marked.
+
+Normalized heights are a commonly used attribute of point cloud data. This can
+also be referred to as *height above ground* (HAG) or *above ground level*
+(AGL) heights. In the end, it is simply a measure of a point's relative height
+as opposed to its raw elevation value.
+
+The filter finds the `count`_ ground points nearest the non-ground point under
+consideration.  It calculates an average ground height weighted by the distance
+of each ground point from the non-ground point.  The ``HeightAboveGround`` is
+the difference between the ``Z`` value of the non-ground point and the
+interpolated ground height.
+
+.. embed::
+
+Example #1
+----------
+
+Using the autzen dataset (here shown colored by elevation), which already has
+points classified as ground
+
+.. image:: ./images/autzen-elevation.png
+   :height: 400px
+
+we execute the following pipeline
+
+.. code-block:: json
+
+  [
+      "autzen.laz",
+      {
+          "type":"filters.hag_nn"
+      },
+      {
+          "type":"writers.laz",
+          "filename":"autzen_hag_nn.laz",
+          "extra_dims":"HeightAboveGround=float32"
+      }
+  ]
+
+which is equivalent to the ``pdal translate`` command
+
+::
+
+    $ pdal translate autzen.laz autzen_hag_nn.laz hag_nn \
+        --writers.las.extra_dims="HeightAboveGround=float32"
+
+In either case, the result, when colored by the normalized height instead of
+elevation is
+
+.. image:: ./images/autzen-hag-nn.png
+   :height: 400px
+
+Example #2
+-------------------------------------------------------------------------------
+
+In the previous example, we chose to write ``HeightAboveGround`` using the
+``extra_dims`` option of :ref:`writers.las`. If you'd instead like to overwrite
+your Z values, then follow the height filter with :ref:`filters.ferry` as shown
+
+.. code-block:: json
+
+  [
+      "autzen.laz",
+      {
+          "type":"filters.hag_nn"
+      },
+      {
+          "type":"filters.ferry",
+          "dimensions":"HeightAboveGround=>Z"
+      },
+      "autzen-height-as-Z.laz"
+  ]
+
+
+Example #3
+-------------------------------------------------------------------------------
+
+If you don't yet have points classified as ground, start with :ref:`filters.pmf`
+or :ref:`filters.smrf` to label ground returns, as shown
+
+.. code-block:: json
+
+  [
+      "autzen.laz",
+      {
+          "type":"filters.smrf"
+      },
+      {
+          "type":"filters.hag_nn"
+      },
+      {
+          "type":"filters.ferry",
+          "dimensions":"HeightAboveGround=>Z"
+      },
+      "autzen-height-as-Z-smrf.laz"
+  ]
+
+Options
+-------------------------------------------------------------------------------
+
+_`count`
+    The number of ground neighbors to consider when determining the height
+    above ground for a non-ground point.  [Default: 1]
+
+max_distance
+    Use only ground points within `max_distance` of non-ground point when
+    performing neighbor interpolation.  [Default: None]
+
+allow_extrapolation
+    If false and a non-ground point lies outside of the bounding box of all
+    ground points, its ``HeightAboveGround`` is set to 0.  If true,
+    extrapolation is used to assign the ``HeightAboveGround`` value.  [Default:
+    false]

doc/stages/filters.rst

@@ -44,7 +44,9 @@ invalidate an existing KD-tree.
    filters.estimaterank
    filters.elm
    filters.ferry
-   filters.hag
+   filters.hag_delaunay
+   filters.hag_dem
+   filters.hag_nn
    filters.info
    filters.lof
    filters.miniball
@@ -106,6 +108,9 @@ invalidate an existing KD-tree.
     Compute pointwise height above ground estimate. Requires points to be
     classified as ground/non-ground prior to estimating.
 
+:ref:`filters.hag_dem`
+    Compute pointwise height above GDAL-readable DEM raster.
+
 :ref:`filters.lof`
     Compute pointwise Local Outlier Factor (along with K-Distance and Local
     Reachability Distance).

filters/DEMFilter.cpp

@@ -102,7 +102,7 @@ void DEMFilter::prepared(PointTableRef table)
     if (m_args->m_dim == Dimension::Id::Unknown)
         throwError("Missing dimension with name '" + m_args->m_range.m_name + "'in input PointView.");
     if (m_args->m_band <= 0)
-        throwError("Band must be greater than 1!");
+        throwError("Band must be greater than 0");
 
 }