Implement coordoperation filter (#2626)

* Implement coordOperation filter * Add filter documentation * Fetch gdal version in cmake using gdal-config and disable coordoperation filter if version is inferior to 3.0.0 * CoordOperation filter modification : - Rename into ProjPipelineFilter - Change a_srs option into out_srs - Delete gdal reference from header - update unit test * udpate documentation
PDAL · Jul 30, 2019 · 47744fd · 47744fd
1 parent 0566c45
commit 47744fd
Show file tree

Hide file tree

Showing 8 changed files with 596 additions and 0 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -240,6 +240,11 @@ if (NOT PDAL_HAVE_LAZPERF)
         ${PDAL_SRC_DIR}/compression/LazPerfVlrCompression.cpp)
     list(REMOVE_ITEM SRCS ${LAZPERF_SRCS})
 endif()
+if (NOT GDAL_3)
+    file(GLOB COORDOP_FILTER_SRCS
+        ${PDAL_FILTERS_DIR}/CoordOperationFilter.cpp)
+    list(REMOVE_ITEM SRCS ${COORDOP_FILTER_SRCS})
+endif()
 PDAL_ADD_LIBRARY(${PDAL_BASE_LIB_NAME} ${SRCS} ${RPLY_SRCS})
 
 #

diff --git a/cmake/gdal.cmake b/cmake/gdal.cmake
@@ -5,6 +5,8 @@ find_package(GDAL 2.2.0)
 set_package_properties(GDAL PROPERTIES TYPE REQUIRED
     PURPOSE "Provides general purpose raster, vector, and reference system support")
 if (GDAL_FOUND)
+    execute_process(COMMAND ${GDAL_CONFIG} "--version" OUTPUT_VARIABLE GDAL_VERSION)
+    string(COMPARE GREATER ${GDAL_VERSION} "3.0.0" GDAL_3)
     mark_as_advanced(CLEAR GDAL_INCLUDE_DIR)
     mark_as_advanced(CLEAR GDAL_LIBRARY)
 else()

diff --git a/doc/stages/filters.projpipeline.rst b/doc/stages/filters.projpipeline.rst
@@ -0,0 +1,85 @@
+.. _filters.projpipeline:
+
+filters.projpipeline
+======================
+
+The projpipeline filter applies a coordinates transformation pipeline. The pipeline could be specified as PROJ string (single step operation or multiple step string starting with +proj=pipeline), a WKT2 string describing a CoordinateOperation, or a “urn:ogc:def:coordinateOperation:EPSG::XXXX” URN.
+
+
+.. note::
+
+    The projpipeline filter does not consider any spatial reference information. 
+    However user could specify an output srs, but no check is done to ensure 
+    the compliance with the provided transformation pipeline.
+
+.. note::
+
+    The projpipeline filter is enabled if the version of GDAL is superior or equal to 3.0
+
+.. streamable::
+
+Example
+-------
+
+This example shift point on the z-axis.
+
+.. code-block:: json
+
+  [
+      "untransformed.las",
+      {
+          "type":"filters.projpipeline",
+          "coord_op":"+proj=affine +zoff=100"
+      },
+      {
+          "type":"writers.las",
+          "filename":"transformed.las"
+      }
+  ]
+
+This example apply a shift on the z-axis then reproject from utm 10
+to WGS84, using the reverse transfo flag. It also set the output srs
+
+.. code-block:: json
+
+
+  [ 
+      "utm10.las", 
+      { 
+          "type":"filters.projpipeline", 
+          "coord_op":"+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=10 +step +proj=affine +zoff=100",
+          "reverse_transfo": "true", 
+          "out_srs": "EPSG:4326" 
+      }, 
+      { 
+          "type":"writers.las", 
+          "filename":"wgs84.las" 
+      } 
+  ]
+
+.. note::
+
+    PDAL use the GDAL `OGRCoordinateTransformation` class to transform coordinates. 
+    By default output angular unit are in radians. To change to degrees we need to 
+    apply a unit conversion step.
+
+
+
+Options
+-------
+
+_`coord_op`
+  The coordinate operation string.
+  Could be specified as PROJ string (single step operation or 
+  multiple step string starting with +proj=pipeline), 
+  a WKT2 string describing a CoordinateOperation, 
+  or a “urn:ogc:def:coordinateOperation:EPSG::XXXX” URN.
+_`reverse_transfo`
+  Boolean, Whether the coordinate operation should be evaluated
+  in the reverse path [Default: false]
+_`out_srs`
+  The spatial reference system of the file to be written. 
+  Can be an EPSG string (e.g. “EPSG:26910”) or a WKT string. 
+  No check is done to ensure the compliance with the specified coordinate
+  operation [Default: Not set]
+
diff --git a/doc/stages/filters.rst b/doc/stages/filters.rst
@@ -158,6 +158,7 @@ PDAL filters that move XYZ coordinates will invalidate an existing KD-tree.
 
    filters.cpd
    filters.icp
