OSGeo
r.in.pdal: fix reported data extent when CRS doesn't match
This PR tries to fix problems with CRS
- should fix #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
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)
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.
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.
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)
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.
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.
