diff --git a/raster/r.in.pdal/grassrasterwriter.h b/raster/r.in.pdal/grassrasterwriter.h index 4588cda7988..add2f332bb3 100644 --- a/raster/r.in.pdal/grassrasterwriter.h +++ b/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_; diff --git a/raster/r.in.pdal/info.cpp b/raster/r.in.pdal/info.cpp index e6e62af587b..54d058dc39c 100644 --- a/raster/r.in.pdal/info.cpp +++ b/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"; diff --git a/raster/r.in.pdal/main.cpp b/raster/r.in.pdal/main.cpp index fe17e80cb08..de4ef5d00c9 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 @@ -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("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); 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.