OSGeo / grass

GRASS GIS - free and open-source geospatial processing engine
https://grass.osgeo.org
Other
842 stars 308 forks source link

r.in.pdal handling the number of returns #3827

Closed mazingaro closed 4 months ago

mazingaro commented 4 months ago

Zero-indexed values of the return number

r.in.pdal reports zero-indexed values of the return number. The first return is idexed as return[0], second as return[1] etc. r.in.pdal -p input=point_cloud.laz

I am appending a pdal info --stats report.

Using PDAL library version '2.7.1 (git-version: Release)'
File: point_cloud.laz
File version = 1.4
File signature: LASF
File source ID: 0
Global encoding: 17
Project UUID: 00000000-0000-0000-0000-000000000000
System ID: AL;
Software ID: TerraScan
Creation DOY: 34
Creation Year: 2024
VLR offset (header size): 375
VLR Count: 2
Point format: 8
Point offset: 729
Point count: 23595745
Point count by return[0]: 11930597
Point count by return[1]: 7817264
Point count by return[2]: 3105005
Point count by return[3]: 660896
Point count by return[4]: 80761
Point count by return[5]: 1200
Point count by return[6]: 22
Point count by return[7]: 0
Point count by return[8]: 0
Point count by return[9]: 0
Point count by return[10]: 0
Point count by return[11]: 0
Point count by return[12]: 0
Point count by return[13]: 0
Point count by return[14]: 0
Scales X/Y/Z: 0.00025/0.00025/0.00025
Offsets X/Y/Z: 500000/-0/-0
Max X/Y/Z: 479000/71000/808.982
Min X/Y/Z: 478000/69999.9/498.449
Ext. VLR offset: 0
Ext. VLR count: 0
Compressed: true
Dimensions: X, Y, Z, Intensity, ReturnNumber, NumberOfReturns, ScanDirectionFlag, EdgeOfFlightLine, Classification, Synthetic, KeyPoint, Withheld, Overlap, ScanAngleRank, UserData, PointSourceId, GpsTime, ScanChannel, Red, Green, Blue, Infrared, Deviation

The code in question seems to be

https://github.com/OSGeo/grass/blob/main/raster/r.in.pdal/info.cpp#L118C1-L118C67

Return numbers have a range from 1 - 5 instead of the full range.

When using r.in.pdal to create a raster of return numbers it imports just the first 5. The produces map is correct but it lacks returns. r.in.pdal -e -o input=point_cloud.laz output=point_cloud_return_number type=CELL resolution=1 dimension=returns Majorly we find 5 reurns from LiDAR scanners. Modern RIEGL scanners can have up to 7 returns. PDAL shows 7 return values.

marisn commented 4 months ago

Could you, please, clarify what exactly is the issue? Zero based indexing in LAS info output? Or is there a problem when importing a LAS file? If so, could you share an example that fails to import correctly?

We do not interpret return count / numbers and just use plain values as provided by PDAL library.

mazingaro commented 4 months ago

