Skip to content

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

Merged
cwhite911 merged 11 commits into
OSGeo:mainfrom
petrasovaa:pdal-extent
Aug 24, 2026
Merged

cwhite911 merged 11 commits into
OSGeo:mainfrom
petrasovaa:pdal-extent

Conversation

@petrasovaa

Copy link
Copy Markdown
Contributor

This PR tries to fix problems with CRS

  • should fix [Bug] r.in.pdal does not read *.las file CRS metadata #5791, it should now read the CRS metadata correctly
  • -e flag should work as expected now. It returns the reprojected extent (if needed) rather than the extent in the las header that might not match the CRS of the project. This is done by extracting the bbox and reprojecting the bbox (using pdal pipeline).
    • it should be able to handle multiple files with different CRS, but I didn't test.
  • I changed the default behavior to always reproject when needed, not only with -w flag. Previously, if there would be mismatch, it would just end, so this should behave more consistently with e.g. r.import.
  • removed message about an option that is not there.
  • added test

@github-actions github-actions Bot added raster Related to raster data processing Python Related code is in Python C Related code is in C C++ Related code is in C++ module tests Related to Test Suite labels Oct 7, 2025
marisn
marisn previously requested changes Oct 7, 2025
@petrasovaa

Copy link
Copy Markdown
Contributor Author

The reprojected extent is slightly different on macOS, not sure what to make out of it:

======================================================================
FAIL: test_reprojection_extent (__main__.ExtentTest.test_reprojection_extent)
Test extent is automatically reprojected from epsg 6346
----------------------------------------------------------------------
Traceback (most recent call last):
  File "raster/r.in.pdal/testsuite/test_r_in_pdal_extent.py", line 97, in test_reprojection_extent
    self.assertAlmostEqual(extent["n"], 222012.882582, places=6)
AssertionError: 222013.002911 != 222012.882582 within 6 places (0.12032899999758229 difference)

----------------------------------------------------------------------
Ran 5 tests in 4.077s

FAILED (failures=1)

Comment thread raster/r.in.pdal/info.cpp Resolved
Comment thread raster/r.in.pdal/info.cpp Resolved
@marisn

marisn commented Oct 10, 2025

Copy link
Copy Markdown
Contributor

The reprojected extent is slightly different on macOS, not sure what to make out of it:

======================================================================
FAIL: test_reprojection_extent (__main__.ExtentTest.test_reprojection_extent)
Test extent is automatically reprojected from epsg 6346
----------------------------------------------------------------------
Traceback (most recent call last):
  File "raster/r.in.pdal/testsuite/test_r_in_pdal_extent.py", line 97, in test_reprojection_extent
    self.assertAlmostEqual(extent["n"], 222012.882582, places=6)
AssertionError: 222013.002911 != 222012.882582 within 6 places (0.12032899999758229 difference)

----------------------------------------------------------------------
Ran 5 tests in 4.077s

FAILED (failures=1)

I wouldn't call 0.12 "slightly". It needs to be investigated.

@petrasovaa

petrasovaa commented Oct 10, 2025 •

Copy link
Copy Markdown
Contributor Author

The reprojected extent is slightly different on macOS, not sure what to make out of it:

======================================================================
FAIL: test_reprojection_extent (__main__.ExtentTest.test_reprojection_extent)
Test extent is automatically reprojected from epsg 6346
----------------------------------------------------------------------
Traceback (most recent call last):
  File "raster/r.in.pdal/testsuite/test_r_in_pdal_extent.py", line 97, in test_reprojection_extent
    self.assertAlmostEqual(extent["n"], 222012.882582, places=6)
AssertionError: 222013.002911 != 222012.882582 within 6 places (0.12032899999758229 difference)

----------------------------------------------------------------------
Ran 5 tests in 4.077s

FAILED (failures=1)

I wouldn't call 0.12 "slightly". It needs to be investigated.

I tried couple different things. I tested proj 9.4.1 on linux which is newer than in the ubuntu CI (8.2.1) and older than the mac CI (9.6.2). The results on ubuntu are the same, so it might not be version problem.

