r.in.pdal: Add COPC support

raster/r.in.pdal/grassrasterwriter.h

@@ -39,7 +39,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"; }
 
@@ -96,17 +103,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;
         }
 
@@ -117,6 +118,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

@@ -140,6 +140,13 @@ 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();
+        pdal::SpatialReference spatial_reference = table.spatialReference();
+        /* The point table SRS may not be populated before execute();
+         * read it from the reader in that case. */
+        if (spatial_reference.empty()) {
+            spatial_reference = las_reader.getSpatialReference();
+        }
+        std::string proj_wkt = spatial_reference.getWKT();
 
         std::cout << "File: " << infile << std::endl;
         std::cout << "File version = "
@@ -169,6 +176,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

@@ -18,6 +18,8 @@
  *****************************************************************************/
 
 #include <cstdio>
+#include <iomanip>
+#include <sstream>
 
 #if defined(__clang__)
 #pragma clang diagnostic push
@@ -371,6 +373,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';
@@ -736,6 +752,15 @@ int main(int argc, char *argv[])
             las_opts.add(nosrs_opt);
         }
 #endif
+        // Let the COPC octree skip nodes outside the region; bounds are in
+        // the file's CRS, the GRASS filter still does the exact clip.
+        if (use_spatial_filter && !reproject_flag->answer &&
+            pdal_read_driver == "readers.copc") {
+            std::ostringstream bounds_str;
+            bounds_str << std::setprecision(17) << "([" << xmin << ", " << xmax
+                       << "], [" << ymin << ", " << ymax << "])";
+            las_opts.add(pdal::Option("bounds", bounds_str.str()));
+        }
         // stages created by factory are destroyed with the factory
         pdal::Stage *reader = factory.createStage(pdal_read_driver);
         if (!reader)
@@ -796,9 +821,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);
@@ -815,7 +841,17 @@ int main(int argc, char *argv[])
     }
     else if (!reproject_flag->answer) {
         pdal::SpatialReference spatial_reference =
-            merge_filter.getSpatialReference();
+            point_table.spatialReference();
+        /* COPC and some other readers only populate the point table SRS after
+         * execute(); fall back to reading it from the reader stages directly.
+         */
+        if (spatial_reference.empty()) {
+            for (pdal::Stage *reader : readers) {
+                spatial_reference = reader->getSpatialReference();
+                if (!spatial_reference.empty())
+                    break;
+            }
+        }
         if (spatial_reference.empty())
             G_fatal_error(_("The input dataset has undefined projection"));
         std::string dataset_wkt = spatial_reference.getWKT();
@@ -961,6 +997,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,9 +359,14 @@ <h3>Format and projection support</h3>
 
 The typical file extensions for the LAS format are .las and .laz
 (compressed). The compressed LAS (.laz) format can be imported only if
-libLAS has been compiled with LASzip support. It is also recommended to
-compile libLAS with GDAL which is used to test if the LAS coordinate reference
-system matches that of the GRASS project (previously called location).
+PDAL has been compiled with LASzip support. 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. The coordinate reference system of the input is
+expected to match that of the GRASS project; use the <b>-o</b> flag to
+override the projection check or the <b>-w</b> flag to reproject the
+points during import.
 
 <!--
 <h3>LAS file import preparations</h3>
@@ -421,6 +426,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,10 +324,14 @@ parameter is provided.
 
 The typical file extensions for the LAS format are .las and .laz
 (compressed). The compressed LAS (.laz) format can be imported only if
-libLAS has been compiled with LASzip support. It is also recommended to
-compile libLAS with GDAL which is used to test if the LAS coordinate
-reference system matches that of the GRASS project (previously called
-location).
+PDAL has been compiled with LASzip support. 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. The coordinate reference system of the input is
+expected to match that of the GRASS project; use the **-o** flag to
+override the projection check or the **-w** flag to reproject the
+points during import.
 
 ### Memory consumption
 
@@ -356,6 +360,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,75 @@
+"""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_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,
+        )