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

[petrasovaa]

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

[petrasovaa]

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)

[wenzeslaus]

Unrelated, but I got totally confused by PDAL_USE_NOSRS searching for it in PDAL doc or code. A less officially looking name would have avoid my confusion, but no I see it was added in #4750.

It is still needed @marisn? I can't really tell quickly...new library...old format combo or what...

[marisn]

Yes, it is. It prevents failures when run with older PDAL versions. We can remove this option when we explicitly state PDAL minimum version being >=2.4.3.

Feel free to rename it, as it is defined in local header file:

grass/raster/r.in.pdal/info.h

Line 34 in a92a3a0

#define PDAL_USE_NOSRS 1

[wenzeslaus]

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).

[marisn]

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]

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]

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]

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]

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.