+   filters.projpipeline
    filters.reprojection
    filters.transformation
 
@@ -169,6 +170,10 @@ PDAL filters that move XYZ coordinates will invalidate an existing KD-tree.
     Compute and apply transformation between two point clouds using the
     Iterative Closest Point algorithm.
 
+:ref:`filters.projpipeline`
+    Apply coordinates operation on point triplets, based on PROJ pipeline string,
+    WKT2 coordinates operations or URN definitions.
+
 :ref:`filters.reprojection`
     Reproject data using GDAL from one coordinate system to another.
 

diff --git a/filters/ProjPipelineFilter.cpp b/filters/ProjPipelineFilter.cpp
@@ -0,0 +1,137 @@
+/******************************************************************************
+* Copyright (c) 2019, Aurelien Vila (aurelien.vila@delair.aero)
+*
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without
+* modification, are permitted provided that the following
+* conditions are met:
+*
+*     * Redistributions of source code must retain the above copyright
+*       notice, this list of conditions and the following disclaimer.
+*     * Redistributions in binary form must reproduce the above copyright
+*       notice, this list of conditions and the following disclaimer in
+*       the documentation and/or other materials provided
+*       with the distribution.
+*     * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
+*       names of its contributors may be used to endorse or promote
+*       products derived from this software without specific prior
+*       written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+* OF SUCH DAMAGE.
+****************************************************************************/
+
+#include "ProjPipelineFilter.hpp"
+
+#include <pdal/PointView.hpp>
+#include <pdal/private/SrsTransform.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include <ogr_spatialref.h>
+
+namespace pdal
+{
+
+static StaticPluginInfo const s_info
+{
+    "filters.projpipeline",
+    "Transform coordinates using Proj pipeline string, WKT2 coordinate operations or URN definition",
+    "http://pdal.io/stages/filters.projpipeline.html"
+};
+
+CREATE_STATIC_STAGE(ProjPipelineFilter, s_info)
+
+std::string ProjPipelineFilter::getName() const { return s_info.name; }
+
+ProjPipelineFilter::ProjPipelineFilter()
+{}
+
+
+ProjPipelineFilter::~ProjPipelineFilter()
+{}
+
+
+void ProjPipelineFilter::addArgs(ProgramArgs& args)
+{
+    args.add("out_srs", "Output spatial reference", m_outSRS);
+    args.add("reverse_transfo", "Wether the coordinate operation should be evaluated in the reverse path",
+             m_reverseTransfo, false);
+    args.add("coord_op", "Coordinate operation (Proj pipeline or WKT2 string or urn definition)",
+             m_coordOperation).setPositional();
+}
+
+
+void ProjPipelineFilter::initialize()
+{
+    setSpatialReference(m_outSRS);
+    createTransform(m_coordOperation, m_reverseTransfo);
+}
+
+
+void ProjPipelineFilter::createTransform(const std::string coordOperation, bool reverseTransfo)
+{
+    m_coordTransform.reset(new CoordTransform(coordOperation, reverseTransfo));
+}
+
+
+PointViewSet ProjPipelineFilter::run(PointViewPtr view)
+{
+    PointViewSet viewSet;
+    PointViewPtr outView = view->makeNew();
+
+    PointRef point(*view, 0);
+    for (PointId id = 0; id < view->size(); ++id)
+    {
+        point.setPointId(id);
+        if (processOne(point))
+            outView->appendPoint(*view, id);
+    }
+
+    viewSet.insert(outView);
+    return viewSet;
+}
+
+
+bool ProjPipelineFilter::processOne(PointRef& point)
+{
+    double x(point.getFieldAs<double>(Dimension::Id::X));
+    double y(point.getFieldAs<double>(Dimension::Id::Y));
+    double z(point.getFieldAs<double>(Dimension::Id::Z));
+
+    bool ok = m_coordTransform->transform(x, y, z);
+    if (ok)
+    {
+        point.setField(Dimension::Id::X, x);
+        point.setField(Dimension::Id::Y, y);
+        point.setField(Dimension::Id::Z, z);
+    }
+    return ok;
+}
+
+ProjPipelineFilter::CoordTransform::CoordTransform(){}
+
+ProjPipelineFilter::CoordTransform::CoordTransform(const std::string coordOperation, bool reverseTransfo){
+    OGRCoordinateTransformationOptions coordTransfoOptions;
+    coordTransfoOptions.SetCoordinateOperation(coordOperation.c_str(), reverseTransfo);
+    OGRSpatialReference nullSrs("");
+    m_transform.reset(OGRCreateCoordinateTransformation(&nullSrs, &nullSrs, coordTransfoOptions));
+}
+
+bool ProjPipelineFilter::CoordTransform::transform(double &x, double &y, double &z){
+
+    return m_transform && m_transform->Transform(1, &x, &y, &z);
+}
+
+} // namespace pdal
+
diff --git a/filters/ProjPipelineFilter.hpp b/filters/ProjPipelineFilter.hpp
@@ -0,0 +1,89 @@
+/******************************************************************************
+* Copyright (c) 2019, Aurelien Vila (aurelien.vila@delair.aero)
+*
+* All rights reserved.
+*
+* Redistribution and use in source and binary forms, with or without
+* modification, are permitted provided that the following
+* conditions are met:
+*
+*     * Redistributions of source code must retain the above copyright
+*       notice, this list of conditions and the following disclaimer.
+*     * Redistributions in binary form must reproduce the above copyright
+*       notice, this list of conditions and the following disclaimer in
+*       the documentation and/or other materials provided
+*       with the distribution.
+*     * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
+*       names of its contributors may be used to endorse or promote
+*       products derived from this software without specific prior
+*       written permission.
+*
+* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+* OF SUCH DAMAGE.
+****************************************************************************/
+
+#pragma once
+
+#include <pdal/Filter.hpp>
+#include <pdal/Streamable.hpp>
+
+#include <memory>
+
+class OGRCoordinateTransformation;
+
+namespace pdal
+{
+
+
+class PDAL_DLL ProjPipelineFilter : public Filter, public Streamable
+{
+public:
+    class CoordTransform;
+
+    ProjPipelineFilter();
+    ~ProjPipelineFilter();
+
+    std::string getName() const;
+
+private:
+    ProjPipelineFilter& operator=(const ProjPipelineFilter&) = delete;
+    ProjPipelineFilter(const ProjPipelineFilter&) = delete;
+
+    virtual void addArgs(ProgramArgs& args);
+    virtual void initialize();
+    virtual PointViewSet run(PointViewPtr view);
+    virtual bool processOne(PointRef& point);
+
+    void createTransform(const std::string coordOperation, bool reverseTransfo);
+
+    SpatialReference m_outSRS;
+    bool m_reverseTransfo;
+    std::string m_coordOperation;
+    std::unique_ptr<CoordTransform> m_coordTransform;
+};
+
+
+class ProjPipelineFilter::CoordTransform
+{
+public:
+    CoordTransform();
+    CoordTransform(const std::string coordOperation, bool reverseTransfo);
+
+    bool transform(double& x, double& y, double& z);
+private:
+    std::unique_ptr<OGRCoordinateTransformation> m_transform;
+
+};
+
+} // namespace pdal
+
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
@@ -263,6 +263,9 @@ PDAL_ADD_TEST(pdal_filters_overlay_test FILES filters/OverlayFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_pmf_test FILES filters/PMFFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_reprojection_test FILES
     filters/ReprojectionFilterTest.cpp)
+if (GDAL_3)
+    PDAL_ADD_TEST(pdal_filters_projpipeline_test FILES filters/ProjPipelineFilterTest.cpp)
+endif()
 PDAL_ADD_TEST(pdal_filters_range_test FILES filters/RangeFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_randomize_test FILES filters/RandomizeFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_reciprocity_test FILES filters/ReciprocityFilterTest.cpp)