feat: add transform_crs() to reproject point clouds

NAMESPACE

@@ -90,6 +90,7 @@ export(templaz)
 export(temppcd)
 export(tempshp)
 export(temptif)
+export(transform_crs)
 export(transform_with)
 export(triangulate)
 export(unset_exec_option)

NEWS.md

@@ -10,6 +10,7 @@
 - Change: `drop_noise()` and `keep_noise()` now filter classes 7 and 18 instead of 18 only (#283)
 - New: record `use_attribute` as a field in `local_maximum()` output (#310)
 - New: operator `=` in `transform_with()` to assign a field (#314)
+- New: `transform_crs()` reprojects the point cloud to a target CRS (EPSG code or WKT). Unlike `set_crs()`, which only assigns a CRS, `transform_crs()` reprojects the X/Y coordinates with PROJ/GDAL, adapts the scale factors and offsets to the target CRS, updates the bounding box and tags downstream stages and writers with the target CRS.
 
 # lasR 0.20.1
 

R/stages.R

@@ -1130,6 +1130,45 @@ set_crs = function(x)
   stop("Invalid argument")
 }
 
+#' Reproject the point cloud
+#'
+#' Reproject (transform) the point cloud from its current CRS to a target CRS. Unlike
+#' \link{set_crs}, which only **assigns** a CRS without modifying the coordinates, this stage
+#' **reprojects** every point using PROJ/GDAL. The horizontal X and Y coordinates are
+#' reprojected, the scale factors and offsets are adapted to the target CRS (e.g. switching
+#' between metres and degrees), the bounding box is updated and the target CRS is assigned to
+#' the data and to all subsequent stages of the pipeline.
+#'
+#' Only the horizontal coordinates are reprojected: **Z (elevation) is preserved unchanged**,
+#' matching the behaviour of `gdaltransform`, `sf` and `terra` for 2D targets. Vertical datum
+#' transformations (compound/vertical CRS) are out of scope.
+#'
+#' The source CRS is the CRS carried by the files being read or the one assigned with
+#' \link{set_crs} earlier in the pipeline. The stage therefore typically appears after the
+#' \link{reader}. A point that falls outside the domain of the transformation is dropped.
+#'
+#' @param crs integer or string. EPSG code or WKT string of the **target** CRS, understood by GDAL.
+#' @export
+#' @md
+#' @examples
+#' \dontrun{
+#' f = system.file("extdata", "Topography.las", package="lasR")
+#'
+#' # Reproject from the file CRS to WGS 84 (lon/lat) and write the result
+#' pipeline = reader_las() + transform_crs(4326) + write_las()
+#' exec(pipeline, on = f)
+#'
+#' # Reproject then rasterize: the output raster is in the target CRS
+#' pipeline = reader_las() + transform_crs(32619) + rasterize(10, "max")
+#' exec(pipeline, on = f)
+#' }
+transform_crs = function(crs)
+{
+  if (is.numeric(crs)) { return(.APISTAGES$transform_crs_epsg(crs)) }
+  if (is.character(crs)) { return(.APISTAGES$transform_crs_wkt(crs)) }
+  stop("Invalid argument")
+}
+
 #' Sample the point cloud
 #'
 #' Sample the point cloud, keeping one random point per pixel or per voxel or perform a poisson sampling.

python/CMakeLists.txt

@@ -248,6 +248,7 @@ pybind11_add_module(pylasr
   ${LASR_SOURCE_DIR}/src/LASRstages/spikefree.cpp
   ${LASR_SOURCE_DIR}/src/LASRstages/svd.cpp
   ${LASR_SOURCE_DIR}/src/LASRstages/transformwith.cpp
+  ${LASR_SOURCE_DIR}/src/LASRstages/transformcrs.cpp
   ${LASR_SOURCE_DIR}/src/LASRstages/triangulate.cpp
   ${LASR_SOURCE_DIR}/src/LASRstages/writelas.cpp
   ${LASR_SOURCE_DIR}/src/LASRstages/writelax.cpp

python/bindings.cpp

@@ -540,6 +540,20 @@ PYBIND11_MODULE(pylasr, m) {
       return api::Pipeline(s);
       }, "Set CRS using EPSG code or WKT string", py::arg("crs"));
 
+      m.def("transform_crs", [](py::object crs) -> api::Pipeline {
+      api::Stage s("transform_crs");
+      if (py::isinstance<py::int_>(crs)) {
+            s.set("epsg", crs.cast<int>());
+            s.set("wkt", "");
+      } else if (py::isinstance<py::str>(crs)) {
+            s.set("epsg", 0);
+            s.set("wkt", crs.cast<std::string>());
+      } else {
+            throw std::invalid_argument("crs must be either EPSG code (int) or WKT string");
+      }
+      return api::Pipeline(s);
+      }, "Reproject the point cloud to the target CRS (EPSG code or WKT string)", py::arg("crs"));
+
     // Stop conditions
     m.def("stop_if_outside", &api::stop_if_outside,
           "Stop processing if outside bounds",

src/LASRapi/api.cpp

@@ -429,6 +429,24 @@ Pipeline set_crs_wkt(std::string wkt)
   return Pipeline(s);
 }
 
+Pipeline transform_crs_epsg(int epsg)
+{
+  Stage s("transform_crs");
+  s.set("epsg", epsg);
+  s.set("wkt", "");
+
+  return Pipeline(s);
+}
+
+Pipeline transform_crs_wkt(std::string wkt)
+{
+  Stage s("transform_crs");
+  s.set("epsg", 0);
+  s.set("wkt", wkt);
+
+  return Pipeline(s);
+}
+
 Pipeline sampling_voxel(double res, std::vector<std::string> filter, std::string method, int shuffle_size)
 {
   static const std::vector<std::string> choices = {"random"};

src/LASRapi/api.h

@@ -220,6 +220,8 @@ Pipeline remove_rgb();
 Pipeline keep_attributes(std::vector<std::string> names);
 Pipeline set_crs_epsg(int epsg);
 Pipeline set_crs_wkt(std::string wkt);
+Pipeline transform_crs_epsg(int epsg);
+Pipeline transform_crs_wkt(std::string wkt);
 Pipeline sampling_voxel(double res = 2, std::vector<std::string> filter = {""}, std::string method = "random", int shuffle_size = std::numeric_limits<int>::max());
 Pipeline sampling_pixel(double res = 2,  std::vector<std::string> filter = {""}, std::string method = "random", std::string use_attribute = "Z", int shuffle_size = std::numeric_limits<int>::max());
 Pipeline sampling_poisson(double distance = 2,  std::vector<std::string> filter = {""}, int shuffle_size = std::numeric_limits<int>::max());

src/LASRcore/CRS.cpp

@@ -1,7 +1,12 @@
 #include "CRS.h"
 #include "print.h"
 
+#include <ogr_spatialref.h>
+
 #include <stdio.h>
+#include <algorithm>
+#include <limits>
+#include <vector>
 
 CRS::CRS()
 {
@@ -96,6 +101,11 @@ bool CRS::is_feets() const
   return std::fabs(value - 0.3048) < 1e-4;
 }
 
+bool CRS::is_geographic() const
+{
+  return valid && oSRS.IsGeographic();
+}
+
 bool CRS::operator==(const CRS& other) const
 {
   return epsg == other.epsg && valid == other.valid && wkt == other.wkt;
@@ -127,3 +137,70 @@ void CRS::dump() const
 
 // # nocov end
 
+bool reproject_bbox(const CRS& source, const CRS& target, double& xmin, double& ymin, double& xmax, double& ymax)
+{
+  if (!source.is_valid() || !target.is_valid()) return false;
+
+  // Nothing to do for an empty/unset extent.
+  if (xmin > xmax || ymin > ymax) return true;
+
+  OGRSpatialReference oSourceSRS = source.get_crs();
+  OGRSpatialReference oTargetSRS = target.get_crs();
+
+  // Use traditional GIS axis order (x = lon/easting, y = lat/northing) so coordinates
+  // are not swapped under modern PROJ authority-compliant axis ordering.
+  oSourceSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+  oTargetSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+
+  CPLPushErrorHandler(CPLQuietErrorHandler);
+  OGRCoordinateTransformation* ct = OGRCreateCoordinateTransformation(&oSourceSRS, &oTargetSRS);
+  CPLPopErrorHandler();
+
+  if (ct == nullptr) return false;
+
+  const int N = 8; // number of samples per edge
+  std::vector<double> xs;
+  std::vector<double> ys;
+  xs.reserve(4 * (N + 1));
+  ys.reserve(4 * (N + 1));
+
+  for (int i = 0; i <= N; ++i)
+  {
+    double tx = xmin + (xmax - xmin) * i / N;
+    double ty = ymin + (ymax - ymin) * i / N;
+
+    xs.push_back(tx);   ys.push_back(ymin); // bottom edge
+    xs.push_back(tx);   ys.push_back(ymax); // top edge
+    xs.push_back(xmin); ys.push_back(ty);   // left edge
+    xs.push_back(xmax); ys.push_back(ty);   // right edge
+  }
+
+  std::vector<int> ok(xs.size(), 0);
+  ct->Transform((int)xs.size(), xs.data(), ys.data(), nullptr, ok.data());
+  OGRCoordinateTransformation::DestroyCT(ct);
+
+  double nxmin = std::numeric_limits<double>::max();
+  double nymin = std::numeric_limits<double>::max();
+  double nxmax = std::numeric_limits<double>::lowest();
+  double nymax = std::numeric_limits<double>::lowest();
+  bool any = false;
+
+  for (size_t i = 0; i < xs.size(); ++i)
+  {
+    if (!ok[i]) continue;
+    any = true;
+    nxmin = std::min(nxmin, xs[i]);
+    nymin = std::min(nymin, ys[i]);
+    nxmax = std::max(nxmax, xs[i]);
+    nymax = std::max(nymax, ys[i]);
+  }
+
+  if (!any) return false;
+
+  xmin = nxmin;
+  ymin = nymin;
+  xmax = nxmax;
+  ymax = nymax;
+  return true;
+}
+

src/LASRcore/CRS.h

@@ -17,6 +17,7 @@ class CRS
   double get_linear_units() const;
   bool is_meters() const;
   bool is_feets() const;
+  bool is_geographic() const;
   void dump() const;
   bool operator==(const CRS& other) const;
   std::string get_wkt() const;
@@ -28,4 +29,11 @@ class CRS
   OGRSpatialReference oSRS;
 };
 
-#endif
+// Reproject an axis-aligned bounding box from `source` to `target` CRS. The boundary is
+// densified before transforming so a curved (non-affine) reprojection is bounded correctly
+// rather than only at the four corners. Returns false if the CRS are invalid, the
+// transformation cannot be built, or no sample reprojects (the box is entirely outside the
+// transform domain). An empty/unset box (min > max) is left unchanged and returns true.
+bool reproject_bbox(const CRS& source, const CRS& target, double& xmin, double& ymin, double& xmax, double& ymax);
+
+#endif

src/LASRcore/Engine.cpp

@@ -454,15 +454,20 @@ bool Engine::need_points() const
 
 double Engine::need_buffer()
 {
-  for (auto&& stage : pipeline)
+  double required = 0;
+
+  for (auto it = pipeline.rbegin(); it != pipeline.rend(); ++it)
   {
+    Stage* stage = it->get();
     if (!stage->is_streamable())
     {
-      buffer = MAX(buffer, stage->need_buffer());
+      required = MAX(required, stage->need_buffer());
     }
+
+    required = stage->translate_buffer_to_input(required);
   }
 
-  return buffer;
+  return MAX(buffer, required);
 }
 
 void Engine::sort()
@@ -549,4 +554,4 @@ nlohmann::json Engine::to_json()
   }
   print("-------------\n");
   // # nocov end
-}*/
+}*/

