Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
fab7218
r.in.pdal: fix reported data extent when CRS doesn't match
petrasovaa Oct 7, 2025
3787811
fix tests
petrasovaa Oct 7, 2025
b7e942f
forgot to add test las files
petrasovaa Oct 7, 2025
4eb4596
address review, improve test
petrasovaa Oct 18, 2025
a325585
set higher delta for macOS
petrasovaa Oct 18, 2025
62eed27
add points along boox edges
petrasovaa Nov 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
177 changes: 148 additions & 29 deletions raster/r.in.pdal/info.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,19 @@
* 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>

#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)
Expand All @@ -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++) {
Expand All @@ -38,7 +47,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);
Expand All @@ -52,43 +61,153 @@ 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 R_IN_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_;
}
}
}
/* Reproject 4 points representing bounding box using PDAL pipeline */
Copy link
Contributor

@metzm metzm Nov 16, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using only the 4 corner points might underestimate the real reprojected extents, see e.g. https://github.com/OSGeo/grass/blob/main/lib/proj/do_proj.c#L124 and https://github.com/OSGeo/gdal/blob/master/ogr/ogrct.cpp#L2110 for suggestions to insert some more vertices along the box edges.

Copy link
Contributor Author

@petrasovaa petrasovaa Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was aware of that but given points clouds are typically local scale, I didn't consider the distortion significant. But sure, I can add that.

Copy link
Contributor

@metzm metzm Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True, at local scale the difference between reprojecting only corner points and some more points along the edges of the bbox should be minimal. But using some more points along the edges comes at no cost and could improve accuracy of the estimate.

void get_reprojected_extent(pdal::SpatialReference &spatial_reference,
Copy link
Member

@wenzeslaus wenzeslaus Oct 7, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice to document this either with function doc or in-code comments just to provide a quick reference on the approach taken here (create bbox as points and reproject it using PDAL functions).

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
// Add single corner point and interior points for each edge
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);
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;
}
}
}

#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);
Expand All @@ -97,7 +216,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)
Expand All @@ -123,7 +242,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);
Expand Down
7 changes: 5 additions & 2 deletions raster/r.in.pdal/info.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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_R_IN_PDAL_USE_NOSRS 1
#endif

extern "C" {
Expand All @@ -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);
Expand All @@ -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
Loading
Loading