diff --git a/raster/r.in.pdal/info.cpp b/raster/r.in.pdal/info.cpp index e6e62af587b..6787aa02094 100644 --- a/raster/r.in.pdal/info.cpp +++ b/raster/r.in.pdal/info.cpp @@ -8,11 +8,18 @@ * Read the COPYING file that comes with GRASS for details. * */ +#include + +#include +#include #include "info.h" -#include -#ifdef PDAL_USE_NOSRS +extern "C" { +#include "projection.h" +} + +#ifdef R_IN_PDAL_USE_NOSRS void get_extent(struct StringList *infiles, double *min_x, double *max_x, double *min_y, double *max_y, double *min_z, double *max_z, bool nosrs) @@ -22,8 +29,9 @@ void get_extent(struct StringList *infiles, double *min_x, double *max_x, #endif { pdal::StageFactory factory; - bool first = 1; - + pdal::SpatialReference spatial_reference; + bool first = true; + double min_x_, max_x_, min_y_, max_y_, min_z_, max_z_; *min_x = *max_x = *min_y = *max_y = *min_z = *max_z = NAN; for (int i = 0; i < infiles->num_items; i++) { @@ -38,7 +46,7 @@ void get_extent(struct StringList *infiles, double *min_x, double *max_x, pdal::Options las_opts; pdal::Option las_opt("filename", infile); las_opts.add(las_opt); -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS if (nosrs) { pdal::Option nosrs_opt("nosrs", true); las_opts.add(nosrs_opt); @@ -52,35 +60,162 @@ void get_extent(struct StringList *infiles, double *min_x, double *max_x, catch (const std::exception &err) { G_fatal_error(_("PDAL error: %s"), err.what()); } + spatial_reference = las_reader.getSpatialReference(); + std::string dataset_wkt = spatial_reference.getWKT(); + bool proj_match = is_wkt_projection_same_as_loc(dataset_wkt.c_str()); + const pdal::LasHeader &las_header = las_reader.header(); + + min_x_ = las_header.minX(); + min_y_ = las_header.minY(); + min_z_ = las_header.minZ(); + max_x_ = las_header.maxX(); + max_y_ = las_header.maxY(); + max_z_ = las_header.maxZ(); + bool srs_known = + spatial_reference.valid() && !spatial_reference.empty(); +#ifdef R_IN_PDAL_USE_NOSRS + bool check_srs = !nosrs; +#else + bool check_srs = true; +#endif + bool need_to_reproject = check_srs && !proj_match && srs_known; + + // The extent of a file without a CRS cannot be reprojected, so its + // header values are reported as they are. Unlike the import, which + // fails in this case, printing the extent continues, because the + // user may still need the numbers to set up the project. + if (check_srs && !srs_known) + G_warning(_("The input dataset <%s> has undefined CRS, " + "its extent is reported as stored in the file"), + infile); + + if (need_to_reproject) { + get_reprojected_extent(spatial_reference, &min_x_, &max_x_, &min_y_, + &max_y_, &min_z_, &max_z_); + } if (first) { - *min_x = las_header.minX(); - *min_y = las_header.minY(); - *min_z = las_header.minZ(); - *max_x = las_header.maxX(); - *max_y = las_header.maxY(); - *max_z = las_header.maxZ(); - - first = 0; + *min_x = min_x_; + *min_y = min_y_; + *min_z = min_z_; + *max_x = max_x_; + *max_y = max_y_; + *max_z = max_z_; + + first = false; + } + else { + if (*min_x > min_x_) + *min_x = min_x_; + if (*min_y > min_y_) + *min_y = min_y_; + if (*min_z > min_z_) + *min_z = min_z_; + if (*max_x < max_x_) + *max_x = max_x_; + if (*max_y < max_y_) + *max_y = max_y_; + if (*max_z < max_z_) + *max_z = max_z_; + } + } +} +/* Reproject bounding box from the input CRS to the project CRS using PDAL. + * Points are sampled along bbox edges (not just corners) to account for + * reprojection curvature. */ +void get_reprojected_extent(pdal::SpatialReference &spatial_reference, + double *min_x, double *max_x, double *min_y, + double *max_y, double *min_z, double *max_z) +{ + pdal::PointTable table; + table.layout()->registerDim(pdal::Dimension::Id::X); + table.layout()->registerDim(pdal::Dimension::Id::Y); + table.layout()->registerDim(pdal::Dimension::Id::Z); + + pdal::PointViewPtr view(new pdal::PointView(table)); + pdal::PointId idx = 0; + int points_per_edge = 10; + // Generate points along the 4 edges of the bounding box + // First point of each edge starts at the corner and + // the last point is right before the corner + for (int i = 0; i < points_per_edge - 1; ++i) { + double t = i / (double)(points_per_edge - 1); + + // South edge (y = min_y) + view->setField(pdal::Dimension::Id::X, idx, + *min_x + t * (*max_x - *min_x)); + view->setField(pdal::Dimension::Id::Y, idx, *min_y); + view->setField(pdal::Dimension::Id::Z, idx, *min_z); + idx++; + + // East edge (x = max_x) + view->setField(pdal::Dimension::Id::X, idx, *max_x); + view->setField(pdal::Dimension::Id::Y, idx, + *min_y + t * (*max_y - *min_y)); + view->setField(pdal::Dimension::Id::Z, idx, *max_z); + idx++; + + // North edge (y = max_y) + view->setField(pdal::Dimension::Id::X, idx, + *max_x - t * (*max_x - *min_x)); + view->setField(pdal::Dimension::Id::Y, idx, *max_y); + view->setField(pdal::Dimension::Id::Z, idx, *min_z); + idx++; + + // West edge (x = min_x) + view->setField(pdal::Dimension::Id::X, idx, *min_x); + view->setField(pdal::Dimension::Id::Y, idx, + *max_y - t * (*max_y - *min_y)); + view->setField(pdal::Dimension::Id::Z, idx, *max_z); + idx++; + } + + pdal::BufferReader reader; + reader.addView(view); + + pdal::StageFactory factory; + pdal::Stage *reproject = factory.createStage("filters.reprojection"); + pdal::Options reproject_options; + reproject_options.add("in_srs", spatial_reference.getWKT()); + reproject_options.add("out_srs", location_projection_as_wkt(false)); + reproject->setOptions(reproject_options); + reproject->setInput(reader); + try { + reproject->prepare(table); + reproject->execute(table); + } + catch (const std::exception &err) { + G_fatal_error(_("Reprojection of the data extent failed: %s"), + err.what()); + } + + for (pdal::PointId i = 0; i < view->size(); ++i) { + double x = view->getFieldAs(pdal::Dimension::Id::X, i); + double y = view->getFieldAs(pdal::Dimension::Id::Y, i); + double z = view->getFieldAs(pdal::Dimension::Id::Z, i); + if (i == 0) { + *min_x = *max_x = x; + *min_y = *max_y = y; + *min_z = *max_z = z; } else { - if (*min_x > las_header.minX()) - *min_x = las_header.minX(); - if (*min_y > las_header.minY()) - *min_y = las_header.minY(); - if (*min_z > las_header.minZ()) - *min_z = las_header.minZ(); - if (*max_x < las_header.maxX()) - *max_x = las_header.maxX(); - if (*max_y < las_header.maxY()) - *max_y = las_header.maxY(); - if (*max_z < las_header.maxZ()) - *max_z = las_header.maxZ(); + if (*min_x > x) + *min_x = x; + if (*min_y > y) + *min_y = y; + if (*min_z > z) + *min_z = z; + if (*max_x < x) + *max_x = x; + if (*max_y < y) + *max_y = y; + if (*max_z < z) + *max_z = z; } } } -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS void print_extent(struct StringList *infiles, bool nosrs) #else void print_extent(struct StringList *infiles) @@ -88,7 +223,7 @@ void print_extent(struct StringList *infiles) { double min_x, max_x, min_y, max_y, min_z, max_z; -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS get_extent(infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z, nosrs); #else get_extent(infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z); @@ -97,7 +232,7 @@ void print_extent(struct StringList *infiles) min_x, min_z, max_z); } -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS void print_lasinfo(struct StringList *infiles, bool nosrs) #else void print_lasinfo(struct StringList *infiles) @@ -123,7 +258,7 @@ void print_lasinfo(struct StringList *infiles) pdal::Options las_opts; pdal::Option las_opt("filename", infile); las_opts.add(las_opt); -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS if (nosrs) { pdal::Option nosrs_opt("nosrs", true); las_opts.add(nosrs_opt); diff --git a/raster/r.in.pdal/info.h b/raster/r.in.pdal/info.h index 09736bab543..47101fd048a 100644 --- a/raster/r.in.pdal/info.h +++ b/raster/r.in.pdal/info.h @@ -23,6 +23,7 @@ #include #include #include +#include #if defined(__clang__) #pragma clang diagnostic pop #endif @@ -31,7 +32,7 @@ #if (PDAL_VERSION_MAJOR >= 2 && PDAL_VERSION_MINOR > 4) || \ (PDAL_VERSION_MAJOR == 2 && PDAL_VERSION_MINOR == 4 && \ PDAL_VERSION_PATCH == 3) -#define PDAL_USE_NOSRS 1 +#define R_IN_PDAL_USE_NOSRS 1 #endif extern "C" { @@ -40,7 +41,7 @@ extern "C" { #include "string_list.h" } -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS void get_extent(struct StringList *, double *, double *, double *, double *, double *, double *, bool); void print_extent(struct StringList *, bool); @@ -51,5 +52,7 @@ void get_extent(struct StringList *, double *, double *, double *, double *, void print_extent(struct StringList *); void print_lasinfo(struct StringList *); #endif +void get_reprojected_extent(pdal::SpatialReference &, double *, double *, + double *, double *, double *, double *); #endif // INFO_H diff --git a/raster/r.in.pdal/main.cpp b/raster/r.in.pdal/main.cpp index fe17e80cb08..32231994cda 100644 --- a/raster/r.in.pdal/main.cpp +++ b/raster/r.in.pdal/main.cpp @@ -250,11 +250,10 @@ int main(int argc, char *argv[]) reproject_flag->key = 'w'; reproject_flag->label = - _("Reproject to project's coordinate system if needed"); + _("Reproject to project's coordinate system if needed [deprecated]"); reproject_flag->description = - _("Reprojects input dataset to the coordinate system of" - " the GRASS project (by default only datasets with" - " matching coordinate system can be imported"); + _("This flag is deprecated and will be removed in a future release. " + "Input dataset will always be reprojected if needed."); reproject_flag->guisection = _("Projection"); // TODO: from the API it seems that also prj file path and proj string will @@ -446,7 +445,7 @@ int main(int argc, char *argv[]) /* If we print extent, there is no need to validate rest of the input */ if (print_extent_flag->answer) { -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS print_extent(&infiles, over_flag->answer); #else print_extent(&infiles); @@ -455,7 +454,7 @@ int main(int argc, char *argv[]) } if (print_info_flag->answer) { -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS print_lasinfo(&infiles, over_flag->answer); #else print_lasinfo(&infiles); @@ -463,6 +462,11 @@ int main(int argc, char *argv[]) exit(EXIT_SUCCESS); } + if (reproject_flag->answer) { + G_verbose_message( + _("Flag 'w' is deprecated and will be removed in a future release. " + "Input dataset will always be reprojected if needed.")); + } /* we could use rules but this gives more info and allows continuing */ if (set_region_flag->answer && !(extents_flag->answer || res_opt->answer || base_rast_res_flag->answer)) { @@ -515,7 +519,7 @@ int main(int argc, char *argv[]) if (extents_flag->answer) { double min_x, max_x, min_y, max_y, min_z, max_z; -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS get_extent(&infiles, &min_x, &max_x, &min_y, &max_y, &min_z, &max_z, over_flag->answer); #else @@ -718,6 +722,7 @@ int main(int argc, char *argv[]) std::vector readers; pdal::StageFactory factory; pdal::MergeFilter merge_filter; + bool need_to_reproject = false; /* loop of input files */ for (int i = 0; i < infiles.num_items; i++) { const char *infile = infiles.items[i]; @@ -730,7 +735,7 @@ int main(int argc, char *argv[]) pdal::Options las_opts; pdal::Option las_opt("filename", infile); las_opts.add(las_opt); -#ifdef PDAL_USE_NOSRS +#ifdef R_IN_PDAL_USE_NOSRS if (over_flag->answer) { pdal::Option nosrs_opt("nosrs", true); las_opts.add(nosrs_opt); @@ -744,8 +749,35 @@ int main(int argc, char *argv[]) infile); reader->setOptions(las_opts); readers.push_back(reader); + + // getting projection is possible only after prepare + if (!over_flag->answer) { + pdal::PointTable table; + try { + reader->prepare(table); + } + catch (const std::exception &err) { + G_fatal_error(_("PDAL error while reading <%s>: %s"), infile, + err.what()); + } + pdal::SpatialReference spatial_reference = + reader->getSpatialReference(); + 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()); + + if (!proj_match) + need_to_reproject = true; + } merge_filter.setInput(*reader); } + if (over_flag->answer) { + G_important_message(_("Overriding projection check and assuming" + " that the CRS of input matches" + " the project's CRS")); + } // we need to keep pointer to the last stage // merge filter puts the n readers into one stage, @@ -753,8 +785,7 @@ int main(int argc, char *argv[]) pdal::Stage *last_stage = &merge_filter; pdal::ReprojectionFilter reprojection_filter; - // we reproject when requested regardless of the input projection - if (reproject_flag->answer) { + if (need_to_reproject) { G_message(_("Reprojecting the input to the project's CRS")); char *proj_wkt = location_projection_as_wkt(false); @@ -807,24 +838,6 @@ int main(int argc, char *argv[]) G_fatal_error(_("PDAL error: %s"), err.what()); } - // getting projection is possible only after prepare - if (over_flag->answer) { - G_important_message(_("Overriding projection check and assuming" - " that the CRS of input matches" - " the project's CRS")); - } - else if (!reproject_flag->answer) { - pdal::SpatialReference spatial_reference = - merge_filter.getSpatialReference(); - 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()); - - if (!proj_match) - wkt_projection_mismatch_report(dataset_wkt.c_str()); - } - G_important_message(_("Running PDAL algorithms...")); // get the layout to see the dimensions diff --git a/raster/r.in.pdal/projection.c b/raster/r.in.pdal/projection.c index a2b17d8d50e..8c177b17f3a 100644 --- a/raster/r.in.pdal/projection.c +++ b/raster/r.in.pdal/projection.c @@ -88,9 +88,6 @@ void projection_mismatch_report(struct Cell_head cellhd, " in the coordinate reference system definitions," " use the -o flag to ignore them and use" " current project definition.\n")); - strcat(error_msg, - _("Consider generating a new project with 'project' parameter" - " from input data set.\n")); G_fatal_error("%s", error_msg); } diff --git a/raster/r.in.pdal/r.in.pdal.html b/raster/r.in.pdal/r.in.pdal.html index d26ea92d4e8..bf2f10727b2 100644 --- a/raster/r.in.pdal/r.in.pdal.html +++ b/raster/r.in.pdal/r.in.pdal.html @@ -357,11 +357,28 @@

NOTES

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). +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). + +

+The coordinate reference system (CRS) of the input is read from the file +metadata and compared with the CRS of the current project (previously called +location). When the two differ, the points are reprojected to the project's +CRS during the import. The -w flag, which used to be required to allow +the reprojection, is deprecated and has no effect. + +

+The -o flag skips the CRS check and assumes that the input is already +in the project's CRS. Use it when the file has no CRS metadata or when the +metadata is known to be wrong; without it, the import of a file without CRS +metadata fails. + +

+The extent printed with the -g flag is reprojected as well, so that it +matches the CRS of the current project and can be passed directly to +g.region. For a file without CRS metadata, the +extent is printed as stored in the file and a warning is issued.