Fix #286

NEWS.md

@@ -24,6 +24,10 @@
 
 8. Fix `lasmergespatial()` with on disk raster [#285](https://github.com/Jean-Romain/lidR/issues/285).
 
+### ENHANCEMENTS
+
+1. `pitfree()` handle more errors and fails more nicely in some cases [#286](https://github.com/Jean-Romain/lidR/issues/286).
+
 ## lidR v2.1.3 (Release date: 2019-09-10)
 
 #### NEW FEATURES

R/algorithm-dsm.r

@@ -284,29 +284,34 @@ pitfree = function(thresholds = c(0,2,5,10,15), max_edge = c(0,1), subcircle = 0
     data.table::setnames(grid, c("x", "y"), c("X", "Y"))
     cloud <- cloud[cloud[, .I[which.max(Z)], by = cells]$V1]
 
+    if (nrow(cloud) < 3)
+      stop("There are not enought points to triangulate.", call. = FALSE)
+
     # Perform the triangulation and the rasterization (1 loop for classical triangulation, several for Khosravipour et al.)
-    i <- 1
+    i <- 0
+    thresholds <- sort(thresholds)
     for (th in thresholds)
     {
-      verbose(glue::glue("Triangulation pass {i} of {length(thresholds)}..."))
-      i <- i + 1
-
-      if (th == 0)
-        edge <- max_edge[1]
-      else
-        edge <- max_edge[2]
+      verbose(glue::glue("Triangulation pass {i+1} of {length(thresholds)}..."))
 
+      edge  <- if (th == thresholds[1]) max_edge[1] else max_edge[2]
       cloud <- cloud[Z >= th]
 
       if (nrow(cloud) >= 3)
       {
         Ztemp <- interpolate_delaunay(cloud, grid, edge)
-        z     <- pmax(z, Ztemp, na.rm = T)
+
+        if (i == 0 && all(is.na(Ztemp)))
+          stop("Interpolation failed in the first layer (NAs everywhere). Maybe there are too few points.")
+
+        z <- pmax(z, Ztemp, na.rm = T)
       }
+
+      i <- i + 1
     }
 
     if (all(is.na(z)))
-      stop("Interpolation failed. Input parameters might be wrong.")
+      stop("Interpolation failed (NAs everywhere). Input parameters might be wrong.", call. = FALSE)
 
     return(z)
   }

R/grid_canopy.r

@@ -79,6 +79,7 @@ grid_canopy.LAS = function(las, res, algorithm)
   if (is_a_number(res)) assert_all_are_non_negative(res)
   assert_is_algorithm(algorithm)
   assert_is_algorithm_dsm(algorithm)
+  assert_las_is_not_empty(las)
 
   # Some algorithm have an extra option 'subscircle' that need to buffer the layout
   # Must be rewritten because it is a hack !

R/utils_assertive.r

@@ -204,6 +204,12 @@ assert_is_valid_context = function(expected_contexts, name = "", null_allowed =
   return(NULL)
 }
 
+assert_las_is_not_empty = function(x)
+{
+  if (is.empty(x))
+    stop("The point cloud contains 0 point", call. = FALSE)
+}
+
 
 stopifnotlas = function(x)
 {

R/utils_misc.r

@@ -109,7 +109,7 @@ subcircled = function(dt, r, n)
   px = r*cos(alpha)
   py = r*sin(alpha)
 
-  return(dt[, f(X, Y, Z, px, py), by = rownames(dt)][, rownames := NULL][])
+  return(dt[, f(X, Y, Z, px, py), by = 1:nrow(dt)][, nrow := NULL][])
 }
 
 coordinates = function(las)