src/LASRcore/Stage.h

@@ -94,6 +94,7 @@ class Stage
   virtual double need_buffer() const { return 0; };
   virtual bool need_points() const { return true; };
   virtual void get_extent(double& xmin, double& ymin, double& xmax, double& ymax) { return; };
+  virtual double translate_buffer_to_input(double downstream_buffer) const { return downstream_buffer; };
 
   virtual bool connect(const std::list<std::unique_ptr<Stage>>&, const std::string& uid) { return true; };
 

src/LASRcore/parser.cpp

@@ -36,6 +36,7 @@
 #include "svd.h"
 #include "triangulate.h"
 #include "transformwith.h"
+#include "transformcrs.h"
 #include "writelas.h"
 #include "writelax.h"
 #include "writevpc.h"
@@ -85,6 +86,11 @@ bool Engine::parse(const nlohmann::json& json, bool progress)
   bool reader = false;
   bool indexer = false;
 
+  // Tracks the CRS as it flows through the stages while parsing. It is needed so that
+  // stages whose get_extent() depends on the CRS (e.g. transform_crs, which reprojects
+  // the coverage extent) can be resolved before the extent is propagated downstream.
+  CRS parse_crs;
+
   parsed = false;
   pipeline.clear();
 
@@ -124,6 +130,7 @@ bool Engine::parse(const nlohmann::json& json, bool progress)
     {"summarise",            create_instance<LASRsummary>},
     {"svd",                  create_instance<LASRsvd>},
     {"transform_with",       create_instance<LASRtransformwith>},
+    {"transform_crs",        create_instance<LASRtransformcrs>},
     {"triangulate",          create_instance<LASRtriangulate>},
     {"write_las",            create_instance<LASRlaswriter>},
     {"write_vpc",            create_instance<LASRvpcwriter>},
@@ -249,6 +256,9 @@ bool Engine::parse(const nlohmann::json& json, bool progress)
           last_error = "Internal error: bad build_catalog stage";
           return false;
         }
+
+        // The catalog now carries the source CRS of the whole pipeline.
+        if (catalog != nullptr) parse_crs = catalog->get_crs();
       }
       else if (name == "reader")
       {
@@ -456,6 +466,12 @@ bool Engine::parse(const nlohmann::json& json, bool progress)
           return false;
         }
 
+        // Propagate the CRS so stages that reproject (transform_crs) know their source
+        // CRS before their extent is computed and pushed to the next stages. The final
+        // CRS assignment (used for outputs) is redone in the second pass below.
+        it->set_crs(parse_crs);
+        parse_crs = it->get_crs();
+
         it->get_extent(xmin, ymin, xmax, ymax);
 
         // If we intend to actually process the point cloud we check that a reader stage is present if needed