Merge remote-tracking branch 'origin/master' into issue-2852

CMakeLists.txt

@@ -5,7 +5,7 @@
 
 cmake_minimum_required(VERSION 3.5)
 
-project(PDAL VERSION 2.0.0 LANGUAGES CXX C)
+project(PDAL VERSION 2.1.0 LANGUAGES CXX C)
 string(TOLOWER ${PROJECT_NAME} PROJECT_NAME_LOWER)
 
 set(ROOT_DIR "${PROJECT_SOURCE_DIR}")

cmake/gtest.cmake

@@ -4,14 +4,3 @@ if (MSVC)
     option(gtest_force_shared_crt "Always use shared Visual C++ run-time DLL" ON)
 endif()
 add_subdirectory(vendor/gtest)
-
-# gtest 1.7.0 has some CMake warnings, so we silence these by setting the
-# following properties.
-set_target_properties(gtest PROPERTIES
-    MACOSX_RPATH ON
-    LIBRARY_OUTPUT_DIRECTORY "${gtest_BINARY_DIR}/src"
-)
-set_target_properties(gtest_main PROPERTIES
-    MACOSX_RPATH ON
-    LIBRARY_OUTPUT_DIRECTORY "${gtest_BINARY_DIR}/src"
-    )

doc/stages/filters.python.rst

@@ -12,6 +12,15 @@ The function must have two `NumPy`_ arrays as arguments, ``ins`` and ``outs``.
 The ``ins`` array represents the points before the ``filters.python``
 filter and the ``outs`` array represents the points after filtering.
 
+.. warning::
+
+    Make sure `NumPy`_ is installed in your `Python`_ environment.
+
+    .. code-block:: shell
+
+        $ python3 -c "import numpy; print(numpy.__version__)"
+        1.18.1
+
 .. warning::
 
     Each array contains all the :ref:`dimensions` of the incoming ``ins``
@@ -102,14 +111,14 @@ classified 1 or 2 to be dropped from the point stream.
   def filter(ins,outs):
      cls = ins['Classification']
 
-     keep_classes = [1,2]
+     keep_classes = [1, 2]
 
      # Use the first test for our base array.
      keep = np.equal(cls, keep_classes[0])
 
      # For 1:n, test each predicate and join back
      # to our existing predicate array
-     for k in range(1,len(keep_classes)):
+     for k in range(1, len(keep_classes)):
          t = np.equal(cls, keep_classes[k])
          keep = keep + t
 
@@ -164,9 +173,9 @@ for the in-scope :ref:`filters.python` :cpp:class:`pdal::Stage`.
 .. code-block:: python
 
    def myfunc(ins,outs):
-       print ('schema: ', schema)
-       print ('srs: ', spatialreference)
-       print ('metadata: ', metadata)
+       print('schema: ', schema)
+       print('srs: ', spatialreference)
+       print('metadata: ', metadata)
        outs = ins
        return True
 
@@ -207,7 +216,7 @@ basic Python function while substituting values as necessary.
           "type":"filters.python",
           "module":"anything",
           "function":"filter",
-          "source":"arguments.py",
+          "script":"arguments.py",
           "pdalargs":"{\"factor\":0.3048,\"an_argument\":42, \"another\": \"a string\"}"
       },
       "output.las"

doc/stages/readers.text.rst

@@ -9,10 +9,24 @@ into a number of fields by a separator.  Each field represents a value for
 a point's dimension.  Each value needs to be `formatted`_ properly for
 C++ language double-precision values.
 
-The text reader expects a header line to 1) indicate the separator character
-for the fields and 2) name the point dimension for each field.  Any
-single non-alphanumeric character can be used as a separator.  The header line
-separator can be overridden by the `separator`_ option.  Each line in the
+The text reader expects a header line to indicate the dimensions are
+in each subsequent line.  There are two types of header lines.
+
+Quoted dimension names
+----------------------
+When the first character of the header is a double quote, each dimension name
+is assumed to be surrounded by double quotes.  Any text following a quoted
+dimension name and the start of the next dimension name is ignored.  The
+`separator`_ option can't be used with quoted dimension names.
+
+Unquoted dimension names
+------------------------
+
+The first non alpha-numeric character encountered is treated as a separator
+between dimension names.  The separator in the header line can be overridden
+by the `separator`_ option.
+
+Each line in the
 file must contain the same number of fields as indicated by
 dimension names in the header.  Spaces are generally ignored in the input
 unless used as a separator.  When a space character is used as a separator,
