r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
607 stars 132 forks source link

rumple_index fails in grid_metrics #143

Closed spono closed 6 years ago

spono commented 6 years ago

ciao JR, rumple_index works fine when used on a whole plot (point-cloud or raster) but fails when used as a function in grid_metrics

r = grid_metrics(las, func=rumple_index(X, Y, Z), res=20)

, giving the following error:

QH6022 qhull input error: 1'th dimension's new bounds [-0.5, 0.5] too wide for
existing bounds [5.1e+06, 5.1e+06]
Error in geometry::delaunayn(X[, 1:2], options = "QbB") : 
  Received error code 1 from qhull.

any idea?

Jean-Romain commented 6 years ago

This is an error from the Qhull library that is wrapped in R by the geometry package. To be honest I have never tried to run rumple_index in grid_canopy. You should be very careful with what you're doing. I may or may not have sense.

This is probably due to a specific set of points that the Qull library doest like. Too few points in a given cells? Let try

r = grid_metrics(las, func=if(.N > 3) rumple_index(X, Y, Z), res=20)

If not we must figure out which pixel failed.

rumple_index_debug = function(X,Y,Z)
{
  tryCatch(
  {
    rumple_index(X,Y,Z)
  },
  error = function(e)
  {
    print(dput(data.frame(X,Y,Z)))
    stop(e)
  })
}

r = grid_metrics(las, func=rumple_index_debug(X, Y, Z), res=20)
spono commented 6 years ago

if(.N > 3) solves the issue.

Is this something that the user should take care of or is it something that should be implemented at the function (i.e. grid_metrics) level?

BTW, thanks!

Jean-Romain commented 6 years ago

Something that should be implemented in rumple_index :wink:

Jean-Romain commented 6 years ago

Fixed. Returns NA if the number of points is < 3. This way it will not fail. However, I guess I should also check for degenerated points...

spono commented 6 years ago

gosh, that might have been just one of the errors. I tried later on a catalog of 4 tiles and then on the single file derived from their merging. This time I get:

Progress: 52%QH6114 qhull precision error: initial simplex is not convex. Distance=0
Error in geometry::delaunayn(X[, 1:2], options = "QbB") : 
  Received error code 2 from qhull.
Jean-Romain commented 6 years ago

I'm not surprised. This one I guess is a pixel with 3 (or more) aligned points or something like that.

rumple_index was not designed to work with grid_metrics so I didn't think about all these issues. But after all why not. We will try to find all the case that should return NA before to fail. Could you run the debug function.

Jean-Romain commented 6 years ago

I wrapped the triangulation in a tryCatch block. It should never fail anymore. It it fails NA is returned, a warning is raised instead of an error. This way it is non blocking for the subsequent pixels.

Also I will add a progress bar in grid_metrics in version 1.6.0 :wink:

spono commented 6 years ago

great! is the trycatchblock already merged into master?

BTW, here the result from the debug function:

Progress: 52%QH6114 qhull precision error: initial simplex is not convex. Distance=0
structure(list(X = c(383739.939999983, 383739.939999983, 383739.909999983, 
383739.999999983), Y = c(5070082.17000001, 5070082.17000001, 
5070083.21000001, 5070083.85000001), Z = c(1754.38003417969, 
1752.19003417969, 1756.07003417969, 1756.24003417969)), .Names = c("X", 
"Y", "Z"), row.names = c(NA, -4L), class = "data.frame")
         X       Y       Z
1 383739.9 5070082 1754.38
2 383739.9 5070082 1752.19
3 383739.9 5070083 1756.07
4 383740.0 5070084 1756.24
Jean-Romain commented 6 years ago

As I guessed, your issue come from degenerated points. Here the issue is that points 1 & 2 have the same X Y coordinates. The qhull library supports that if the good options are passed. This is fixable but so far degenerated points are tested only when interpolating the terrain. Not in rumple index.

The current fix is available in the master branch and will return NA for all the problematic cases.

Jean-Romain commented 6 years ago

I put that here

lassp = lasfiltersurfacepoints(las, 2)
r = grid_metrics(lassp, func=rumple_index(X, Y, Z), res=20)

Is probably a better way to compute an almost fair rumple index