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

raster/r.in.pdal/info.cpp

@@ -8,11 +8,18 @@
  *   Read the COPYING file that comes with GRASS for details.
  *
  */
+#include <cmath>
+
+#include <pdal/filters/ReprojectionFilter.hpp>
+#include <pdal/io/BufferReader.hpp>
 
 #include "info.h"
-#include <cmath>
 
-#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,43 +60,170 @@ 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<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;
         }
     }
 }
 
-#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)
 #endif
 {
     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);

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
@@ -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