In the ncspm test dataset (running in CI on ubuntu) the reprojected bbox values of the points_6346.las are
n=222012.882582 s=222000.743357 e=637511.895399 w=637499.697426 b=100.000000 t=109.000000

when I created a new dataset from EPSG:3358, the results change, they are somewhat closer to the original dataset, approx difference of 0.1 meter:

n=222012.100568 s=221999.961347 e=637512.123476 w=637499.925507 b=100.000000 t=109.000000

The new dataset has PROJ_WKT unlike the one in CI.

Then I tried to compare directly the coordinates of the points once they are reprojected (inside GrassLidarFilter::processOne), they are very close to what they should be (millimeters difference):

637510.0011 221999.9998 100
637504.9966 222012.0019 101
637512.0037 222005.9987 102
637500.0032 222004.0041 103
637502.9953 222009.0028 104
637502.9953 222009.0028 106
637504.0048 222008.0022 106
637511.0015 222000.9993 107
637510.002 222001.9997 108
637504.0039 222006.0023 109

But for some reason the reprojection pipeline transforming the bbox points gives different values.

And none of this seems related to the difference with macOS.

@echoix

echoix commented Oct 10, 2025

Copy link
Copy Markdown
Member

I remember reading something, probably gdal-dev mailing list, and I think you need to look about what proj-data is available in each, not only the version of the software. (And maybe which one is used)

@petrasovaa

Copy link
Copy Markdown
Contributor Author

I remember reading something, probably gdal-dev mailing list, and I think you need to look about what proj-data is available in each, not only the version of the software. (And maybe which one is used)

Right, that could be the difference between platforms. But I am struggling to get the reprojected bbox match exactly the reprojected points.

@metzm

metzm commented Oct 11, 2025

Copy link
Copy Markdown
Contributor

You could try projinfo -o PROJ -s {srs_def} -t {srs_def} on the different systems to check if the exact same projection pipeline is used on the different systems.

Comment thread raster/r.in.pdal/info.cpp Resolved
@petrasovaa
petrasovaa requested a review from metzm November 25, 2025 21:27
Comment thread raster/r.in.pdal/info.cpp Resolved
Comment thread raster/r.in.pdal/info.h Resolved

@cwhite911 cwhite911 left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Looks good. Here are a few additional things Claude picked up.

  1. Add fallback for missing SRS in get_extent.
    When a file's SRS is invalid or empty, the raw header extent gets unioned with reprojected extents from other files, producing a meaningless combined bbox. A G_warning when SRS is missing on the -e path would help; note the import path fatals in this situation, so the two paths disagree quietly.

  2. Documentation not updated. This changes user-facing behavior (auto-reprojection is now the default, -w deprecated) but r.in.pdal.md/.html are untouched. The manual's description of projection handling and -w should be updated.

Comment thread raster/r.in.pdal/info.cpp Resolved
Comment thread raster/r.in.pdal/main.cpp Resolved
@petrasovaa
petrasovaa requested a review from cwhite911 August 24, 2026 18:15
@cwhite911
cwhite911 enabled auto-merge (squash) August 24, 2026 18:28
@github-actions github-actions Bot added HTML Related code is in HTML docs markdown Related to markdown, markdown files labels Aug 24, 2026
@petrasovaa
petrasovaa dismissed stale reviews from marisn and wenzeslaus August 24, 2026 20:51

Adressed

@cwhite911
cwhite911 merged commit 1c68432 into OSGeo:main Aug 24, 2026
26 of 27 checks passed
@github-actions github-actions Bot added this to the 8.6.0 milestone Aug 24, 2026
@cwhite911

Copy link
Copy Markdown
Contributor

@petrasovaa do you want to open an issue to depracte the -w flag in v.in.pdal to keep the behavior consistent?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

C Related code is in C C++ Related code is in C++ docs HTML Related code is in HTML markdown Related to markdown, markdown files module Python Related code is in Python raster Related to raster data processing tests Related to Test Suite

Projects

None yet

Development

Successfully merging this pull request may close these issues.

[Bug] r.in.pdal does not read *.las file CRS metadata

6 participants