1

grid_canopy() fails in catalog_apply() when chunk has only 1 vegetation return. (In the case of this file (see code), not dropping the overlap will allow the script to complete, but I don't think this is relevant.) There's a reference to a fix in lidR 2.0 so it must be a different problem than that one.

error message: Interpolation failed (NAs everywhere). Input parameters might be wrong.

I thought I could write code to ignore the error: chm[chm < 0] <- NaN but this didn't work.

Can anyone see how I can ignore the NAs?

# In Nova Scotia elevation explorer, search for:
# 595_4993_201801.laz

library(lidR)
las <- readLAS('d:/working/test/595_4993_201801.laz', 
               select = "xyzcr", 
               filter = "-drop_withheld -drop_overlap -drop_z_below -5")
               # this works: filter = "-drop_withheld -drop_z_below -5")

interpolationMethod <- tin('', extrapolate = knnidw(k = 10, p = 2, rmax = 50))
lasNorm <- normalize_height(las, algorithm = interpolationMethod, na.rm = TRUE, use_class = c(2L, 9L))
lasNormFltr <- filter_poi(lasNorm, ((Classification == 3L | Classification == 4L | Classification ==5L) & (Z >= 0 & Z <= 100)) | (Classification == 2L | Classification == 9L)) #, buffer == 0)  #! LASFilter is depracated. use Filter_POi

dir.create('d:/working/test/norm/')

writeLAS(lasNormFltr, 'd:/working/test/norm/595_4993_201801_NORM.las')

ctgNorm <- readLAScatalog('d:/working/test/norm/',
                          progress = TRUE,
                          select = "xyzcr",
                          filter = "-drop_withheld -drop_overlap -drop_z_below -1") # For las files normalized

ctgNorm@output_options$drivers$Raster$param$overwrite <- TRUE

makeCHM_function = function(chunk) {
  las = readLAS(chunk)
  if (is.empty(las))
    return(NULL)
  
  # Since chunks are smaller than tiles it is possible not to have vegetation.
  # Ensure at least some vegetation
  vegetation <- filter_poi(las, (Classification == 3L | Classification == 4L | Classification == 5L) & Z > 0) #tile fails if veg 0
  if (is.empty(vegetation))
    return(NULL)

  chm <- grid_canopy(las, 1, pitfree(thresholds = c(0,2,5,10,15), max_edge = c(4, 2.5), subcircle = 0.3))
  #chm[chm < 0] <- NaN   # this isn't a solution in this case.
  #NAvalue(chm) <- NaN
  return(chm)
}

opt_chunk_size(ctgNorm) <- 500
opt_chunk_alignment(ctgNorm) <- c(500, 500)
opt_progress(ctgNorm) <- TRUE
opt_output_files(ctgNorm) <- ('d:/working/test/norm/chmNorm_tile_{XLEFT}_{YBOTTOM}')
myCHM <- catalog_apply(ctgNorm, makeCHM_function)
2
  • 3
    You could try something like if (npoints(vegetation) < 3) return(NULL)
    – JRR
    Commented May 13, 2022 at 18:34
  • This line works, thanks. Because it's rare (probably a mis-classification or a lone piece of seaweed) I think it's safe to use this approach.
    – Ray J
    Commented May 13, 2022 at 19:02

1 Answer 1

1

A good work around is the one provided by JRR in the comments.

if (npoints(vegetation) < 3) return(NULL)

This is because it's not the single return that's causing the problem, but the fact that the whole chunk doesn't have any realistic vegetation return in it at all. I.e., when I used a bigger chunk of the same file, which did include more vegetation, it worked even though that isolated return was still there.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.