@@ -44,8 +58,8 @@ This input file contains X, Y and Z value for 10 points.
     289818.01,4320980.38,170.61
     289818.50,4320980.59,170.58
 
-Example Pipeline
-----------------
+Example #1
+----------
 
 .. code-block:: json
 
@@ -60,6 +74,27 @@ Example Pipeline
       }
   ]
 
+Example #2
+----------
+
+This reads the data in the input file as Red, Green and Blue instead of
+as X, Y and Z.
+
+.. code-block:: json
+
+  [
+      {
+          "type":"readers.text",
+          "filename":"inputfile.txt",
+          "header":"Red, Green, Blue",
+          "skip":1
+      },
+      {
+          "type":"writers.text",
+          "filename":"outputfile.txt"
+      }
+  ]
+
 
 Options
 -------

doc/stages/writers.las.rst

@@ -177,7 +177,7 @@ offset_x, offset_y, offset_z
 
 filesource_id
   The file source id number to use for this file (a value between
-  1 and 65535) [Default: 0]
+  0 and 65535 - 0 implies "unassigned") [Default: 0]
 
 discard_high_return_numbers
   If true, discard all points with a return number greater than the maximum

examples/writing-reader/MyReader.cpp

@@ -31,6 +31,7 @@ namespace pdal
 
   void MyReader::ready(PointTableRef)
   {
+    m_index = 0;
     SpatialReference ref("EPSG:4385");
     setSpatialReference(ref);
   }

filters/ELMFilter.cpp

@@ -64,7 +64,7 @@ std::string ELMFilter::getName() const
 void ELMFilter::addArgs(ProgramArgs& args)
 {
     args.add("cell", "Cell size", m_cell, 10.0);
-    args.add("class", "Class to use for noise points", m_class, uint8_t(7));
+    args.add("class", "Class to use for noise points", m_class, ClassLabel::LowPoint);
     args.add("threshold", "Threshold value", m_threshold, 1.0);
 }
 

filters/HAGFilter.cpp

@@ -235,7 +235,7 @@ void HAGFilter::filter(PointView& view)
     // Separate into ground and non-ground views.
     for (PointId i = 0; i < view.size(); ++i)
     {
-        if (view.getFieldAs<short>(Id::Classification, i) == 2)
+        if (view.getFieldAs<uint8_t>(Id::Classification, i) == ClassLabel::Ground)
         {
             view.setField(Id::HeightAboveGround, i, 0);
             gView->appendPoint(view, i);

filters/HexBinFilter.cpp

@@ -83,7 +83,8 @@ void HexBin::addArgs(ProgramArgs& args)
     m_cullArg = &args.add("hole_cull_area_tolerance", "Tolerance area to "
         "apply to holes before cull", m_cullArea);
     args.add("smooth", "Smooth boundary output", m_doSmooth, true);
-    args.add("preserve_topology", "Preserve topology when smoothing", m_preserve_topology, true);
+    args.add("preserve_topology", "Preserve topology when smoothing",
+        m_preserve_topology, true);
 }
 
 
@@ -189,6 +190,8 @@ void HexBin::done(PointTableRef table)
             "Boundary MULTIPOLYGON of domain");
     }
 
+    SpatialReference srs(table.anySpatialReference());
+    pdal::Polygon p(polygon.str(), srs);
 
     /***
       We want to make these bumps on edges go away, which means that
@@ -205,49 +208,34 @@ void HexBin::done(PointTableRef table)
        A /        \ D
 
     ***/
