r.in.pdal: fix reported data extent when CRS doesn't match

raster/r.in.pdal/info.cpp

@@ -8,9 +8,17 @@
  *   Read the COPYING file that comes with GRASS for details.
  *
  */
+#include <cmath>
+#include <stdbool.h>
+
+#include <pdal/filters/ReprojectionFilter.hpp>
+#include <pdal/io/BufferReader.hpp>
 
 #include "info.h"
-#include <cmath>
+
+extern "C" {
+#include "projection.h"
+}
 
 #ifdef PDAL_USE_NOSRS
 void get_extent(struct StringList *infiles, double *min_x, double *max_x,
@@ -22,8 +30,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++) {
@@ -52,30 +61,118 @@ 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();
+#ifdef PDAL_USE_NOSRS
+        bool need_to_reproject = !nosrs && !proj_match &&
+                                 spatial_reference.valid() &&
+                                 !spatial_reference.empty();
+#else
+        bool need_to_reproject = !proj_match && spatial_reference.valid() &&
+                                 !spatial_reference.empty();
+#endif
+
+        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_;
+        }
+    }
+}
+
+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));
+
+    view->setField(pdal::Dimension::Id::X, 0, *min_x);
+    view->setField(pdal::Dimension::Id::Y, 0, *min_y);
+    view->setField(pdal::Dimension::Id::Z, 0, *min_z);
+    view->setField(pdal::Dimension::Id::X, 1, *min_x);
+    view->setField(pdal::Dimension::Id::Y, 1, *max_y);
+    view->setField(pdal::Dimension::Id::Z, 1, *min_z);
+    view->setField(pdal::Dimension::Id::X, 2, *max_x);
+    view->setField(pdal::Dimension::Id::Y, 2, *min_y);
+    view->setField(pdal::Dimension::Id::Z, 2, *max_z);
+    view->setField(pdal::Dimension::Id::X, 3, *max_x);
+    view->setField(pdal::Dimension::Id::Y, 3, *max_y);
+    view->setField(pdal::Dimension::Id::Z, 3, *max_z);
+
+    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);
+    reproject->prepare(table);
+    reproject->execute(table);
+
+    for (pdal::PointId i = 0; i < view->size(); ++i) {
+
+        double x = view->getFieldAs<double>(pdal::Dimension::Id::X, i);
+        double y = view->getFieldAs<double>(pdal::Dimension::Id::Y, i);
+        double z = view->getFieldAs<double>(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;
         }
     }
 }

raster/r.in.pdal/info.h

@@ -23,6 +23,7 @@
 #include <pdal/Options.hpp>
 #include <pdal/io/LasReader.hpp>
 #include <pdal/io/LasHeader.hpp>
+#include <pdal/SpatialReference.hpp>
 #if defined(__clang__)
 #pragma clang diagnostic pop
 #endif
@@ -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

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
@@ -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)) {
@@ -718,6 +722,8 @@ int main(int argc, char *argv[])
     std::vector<pdal::Stage *> readers;
     pdal::StageFactory factory;
     pdal::MergeFilter merge_filter;
+    pdal::SpatialReference spatial_reference;
+    bool need_to_reproject = false;
     /* loop of input files */
     for (int i = 0; i < infiles.num_items; i++) {
         const char *infile = infiles.items[i];
@@ -744,8 +750,28 @@ 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;
+            reader->prepare(table);
+            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,
@@ -754,7 +780,7 @@ int main(int argc, char *argv[])
     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 +833,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

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);
 }