Hi @marisn. The issue is that when you want to check (-p flag, not import ) your point cloud (header) with r.in.pdal, you get the back the number of returns from 0 onward. As you can see in the firs post, the first entry of the returns is: Point count by return[0]: 11930597 This is the information for the first return, not the zeroth return. If a user wants to check the number of fourth returns, it will probably focus on reading Point count by return[4]: 80761, which in reality shows the number of fifth returns. The PDAL library returns JSON outputs, which means, that this is parsed? I am appending a part of pdal info, which doesnt show any return of zero value `pdal info --all 'point_cloud.laz'


      {
        "average": 1.326737844,
        "count": 22632487,
        "maximum": 7,
        "minimum": 1,
        "name": "ReturnNumber",
        "position": 4,
        "stddev": 0.6448665976,
        "variance": 0.4158529287
      }
rdzur commented 4 months ago

Hello @mazingaro, you mentioned that the "point_cloud_return_number" map contains values 1 through 5, so I'm wondering if those are still the correct returns for return numbers 1 through 5 in your return number map; and I'm curious if it's possible that the 6th and 7th returns (return[5]: 1200 & return[6]: 22) may not be represented in the map, potentially due to averaging at (I'm assuming) 1 meter resolution given their particular spatial representation with the other points?

I also tried to replicate this with a small test point cloud (640k points) that I pulled down from https://usgs.entwine.io

grass -c EPSG:6531 -e ~/grassdata/GRASS_6531/ pdal translate ept://https://s3-us-west-2.amazonaws.com/usgs-lidar-public/USGS_LPC_NM_Roosevelt_Curry_2015_LAS_2017 clovis.laz --readers.ept.bounds="([-11487004.214643, -11486566.316047], [4088336.255706, 4088874.486840])" reprojection --filters.reprojection.in_srs="EPSG:3857+5703" --filters.reprojection.out_srs="EPSG:6531+6360" --writers.las.minor_version=4 -v 8

This point cloud has 6 returns beginning with return[0] when running r.in.pdal -p clovis.laz on GRASS version: 8.3.2 on Platform: macOS-14.2.1-arm64-arm-64bit ; however the sixth return is only one point in this case so I just assumed it would not be represented when importing with r.in.pdal -e -o -n input=clovis.laz output=clovis_return_number type=CELL resolution=3.2808333333 dimension=returns --overwrite

r.info -r map=clovis_return_number produces min=1 max=5 which seemed like a reasonable result portraying what appear to be the correct return numbers albeit for 1 through 5.

Incidentally, in running this sequence on another machine with the same version of GRASS but a different Mac OS (macOS-14.5-arm64-arm-64bit), the above r.in.pdal import command returns a "Segmentation fault: 11". I realize this may be an entirely separate issue, but I was wondering if anyone else may be having trouble with r.in.pdal on Sonoma 14.5 as it seems I am have the same segmentation fault issue on other point clouds with Mac OS 14.5, but no issue with r.in.pdal in 14.2.1.

mazingaro commented 4 months ago

Hi @rdzur, yes: the maps are correct. As mentioned above: the produced map is correct but it lacks returns values above 5. So the first is a report problem and the second a value problem. I probably should opened two issues to avoid confusion. When creating raster maps, lasgrid from LAStools also cuts the values. It is maybe a "historical" version problem which could mean that there were not algorithms updates for such things from point record version changes from LAS 1.4 v6 onward:

  1. The Return Number is the pulse return number for a given output pulse. A given output laser pulse can have many returns, and they must be marked in sequence of return. The first return will have a Return Number of one, the second a Return Number of two, and so on up to fifteen returns. This is for the LAS 1.4 V6. Page 27.
  2. The Return Number is the pulse return number for a given output pulse. A given output laser pulse can have many returns, and they must be marked in sequence of return. The first return will have a Return Number of one, the second a Return Number of two, and so on up to five returns. The Return Number must be between 1 and the Number of Returns, inclusive. This is for the LAS 1.4 V0. Page 17.

PDAL reports 7 return numbers in my data, which is correct.

rdzur commented 4 months ago

Hi @mazingaro, if you reduce your raster resolution, to say 0.33, can you get all seven return values represented in your output return_number raster?

In testing with the clovis.laz file, for example, I modified the r.in.pdal command resolution to 1 ftUS r.in.pdal -e -o -n input=clovis.laz output=clovis_return_number type=CELL resolution=1 dimension=returns --overwrite

After making this modification, the output return number raster represented all 6 returns from the input point cloud with r.info -r map=clovis_return_number providing min=1 and max=6. This change, however, does result in a map with sparse data interspersed with nodata values given the point cloud density versus raster resolution. If you need to approximate a more visually continuous map albeit with interpolated return number counts which may not work for your application, you could also run r.fill.stats on the return raster, e.g. r.fill.stats -k --overwrite input=clovis_return_number output=clovis_return_number_1ftUS_fill distance=4 mode=mode power=2.0 cells=8

marisn commented 4 months ago
  1. The Return Number is the pulse return number for a given output pulse. A given output laser pulse can have many returns, and they must be marked in sequence of return. The first return will have a Return Number of one, the second a Return Number of two, and so on up to fifteen returns. This is for the LAS 1.4 V6. Page 27.

OK, I see now. In the code is a simple iterator and it is 0-based. It will be an easy change to print out n+1.

@rdzur MacOS issue is unrelated to this bug report. Please check that GRASS has been compiled with the PDAL library in use (and all of their dependencies are recompiled to correct versions too).

neteler commented 3 months ago

@mazingaro The second part above "Return numbers have a range from 1 - 5 instead of the full range." should go into a new issue.

marisn commented 3 months ago

@mazingaro The second part above "Return numbers have a range from 1 - 5 instead of the full range." should go into a new issue.

I'm not certain if this is a bug (or do we want to tackle this). Looking at the PDAL code, it assumes for LAS files with version being >1.4 to always follow 15 return definition. As GRASS just reads the data provided by PDAL, doing otherwise would require to implement some kind of logic to what to display if there are no 15 returns in the data. Not worth it.

mazingaro commented 3 months ago

@rdzur for getting a whole, smooth map it would be maybe better to get more return maps (different resolutions) and patch them? But anyway, the information where there is no returns is also important when inspecting.

@rdzur ,@neteler and @marisn I agree that the second part should be another issue. There we can continue this talk and if we find it not worthy also close it soon.