-    double tolerance = 1.1 * m_grid->height() / 2;
-
-    double cull = m_cullArg->set() ? m_cullArea : (6 * tolerance * tolerance);
-
-    SpatialReference srs(table.anySpatialReference());
-    pdal::Polygon p(polygon.str(), srs);
-    pdal::Polygon density_p(polygon.str(), srs);
+    if (m_doSmooth)
+    {
+        double tolerance = 1.1 * m_grid->height() / 2;
+        double cull = m_cullArg->set() ?
+            m_cullArea : (6 * tolerance * tolerance);
+        p.simplify(tolerance, cull, m_preserve_topology);
+    }
 
     // If the SRS was geographic, use relevant
     // UTM for area and density computation
+    Polygon density_p(p);
     if (srs.isGeographic())
     {
         // Compute a UTM polygon
         BOX3D box = p.bounds();
         int zone = SpatialReference::calculateZone(box.minx, box.miny);
 
-        auto makezone = [] (int zone) -> std::string
-        {
-
-            std::ostringstream z;
-
-            // Use WGS84 UTM zones
-            z << "EPSG:327" << abs(zone);
-            return z.str();
-        };
-
-        SpatialReference utm(makezone(zone));
-        density_p = p;
+        SpatialReference utm("EPSG:327" + std::to_string(std::abs(zone)));
         density_p.transform(utm);
     }
 
-    if (m_doSmooth)
-        p.simplify(tolerance, cull, m_preserve_topology);
-
-    std::string boundary_text = p.wkt(m_precision);
-
-    m_metadata.add("boundary", boundary_text,
+    m_metadata.add("boundary", p.wkt(m_precision),
         "Approximated MULTIPOLYGON of domain");
     m_metadata.addWithType("boundary_json", p.json(), "json",
         "Approximated MULTIPOLYGON of domain");
-    double area = density_p.area();
 
-    double density = (double) m_count/ area ;
+    double area = density_p.area();
+    double density = m_count / area;
     if (std::isinf(density))
     {
         density = -1.0;
@@ -257,29 +245,22 @@ void HexBin::done(PointTableRef table)
     m_metadata.add("density", density,
         "Number of points per square unit (total area)");
     m_metadata.add("area", area, "Area in square units of tessellated polygon");
+    m_metadata.add("avg_pt_spacing", std::sqrt(1 / density),
+        "Avg point spacing (x/y units)");
 
-    double moving_avg(0.0);
-    double avg_count(0.0);
-
-    double hex_area(((3 * SQRT_3)/2.0) * (m_grid->height() * m_grid->height()));
     int n(0);
     point_count_t totalCount(0);
-    double totalArea(0.0);
     for (HexIter hi = m_grid->hexBegin(); hi != m_grid->hexEnd(); ++hi)
     {
         HexInfo h = *hi;
         totalCount += h.density();
-        totalArea += hex_area;
         ++n;
     }
 
-    double avg_density = totalArea /(double) totalCount;
+    double hexArea(((3 * SQRT_3)/2.0) * (m_grid->height() * m_grid->height()));
+    double avg_density = (n * hexArea) / totalCount;
     m_metadata.add("avg_pt_per_sq_unit", avg_density, "Number of points "
         "per square unit (tessellated area within inclusions)");
-
-    double avg_spacing = std::sqrt(1/density);
-    m_metadata.add("avg_pt_spacing", avg_spacing,
-        "Avg point spacing (x/y units)");
 }
 
 } // namespace pdal

filters/OutlierFilter.cpp

@@ -66,7 +66,7 @@ void OutlierFilter::addArgs(ProgramArgs& args)
     args.add("radius", "Radius", m_radius, 1.0);
     args.add("mean_k", "Mean number of neighbors", m_meanK, 8);
     args.add("multiplier", "Standard deviation threshold", m_multiplier, 2.0);
-    args.add("class", "Class to use for noise points", m_class, uint8_t(7));
+    args.add("class", "Class to use for noise points", m_class, ClassLabel::LowPoint);
 }
 
 void OutlierFilter::addDimensions(PointLayoutPtr layout)

