From 7732018c633c81bded3c5fb6edef1803f9f3250a Mon Sep 17 00:00:00 2001 From: cwhite911 Date: Tue, 24 Mar 2026 13:43:10 -0400 Subject: [PATCH 1/3] r.in.pdal: Added ability to read COPC data --- raster/r.in.pdal/grassrasterwriter.h | 12 ++++++--- raster/r.in.pdal/info.cpp | 12 +++++++++ raster/r.in.pdal/main.cpp | 37 +++++++++++++++++++++++++--- 3 files changed, 54 insertions(+), 7 deletions(-) diff --git a/raster/r.in.pdal/grassrasterwriter.h b/raster/r.in.pdal/grassrasterwriter.h index 4588cda7988..4e0155a85e0 100644 --- a/raster/r.in.pdal/grassrasterwriter.h +++ b/raster/r.in.pdal/grassrasterwriter.h @@ -39,7 +39,13 @@ class GrassRasterWriter : public pdal::NoFilenameWriter, class GrassRasterWriter : public pdal::Writer, public pdal::Streamable { #endif public: - GrassRasterWriter() : n_processed(0) {} + GrassRasterWriter() + : n_processed(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"; } @@ -105,8 +111,8 @@ class GrassRasterWriter : public pdal::Writer, public pdal::Streamable { 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.")); + G_debug(3, "A point on the edge of computational region detected. " + "Ignoring."); return false; } diff --git a/raster/r.in.pdal/info.cpp b/raster/r.in.pdal/info.cpp index e6e62af587b..0687f1374b5 100644 --- a/raster/r.in.pdal/info.cpp +++ b/raster/r.in.pdal/info.cpp @@ -140,6 +140,14 @@ 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(); + /* 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()) { + spatial_reference = las_reader.getSpatialReference(); + } + std::string proj_wkt = spatial_reference.getWKT(); std::cout << "File: " << infile << std::endl; std::cout << "File version = " @@ -153,6 +161,10 @@ void print_lasinfo(struct StringList *infiles) std::cout << "Creation DOY: " << h.creationDOY() << "\n"; std::cout << "Creation Year: " << h.creationYear() << "\n"; std::cout << "VLR offset (header size): " << h.vlrOffset() << "\n"; + if (!proj_wkt.empty()) + std::cout << "Projection (WKT): " << proj_wkt << "\n"; + else + std::cout << "Projection: (undefined)" << "\n"; std::cout << "VLR Count: " << h.vlrCount() << "\n"; std::cout << "Point format: " << (int)h.pointFormat() << "\n"; std::cout << "Point offset: " << h.pointOffset() << "\n"; diff --git a/raster/r.in.pdal/main.cpp b/raster/r.in.pdal/main.cpp index fe17e80cb08..7968e3a72ab 100644 --- a/raster/r.in.pdal/main.cpp +++ b/raster/r.in.pdal/main.cpp @@ -371,6 +371,21 @@ 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("10000"); + + point_table_capacity_opt->description = + _("Set capacity of point table used for buffering points during " + "processing. " + "Increasing this value may improve performance for large datasets, " + "but also increases memory usage. Default is 10,000 points, which " + "should be sufficient for most cases."); + point_table_capacity_opt->guisection = _("Performance"); + Flag *extents_flag = G_define_flag(); extents_flag->key = 'e'; @@ -796,9 +811,13 @@ 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; + // there is no difference between 1 and 10k points in memory + // consumption, so using 10k in case it is faster for some cases + // Added user option to set point table capacity to a higher value, which + // may improve performance for large datasets but also increases memory + // usage + 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 +834,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(); From 927a399a85c8225263695b6efcd798d02e3a3188 Mon Sep 17 00:00:00 2001 From: Corey White Date: Fri, 31 Jul 2026 21:05:13 -0400 Subject: [PATCH 2/3] r.in.pdal: Push region bounds to COPC reader When importing into a region smaller than the point cloud, pass the region bounds to readers.copc so its octree index skips nodes outside the region instead of decoding every point. Bounds are evaluated in the file's CRS, which is expected to match the project's CRS, so the pushdown is skipped when reprojection is requested. The GRASS spatial filter is kept for the exact per-point clip. --- raster/r.in.pdal/main.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/raster/r.in.pdal/main.cpp b/raster/r.in.pdal/main.cpp index 7968e3a72ab..b4345ceb6b2 100644 --- a/raster/r.in.pdal/main.cpp +++ b/raster/r.in.pdal/main.cpp @@ -18,6 +18,8 @@ *****************************************************************************/ #include +#include +#include #if defined(__clang__) #pragma clang diagnostic push @@ -751,6 +753,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) From 8fae4dac6881fbab4c3e0aec4708a3282a20c274 Mon Sep 17 00:00:00 2001 From: Corey White Date: Fri, 31 Jul 2026 21:52:33 -0400 Subject: [PATCH 3/3] r.in.pdal: Validate point_table_capacity, report edge points, add COPC tests and docs --- raster/r.in.pdal/grassrasterwriter.h | 14 ++-- raster/r.in.pdal/info.cpp | 13 ++-- raster/r.in.pdal/main.cpp | 22 +++--- raster/r.in.pdal/r.in.pdal.html | 16 +++- raster/r.in.pdal/r.in.pdal.md | 16 +++- raster/r.in.pdal/tests/conftest.py | 67 +++++++++++++++++ raster/r.in.pdal/tests/r_in_pdal_copc_test.py | 75 +++++++++++++++++++ 7 files changed, 189 insertions(+), 34 deletions(-) create mode 100644 raster/r.in.pdal/tests/conftest.py create mode 100644 raster/r.in.pdal/tests/r_in_pdal_copc_test.py diff --git a/raster/r.in.pdal/grassrasterwriter.h b/raster/r.in.pdal/grassrasterwriter.h index 4e0155a85e0..add2f332bb3 100644 --- a/raster/r.in.pdal/grassrasterwriter.h +++ b/raster/r.in.pdal/grassrasterwriter.h @@ -40,8 +40,9 @@ class GrassRasterWriter : public pdal::Writer, public pdal::Streamable { #endif public: GrassRasterWriter() - : n_processed(0), region_(nullptr), point_binning_(nullptr), - bin_index_nodes_(nullptr), rtype_(FCELL_TYPE), cols_(0), scale_(1.0), + : 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) { @@ -102,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_debug(3, "A point on the edge of computational region detected. " - "Ignoring."); + n_on_edge++; return false; } @@ -123,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_; diff --git a/raster/r.in.pdal/info.cpp b/raster/r.in.pdal/info.cpp index 0687f1374b5..54d058dc39c 100644 --- a/raster/r.in.pdal/info.cpp +++ b/raster/r.in.pdal/info.cpp @@ -141,9 +141,8 @@ void print_lasinfo(struct StringList *infiles) pdal::PointLayoutPtr point_layout = table.layout(); const pdal::Dimension::IdList &dims = point_layout->dims(); pdal::SpatialReference spatial_reference = 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. - */ + /* 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(); } @@ -161,10 +160,6 @@ void print_lasinfo(struct StringList *infiles) std::cout << "Creation DOY: " << h.creationDOY() << "\n"; std::cout << "Creation Year: " << h.creationYear() << "\n"; std::cout << "VLR offset (header size): " << h.vlrOffset() << "\n"; - if (!proj_wkt.empty()) - std::cout << "Projection (WKT): " << proj_wkt << "\n"; - else - std::cout << "Projection: (undefined)" << "\n"; std::cout << "VLR Count: " << h.vlrCount() << "\n"; std::cout << "Point format: " << (int)h.pointFormat() << "\n"; std::cout << "Point offset: " << h.pointOffset() << "\n"; @@ -181,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"; diff --git a/raster/r.in.pdal/main.cpp b/raster/r.in.pdal/main.cpp index b4345ceb6b2..de4ef5d00c9 100644 --- a/raster/r.in.pdal/main.cpp +++ b/raster/r.in.pdal/main.cpp @@ -379,13 +379,12 @@ int main(int argc, char *argv[]) point_table_capacity_opt->type = TYPE_INTEGER; point_table_capacity_opt->required = NO; point_table_capacity_opt->answer = const_cast("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 = - _("Set capacity of point table used for buffering points during " - "processing. " - "Increasing this value may improve performance for large datasets, " - "but also increases memory usage. Default is 10,000 points, which " - "should be sufficient for most cases."); + _("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(); @@ -822,11 +821,8 @@ 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 - // Added user option to set point table capacity to a higher value, which - // may improve performance for large datasets but also increases memory - // usage + // 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); @@ -1001,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); diff --git a/raster/r.in.pdal/r.in.pdal.html b/raster/r.in.pdal/r.in.pdal.html index d26ea92d4e8..503871852d7 100644 --- a/raster/r.in.pdal/r.in.pdal.html +++ b/raster/r.in.pdal/r.in.pdal.html @@ -359,9 +359,14 @@

Format and projection support

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.