r.in.pdal: Add COPC support

raster/r.in.pdal/grassrasterwriter.h

@@ -37,7 +37,14 @@ class GrassRasterWriter : public pdal::NoFilenameWriter,
 class GrassRasterWriter : public pdal::Writer, public pdal::Streamable {
 #endif
 public:
-    GrassRasterWriter() : n_processed(0) {}
+    GrassRasterWriter()
+        : n_processed(0), n_on_edge(0), region_(nullptr),
+          point_binning_(nullptr), bin_index_nodes_(nullptr),
+          rtype_(FCELL_TYPE), cols_(0), scale_(1.0),
+          dim_to_import_(pdal::Dimension::Id::Z), base_segment_(nullptr),
+          input_region_(nullptr), base_raster_data_type_(FCELL_TYPE)
+    {
+    }
 
     std::string getName() const { return "writers.grassbinning"; }
 
@@ -94,17 +101,11 @@ class GrassRasterWriter : public pdal::Writer, public pdal::Streamable {
             z -= base_z;
         }
 
-        // TODO: check the bounds and report discrepancies in
-        // number of filtered out vs processed to the user
-        // (alternatively, change the spatial bounds test to
-        // give same results as this, but it might be actually helpful
-        // to tell user that they have points right on the border)
         int arr_row = (int)((region_->north - y) / region_->ns_res);
         int arr_col = (int)((x - region_->west) / region_->ew_res);
 
         if (arr_row >= region_->rows || arr_col >= region_->cols) {
-            G_message(_("A point on the edge of computational region detected. "
-                        "Ignoring."));
+            n_on_edge++;
             return false;
         }
 
@@ -115,6 +116,7 @@ class GrassRasterWriter : public pdal::Writer, public pdal::Streamable {
     }
 
     gpoint_count n_processed;
+    gpoint_count n_on_edge;
 
 private:
     struct Cell_head *region_;

raster/r.in.pdal/info.cpp

@@ -275,6 +275,7 @@ void print_lasinfo(struct StringList *infiles)
         const pdal::LasHeader &h = las_reader.header();
         pdal::PointLayoutPtr point_layout = table.layout();
         const pdal::Dimension::IdList &dims = point_layout->dims();
+        std::string proj_wkt = las_reader.getSpatialReference().getWKT();
 
         std::cout << "File: " << infile << std::endl;
         std::cout << "File version = "
@@ -304,6 +305,10 @@ void print_lasinfo(struct StringList *infiles)
                   << h.maxZ() << "\n";
         std::cout << "Min X/Y/Z: " << h.minX() << "/" << h.minY() << "/"
                   << h.minZ() << "\n";
+        if (!proj_wkt.empty())
+            std::cout << "Projection (WKT): " << proj_wkt << "\n";
+        else
+            std::cout << "Projection: (undefined)\n";
         if (h.versionAtLeast(1, 4)) {
             std::cout << "Ext. VLR offset: " << h.eVlrOffset() << "\n";
             std::cout << "Ext. VLR count: " << h.eVlrCount() << "\n";

raster/r.in.pdal/main.cpp

@@ -16,6 +16,8 @@
  *****************************************************************************/
 
 #include <cstdio>
+#include <iomanip>
+#include <sstream>
 
 #if defined(__clang__)
 #pragma clang diagnostic push
@@ -368,6 +370,20 @@ int main(int argc, char *argv[])
     user_dimension_opt->description = _("PDAL dimension name");
     user_dimension_opt->guisection = _("Selection");
 
+    Option *point_table_capacity_opt = G_define_option();
+
+    point_table_capacity_opt->key = "point_table_capacity";
+    point_table_capacity_opt->type = TYPE_INTEGER;
+    point_table_capacity_opt->required = NO;
+    point_table_capacity_opt->answer = const_cast<char *>("10000");
+    point_table_capacity_opt->options = "1-";
+    point_table_capacity_opt->label =
+        _("Number of points buffered at once during processing");
+    point_table_capacity_opt->description =
+        _("Larger values may improve performance for large datasets at the "
+          "cost of memory");
+    point_table_capacity_opt->guisection = _("Performance");
+
     Flag *extents_flag = G_define_flag();
 
     extents_flag->key = 'e';
@@ -721,6 +737,7 @@ int main(int argc, char *argv[])
     pdal::StageFactory factory;
     pdal::MergeFilter merge_filter;
     bool need_to_reproject = false;
+    bool bounds_pushed_to_reader = false;
     /* loop of input files */
     for (int i = 0; i < infiles.num_items; i++) {
         const char *infile = infiles.items[i];
@@ -748,6 +765,8 @@ int main(int argc, char *argv[])
         reader->setOptions(las_opts);
         readers.push_back(reader);
 
+        // With the -o flag the CRS is assumed to match the project's CRS.
+        bool proj_match = true;
         // getting projection is possible only after prepare
         if (!over_flag->answer) {
             pdal::PointTable table;
@@ -763,12 +782,27 @@ int main(int argc, char *argv[])
             if (spatial_reference.empty())
                 G_fatal_error(_("The input dataset has undefined projection"));
             std::string dataset_wkt = spatial_reference.getWKT();
-            bool proj_match =
-                is_wkt_projection_same_as_loc(dataset_wkt.c_str());
+            proj_match = is_wkt_projection_same_as_loc(dataset_wkt.c_str());
 
             if (!proj_match)
                 need_to_reproject = true;
         }
+        // Let the COPC octree skip nodes outside the region; the region
+        // bounds are in the project's CRS, so this applies only when the
+        // file's CRS matches. The GRASS filter still does the exact clip.
+        // The option takes effect when the whole pipeline is prepared later.
+        // With -e the region is the extent of the data, so there is nothing
+        // outside it to skip.
+        if (use_spatial_filter && !extents_flag->answer && proj_match &&
+            pdal_read_driver == "readers.copc") {
+            std::ostringstream bounds_str;
+            bounds_str << std::setprecision(17) << "([" << xmin << ", " << xmax
+                       << "], [" << ymin << ", " << ymax << "])";
+            pdal::Options bounds_opts;
+            bounds_opts.add("bounds", bounds_str.str());
+            reader->addOptions(bounds_opts);
+            bounds_pushed_to_reader = true;
+        }
         merge_filter.setInput(*reader);
     }
     if (over_flag->answer) {
@@ -825,9 +859,10 @@ int main(int argc, char *argv[])
     binning_writer.set_output_scale(output_scale);
     binning_writer.setInput(grass_filter);
     // stream_filter.setInput(*last_stage);
-    //  there is no difference between 1 and 10k points in memory
-    //  consumption, so using 10k in case it is faster for some cases
-    pdal::point_count_t point_table_capacity = 10000;
+    // The default capacity of 10k points takes no more memory than 1, but
+    // can be faster; larger values trade memory for speed.
+    pdal::point_count_t point_table_capacity =
+        atoi(point_table_capacity_opt->answer);
     pdal::FixedPointTable point_table(point_table_capacity);
     try {
         binning_writer.prepare(point_table);
@@ -959,8 +994,16 @@ int main(int argc, char *argv[])
     }
 
     G_done_msg("%s", buff);
-    G_message("Filtered spatially " GPOINT_COUNT_FORMAT " points.",
-              grass_filter.num_spatially_filtered());
+    // Points pruned by the reader's spatial index never reach the filter, so
+    // the count below is not comparable to a run without the index.
+    if (bounds_pushed_to_reader)
+        G_message("Filtered spatially " GPOINT_COUNT_FORMAT
+                  " points (excluding those skipped by the spatial index of "
+                  "the input).",
+                  grass_filter.num_spatially_filtered());
+    else
+        G_message("Filtered spatially " GPOINT_COUNT_FORMAT " points.",
+                  grass_filter.num_spatially_filtered());
     G_message("Filtered z range " GPOINT_COUNT_FORMAT " points.",
               grass_filter.num_zrange_filtered());
     G_message("Filtered i range " GPOINT_COUNT_FORMAT " points.",
@@ -972,6 +1015,10 @@ int main(int argc, char *argv[])
     G_message("Filtered return " GPOINT_COUNT_FORMAT " points.",
               grass_filter.num_return_filtered());
 
+    if (binning_writer.n_on_edge)
+        G_message("Skipped " GPOINT_COUNT_FORMAT
+                  " points on the edge of the computational region.",
+                  binning_writer.n_on_edge);
     G_message("Processed into raster " GPOINT_COUNT_FORMAT " points.",
               binning_writer.n_processed);
 

raster/r.in.pdal/r.in.pdal.html

@@ -359,7 +359,13 @@ <h3>Format and projection support</h3>
 
 <em>r.in.pdal</em> reads all point cloud formats supported by the PDAL library.
 The typical file extensions for the LAS format are <code>.las</code> and
-<code>.laz</code> (compressed).
+<code>.laz</code> (compressed). Reading compressed LAS (<code>.laz</code>)
+requires PDAL built with either LASzip or lazperf. Cloud
+Optimized Point Clouds (COPC, <code>.copc.laz</code>) are supported as well.
+When importing a COPC file into a computational region smaller than the point
+cloud extent, the COPC spatial index is used to skip data outside the region,
+which speeds up the import; this optimization does not apply when the points
+are reprojected.
 
 <p>
 The coordinate reference system (CRS) of the input is read from the file
@@ -438,6 +444,11 @@ <h3>Memory consumption</h3>
 The default map <b>type</b>=<code>FCELL</code> is intended as compromise between
 preserving data precision and limiting system resource consumption.
 
+<p>
+The <b>point_table_capacity</b> option controls how many points PDAL
+buffers at once while streaming the input. Larger values may speed up
+the import of large files at the cost of additional memory.
+
 <h3>Trim option</h3>
 <p>
 Trim option value is used only when calculating trimmed mean values.

raster/r.in.pdal/r.in.pdal.md

@@ -324,7 +324,13 @@ parameter is provided.
 
 *r.in.pdal* reads all point cloud formats supported by the PDAL library.
 The typical file extensions for the LAS format are `.las` and `.laz`
-(compressed).
+(compressed). Reading compressed LAS (`.laz`) requires PDAL built with
+either LASzip or lazperf. Cloud Optimized Point Clouds
+(COPC, `.copc.laz`) are supported as well. When importing a COPC file
+into a computational region smaller than the point cloud extent, the
+COPC spatial index is used to skip data outside the region, which speeds
+up the import; this optimization does not apply when the points are
+reprojected.
 
 The coordinate reference system (CRS) of the input is read from the file
 metadata and compared with the CRS of the current project (previously
@@ -369,6 +375,10 @@ memory use for these also depends on the number of data points.
 The default map **type**=`FCELL` is intended as compromise between
 preserving data precision and limiting system resource consumption.
 
+The **point_table_capacity** option controls how many points PDAL
+buffers at once while streaming the input. Larger values may speed up
+the import of large files at the cost of additional memory.
+
 ### Trim option
 
 Trim option value is used only when calculating trimmed mean values.

raster/r.in.pdal/tests/conftest.py

@@ -0,0 +1,67 @@
+import os
+import shutil
+import subprocess
+
+import pytest
+
+import grass.script as gs
+
+
+@pytest.fixture(scope="module")
+def point_cloud_files(tmp_path_factory):
+    """Create LAS and COPC versions of a synthetic point cloud.
+
+    Points form an 18x18 grid with 1 unit spacing starting at 0.5 and
+    z = x + y, georeferenced as EPSG:3358.
+    """
+    if shutil.which("pdal") is None:
+        pytest.skip("pdal command line tool not available")
+    tmp_path = tmp_path_factory.mktemp("point_data")
+    csv_file = tmp_path / "points.csv"
+    las_file = tmp_path / "points.las"
+    copc_file = tmp_path / "points.copc.laz"
+    lines = ["X,Y,Z"]
+    for i in range(18):
+        for j in range(18):
+            x = i + 0.5
+            y = j + 0.5
+            lines.append(f"{x},{y},{x + y}")
+    csv_file.write_text("\n".join(lines) + "\n")
+    subprocess.run(
+        [
+            "pdal",
+            "translate",
+            "-i",
+            str(csv_file),
+            "-o",
+            str(las_file),
+            "-r",
+            "text",
+            "-w",
+            "las",
+            "--writers.las.a_srs=EPSG:3358",
+        ],
+        check=True,
+        capture_output=True,
+    )
+    try:
+        subprocess.run(
+            ["pdal", "translate", "-i", str(las_file), "-o", str(copc_file)],
+            check=True,
+            capture_output=True,
+        )
+    except subprocess.CalledProcessError:
+        pytest.skip("PDAL does not support writing COPC")
+    return {"las": las_file, "copc": copc_file}
+
+
+@pytest.fixture(scope="module")
+def session(tmp_path_factory):
+    """Active session in a project with the CRS of the test data"""
+    tmp_path = tmp_path_factory.mktemp("r_in_pdal_project")
+    project = tmp_path / "test"
+    gs.create_project(project, epsg="3358")
+    with gs.setup.init(project, env=os.environ.copy()) as session:
+        if shutil.which("r.in.pdal", path=session.env["PATH"]) is None:
+            pytest.skip("r.in.pdal not available (GRASS built without PDAL)")
+        yield session

raster/r.in.pdal/tests/r_in_pdal_copc_test.py

@@ -0,0 +1,94 @@
+"""Tests of r.in.pdal COPC support and the point_table_capacity option"""
+
+import pytest
+
+from grass.exceptions import CalledModuleError
+from grass.tools import Tools
+
+
+def test_copc_import(session, point_cloud_files):
+    """SRS of a COPC file is detected and all points are binned"""
+    tools = Tools(session=session)
+    tools.g_region(n=18, s=0, e=18, w=0, res=6)
+    tools.r_in_pdal(input=point_cloud_files["copc"], output="copc_n", method="n")
+    stats = tools.r_univar(map="copc_n", flags="g").keyval
+    assert stats["n"] == 9
+    assert stats["min"] == 36
+    assert stats["max"] == 36
+
+
+def test_copc_mean_z(session, point_cloud_files):
+    """Z values pass through COPC binning"""
+    tools = Tools(session=session)
+    tools.g_region(n=18, s=0, e=18, w=0, res=6)
+    tools.r_in_pdal(input=point_cloud_files["copc"], output="copc_mean", method="mean")
+    stats = tools.r_univar(map="copc_mean", flags="g").keyval
+    assert stats["min"] == pytest.approx(6)
+    assert stats["max"] == pytest.approx(30)
+    assert stats["sum"] == pytest.approx(162)
+
+
+def test_copc_subregion_matches_las(session, point_cloud_files):
+    """Import into a region smaller than the point cloud extent.
+
+    The COPC reader prunes by the region bounds, so the result must be
+    identical to importing the same points from a plain LAS file.
+    """
+    tools = Tools(session=session)
+    tools.g_region(n=12, s=0, e=12, w=0, res=6)
+    tools.r_in_pdal(input=point_cloud_files["copc"], output="copc_sub", method="n")
+    tools.r_in_pdal(input=point_cloud_files["las"], output="las_sub", method="n")
+    tools.r_mapcalc(expression="sub_diff = abs(copc_sub - las_sub)")
+    stats = tools.r_univar(map="sub_diff", flags="g").keyval
+    assert stats["n"] == 4
+    assert stats["max"] == 0
+    stats = tools.r_univar(map="copc_sub", flags="g").keyval
+    assert stats["sum"] == 144
+
+
+def test_copc_full_extent(session, point_cloud_files):
+    """The -e flag imports every point of a COPC file.
+
+    The region is the extent of the data, so the reader is not given
+    bounds to prune by and nothing may be dropped at the edges.
+    """
+    tools = Tools(session=session)
+    tools.r_in_pdal(
+        input=point_cloud_files["copc"],
+        output="copc_full",
+        method="n",
+        resolution=6,
+        flags="e",
+    )
+    tools.g_region(raster="copc_full")
+    stats = tools.r_univar(map="copc_full", flags="g").keyval
+    assert stats["sum"] == 324
+
+
+def test_point_table_capacity(session, point_cloud_files):
+    """Point table capacity does not change the result"""
+    tools = Tools(session=session)
+    tools.g_region(n=18, s=0, e=18, w=0, res=6)
+    tools.r_in_pdal(input=point_cloud_files["las"], output="cap_default", method="n")
+    tools.r_in_pdal(
+        input=point_cloud_files["las"],
+        output="cap_one",
+        method="n",
+        point_table_capacity=1,
+    )
+    tools.r_mapcalc(expression="cap_diff = abs(cap_default - cap_one)")
+    stats = tools.r_univar(map="cap_diff", flags="g").keyval
+    assert stats["max"] == 0
+
+
+def test_point_table_capacity_rejects_zero(session, point_cloud_files):
+    """Zero capacity is rejected by the parser"""
+    tools = Tools(session=session)
+    tools.g_region(n=18, s=0, e=18, w=0, res=6)
+    with pytest.raises(CalledModuleError):
+        tools.r_in_pdal(
+            input=point_cloud_files["las"],
+            output="cap_zero",
+            method="n",
+            point_table_capacity=0,
+        )