lasrangecorrection() and sensor_tracking() Errors

[lucas-johnson]

Initially posted on GIS Stack Exchange: https://gis.stackexchange.com/questions/354527/lidr-range-correction-errors

Description

• sensor_tracking() producing warning messages:

> sensor_tracking(lastile)
class       : SpatialPointsDataFrame 
features    : 97 
extent      : 536890.6, 539228.8, 4886784, 4888532  (xmin, xmax, ymin, ymax)
crs         : NA 
variables   : 4
names       :                Z,     gpstime, PointSourceID, npulses 
min values  : 2458.68847461084, 115043501.5,           236,      53 
max values  : 2756.27182237682,   115044948,           238,     297 
Warning message:
7991 pulses with multiple returns were not actually paired. The point cloud is likely to be wrongly populated. These pulses were removed 

• lasrangecorrection() failing due to unrealistic range:

> lasrangecorrection(lastile, sensor, 2000)
An high range R has been computed relatively to the expected average range Rm = 2106
Point number 851196 at (x,y,z,t) = (538442.61, 4888494.26, 548.58, 115043523.41)
Matched with sensor between (538020.25, 4888520.48, 2601.76, 115044213.00) and (538020.03, 4888486.40, 2598.99, 115044213.50)
The range computed was R = 47387.71
Check the correctness of the sensor positions and the correctness of the gpstime either in the point cloud or in the sensor positions.
Error: Unrealistic range: see message above

Data

• The lastile used to produce the errors can be downloaded from the New York State GIS Clearinghouse FTP server: ftp://ftp.gis.ny.gov/elevation/LIDAR/NYSGPO_WarrenWashingtonEssex_2015/u_5370088700_2015.las
• lascheck() yields no immediate redflags:

> lascheck(lastile)

 Checking the data
  - Checking coordinates... ✓
  - Checking coordinates type... ✓
  - Checking attributes type... ✓
  - Checking ReturnNumber validity... ✓
  - Checking NumberOfReturns validity... ✓
  - Checking ReturnNumber vs. NumberOfReturns... ✓
  - Checking RGB validity... ✓
  - Checking absence of NAs... ✓
  - Checking duplicated points...
   ⚠ 1474 points are duplicated and share XYZ coordinates with other points
  - Checking degenerated ground points...
   ⚠ There were 4 degenerated ground points. Some X Y Z coordinates were repeated.
  - Checking attribute population... ✓
  - Checking flag attributes... ✓
 Checking the header
  - Checking header completeness... ✓
  - Checking scale factor validity... ✓
  - Checking point data format ID validity... ✓
  - Checking extra bytes attributes validity... ✓
  - Checking coordinate reference sytem... ✓
 Checking header vs data adequacy
  - Checking attributes vs. point format... ✓
  - Checking header bbox vs. actual content... ✓
  - Checking header number of points vs. actual content... ✓
  - Checking header return number vs. actual content... ✓
 Checking preprocessing already done 
  - Checking ground classification... yes
  - Checking normalization... no
  - Checking negative outliers... ✓
  - Checking flightline classification... yes

[Jean-Romain]

You have the most weirdly invalid file I've ever seen. The good new is that every invalid point is tagged. I don't know how to explain.

First look at

las = readLAS("u_5370088700_2015.las")
las = lasfilterduplicates(las)
fl236 = lasfilter(las, PointSourceID == 236)
plot(fl236, color = "gpstime")

It is a single flightline so the spatial order is the gpstime order. You should have a smooth gradient of color. But you can see some granularity in the gradient.

If we keep only the first and last return of multiple return we should have always two point per pulse

# Get only the first and last returns of multiple returns
fl236fl <- lasfilter(fl236, (ReturnNumber == NumberOfReturns | ReturnNumber == 1) & NumberOfReturns > 1)
fl236fl <- laspulse(fl236fl)

# Should be 2 for each pulse because we kept only first and last return.
# Sometime it is 1 for some point on the very edge
tab = tabulate(fl236fl$pulseID)
incorrect = which(tab > 2)

We can see that after keeping only first and last return (2 points) we still have some pulses with 3 or more points. Let filter these invalid points