filters/PMFFilter.cpp

@@ -204,7 +204,7 @@ PointViewSet PMFFilter::run(PointViewPtr input)
     // Classify remaining points with value of 1. processGround will mark ground
     // returns as 2.
     for (PointId i = 0; i < secondView->size(); ++i)
-        secondView->setField(Dimension::Id::Classification, i, 1);
+        secondView->setField(Dimension::Id::Classification, i, ClassLabel::Unclassified);
 
     // Run the actual PMF algorithm.
     processGround(firstView);
@@ -369,7 +369,7 @@ void PMFFilter::processGround(PointViewPtr view)
     // set the classification label of ground returns as 2
     // (corresponding to ASPRS LAS specification)
     for (const auto& i : groundIdx)
-        view->setField(Dimension::Id::Classification, i, 2);
+        view->setField(Dimension::Id::Classification, i, ClassLabel::Ground);
 }
 
 } // namespace pdal

filters/SMRFilter.cpp

@@ -133,7 +133,7 @@ void SMRFilter::prepared(PointTableRef table)
     for (auto& r : m_args->m_ignored)
     {
         r.m_id = layout->findDim(r.m_name);
-        if (r.m_id == Dimension::Id::Unknown)
+        if (r.m_id == Id::Unknown)
             throwError("Invalid dimension name in 'ignored' option: '" +
                        r.m_name + "'.");
     }
@@ -149,8 +149,8 @@ void SMRFilter::prepared(PointTableRef table)
             }
         }
 
-        if (!layout->hasDim(Dimension::Id::ReturnNumber) ||
-            !layout->hasDim(Dimension::Id::NumberOfReturns))
+        if (!layout->hasDim(Id::ReturnNumber) ||
+            !layout->hasDim(Id::NumberOfReturns))
         {
             log()->get(LogLevel::Warning) << "Could not find ReturnNumber and "
                                              "NumberOfReturns. Skipping "
@@ -193,9 +193,9 @@ PointViewSet SMRFilter::run(PointViewPtr view)
     for (PointId i = 0; i < keptView->size(); ++i)
     {
         uint8_t nr =
-            keptView->getFieldAs<uint8_t>(Dimension::Id::NumberOfReturns, i);
+            keptView->getFieldAs<uint8_t>(Id::NumberOfReturns, i);
         uint8_t rn =
-            keptView->getFieldAs<uint8_t>(Dimension::Id::ReturnNumber, i);
+            keptView->getFieldAs<uint8_t>(Id::ReturnNumber, i);
         if ((nr == 0) && !nrOneZero)
             nrOneZero = true;
         if ((rn == 0) && !rnOneZero)
@@ -233,7 +233,7 @@ PointViewSet SMRFilter::run(PointViewPtr view)
     }
 
     for (PointId i = 0; i < secondView->size(); ++i)
-        secondView->setField(Dimension::Id::Classification, i, 1);
+        secondView->setField(Id::Classification, i, ClassLabel::Unclassified);
 
     m_srs = firstView->spatialReference();
 
@@ -362,9 +362,9 @@ void SMRFilter::classifyGround(PointViewPtr view, std::vector<double>& ZIpro)
         // vertical distance between each LIDAR point and the provisional
         // DEM, and applying a threshold calculation."
         if (std::fabs(ZIpro[c * m_rows + r] - z) > thresh(r, c))
-            view->setField(Id::Classification, i, 1);
+            view->setField(Id::Classification, i, ClassLabel::Unclassified);
         else
-            view->setField(Id::Classification, i, 2);
+            view->setField(Id::Classification, i, ClassLabel::Ground);
     }
 }
 
@@ -448,8 +448,6 @@ std::vector<int> SMRFilter::createObjMask(std::vector<double> const& ZImin)
 
 std::vector<double> SMRFilter::createZImin(PointViewPtr view)
 {
-    using namespace Dimension;
-
     // "As with many other ground filtering algorithms, the first step is
     // generation of ZImin from the cell size parameter and the extent of the
     // data."