diff --git a/raster/r.in.pdal/grassrasterwriter.h b/raster/r.in.pdal/grassrasterwriter.h index caf6d46e0ce..a8d091ff936 100644 --- a/raster/r.in.pdal/grassrasterwriter.h +++ b/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_; diff --git a/raster/r.in.pdal/info.cpp b/raster/r.in.pdal/info.cpp index 00fad0d9cfd..99f61fa27fa 100644 --- a/raster/r.in.pdal/info.cpp +++ b/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"; diff --git a/raster/r.in.pdal/main.cpp b/raster/r.in.pdal/main.cpp index 482346fc98c..2e88ff55f29 100644 --- a/raster/r.in.pdal/main.cpp +++ b/raster/r.in.pdal/main.cpp @@ -16,6 +16,8 @@ *****************************************************************************/ #include +#include +#include #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("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); diff --git a/raster/r.in.pdal/r.in.pdal.html b/raster/r.in.pdal/r.in.pdal.html index bf2f10727b2..7b4b86016f9 100644 --- a/raster/r.in.pdal/r.in.pdal.html +++ b/raster/r.in.pdal/r.in.pdal.html @@ -359,7 +359,13 @@

Format and projection support

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). +.laz (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 @@ -438,6 +444,11 @@

Memory consumption

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. diff --git a/raster/r.in.pdal/r.in.pdal.md b/raster/r.in.pdal/r.in.pdal.md index 1684fd2d09c..2c6a9f81fec 100644 --- a/raster/r.in.pdal/r.in.pdal.md +++ b/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. diff --git a/raster/r.in.pdal/tests/conftest.py b/raster/r.in.pdal/tests/conftest.py new file mode 100644 index 00000000000..19cafc6c259 --- /dev/null +++ b/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 diff --git a/raster/r.in.pdal/tests/r_in_pdal_copc_test.py b/raster/r.in.pdal/tests/r_in_pdal_copc_test.py new file mode 100644 index 00000000000..2e0d5b42f71 --- /dev/null +++ b/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, + )