invalid = lasfilter(fl236fl, pulseID %in% incorrect)
invalid@data[, .(X,Y,Z, gpstime, UserData,pulseID)]
#>               X       Y      Z   gpstime UserData pulseID
#>     1: 538486.9 4887087 493.93 115043503        0   11427
#>     2: 538481.9 4887087 478.89 115043503        0   11427
#>     3: 538497.9 4887000 480.84 115043503        1   11427
#>     4: 538485.8 4887089 493.67 115043503        0   11449
#>     5: 538481.0 4887088 478.82 115043503        0   11449
#>   ---                                                   
#>199608: 538414.4 4888414 531.37 115043522        1  183387
#>199609: 538409.5 4888499 554.51 115043522        0  183391
#>199610: 538407.4 4888499 549.00 115043522        0  183391
#>199611: 538420.6 4888414 547.65 115043522        1  183391
#>199612: 538414.7 4888414 531.45 115043522        1  183391

We can see that UserData is set to 1 for many of them. Here is the trouble. Points tagged UserData are invalid points with a weird gpstime. Look at the following. The color gradient is smooth

las = readLAS("~/Téléchargements/u_5370088700_2015.las")
las = lasfilterduplicates(las)
las = lasfilter(las, UserData == 0)
fl236 = lasfilter(las, PointSourceID == 236)
plot(fl236, color = "gpstime")

And the sensor tracking runs well

sensor = sensor_tracking(las)
plot(las) %>% add_flightlines3d(sensor)

[Jean-Romain]

sensor_tracking() now handle your case

sensor_tracking(las)
#> Error : After keeping only first and last returns of multiple returns pulses, 7991 pulses still have more than 2 points. This dataset is corrupted and gpstime is likely to be invalid. 

I also updated lascheck() to add a test of gpstime coherence. Sadly this test takes a lot of time and slow down the check. I will try to improve that.

 Checking the data
  - Checking coordinates... ✓
  - Checking coordinates type... ✓
  - Checking attributes type... ✓
  - Checking ReturnNumber validity... ✓
  - Checking NumberOfReturns validity... ✓
  - Checking ReturnNumber vs. NumberOfReturns... ✓
  - Checking RGB validity... ✓
  - Checking absence of NAs... ✓
  - Checking duplicated points...
   ⚠ 1474 points are duplicated and share XYZ coordinates with other points
  - Checking degenerated ground points...
   ⚠ There were 4 degenerated ground points. Some X Y Z coordinates were repeated.
  - Checking attribute population... ✓
  - Checking gpstime incoherances
    ✗ 1639365 pulses (points with the same gpstime) have points with identical ReturnNumber
  - Checking flag attributes... ✓
  - Checking user data attribute...
   ⚠ 3266698 points have a non 0 UserData attribute. This probably has a meaning.
 Checking the header
  - Checking header completeness... ✓
  - Checking scale factor validity... ✓
  - Checking point data format ID validity... ✓
  - Checking extra bytes attributes validity... ✓
  - Checking coordinate reference sytem... ✓
 Checking header vs data adequacy
  - Checking attributes vs. point format... ✓
  - Checking header bbox vs. actual content... ✓
  - Checking header number of points vs. actual content... ✓
  - Checking header return number vs. actual content... ✓
 Checking preprocessing already done 
  - Checking ground classification... yes
  - Checking normalization... no
  - Checking negative outliers... ✓
  - Checking flightline classification... yes

[Jean-Romain]

For lasrangecorrection() I confirm it was a bug with a bad matching at the very edge of a flightline. Hopefully I added safety guards for very weird results.

You must know that the computation is necessarily weaker at the edges of the point cloud because of the lower number of points used to estimate the sensor positionning. Before to compute the range correction check if the positions look good and well aligned especially on the edges. But the good way is to use a large buffer to compute the sensor positions.

[Jean-Romain]

I created a fast version of the gpstime coherence test. Issue closed.

[lucas-johnson]

@Jean-Romain Thank you for investigating this and taking care of it so quickly. This is super helpful. I would like to propose a feature that would help to determine if the intensity values require range-correction. Plotting something like the following (taken from http://www.geogra.uah.es/images/Documentos/emilio/PDF/garcia2010.pdf) would be helpful to diagnose if the intensity values have already been range corrected:

I'm going to try and produce this myself with the help of the output sensor_tracking()

[Jean-Romain]

If you need a function that return the range of each point ask me. This function does not exist because the range is discarded after intensity normalization. But it is few lines of code to change to return the range of each point. Do not try it yourself I already spent a lot of time thinking about special cases (and yet I missed one that was the source of your bug).

[lucas-johnson]

@Jean-Romain Ok - yes it would be very helpful to have a function that computes the range of each point. I'm happy to create a feature request issue if you would like.

[Jean-Romain]

Yes please open a feature request so I can keep tidy logs of what I'm doing. Thanks

[lucas-johnson]

Feature request submitted: #328