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

Define multiple outputs for catalog apply? #404

Closed azh2 closed 3 years ago

azh2 commented 3 years ago

I have the following user-defined catalog_apply function which creates both tree tops and tree crowns, I want to return both, how might I go about doing that?

TreeCrowns <- function(chunk, res, OSGeoPath)
{
las <- readLAS(chunk) #Read the LAScluster
if (is.empty(las)) return(NULL) #Exit early 

Z <- lidR::grid_canopy(las, res, algorithm = lidR::pitfree(thresholds = c(0, 6.56, 16.40, 32.81, 49.21, 65.62, 82.02, 98.42, 114.83, 131.24, 147.64), max_edge = c(9.84252, 6.56168), subcircle = 0.66))

X <- raster::as.matrix(Z)
X <- round(X, digits = 0)
X[is.na(X)] <- 0
Y <- X

for (i in 1:nrow(X))
{ 
   for (j in 1:ncol(X))
   { 
      N = fitlwr(X[i,j])
      M = X[max(1, i-N):min(nrow(X), i+N), max(1, j-N):min(ncol(X), j+N)]
      W = reshape2::melt(M)
      W$d2 = sqrt((W$Var1-mean(W$Var1))^2 + (W$Var2-mean(W$Var2))^2)
      Y[i,j] = max(X[i,j], max(subset(W, d2 <= N, select = value)))
   }
}

raster::values(Z) <- Y
rm(Y)
rm(X)

#Gaussian Smoothing per Qi Chen et al. (2006)
w <- raster::focalWeight(Z, 1, "Gauss")
Z <- raster::focal(Z, w = w, na.rm = TRUE, pad = TRUE)

ttops <- lidR::find_trees(Z, algorithm = alg)

segments <- ForestTools::mcws(treetops = ttops, CHM = Z, minHeight = 1, format = 'raster', OSGeoPath = OSGeoPath)

rm(Z)

Segments <- lidR::merge_spatial(las = las, source = segments, attribute = "WS_ID")
Segments <- lidR::filter_poi(las = Segments, WS_ID != 0 | is.na(WS_ID))
Ras <- lidR::grid_metrics(las = Segments, ~getmode(WS_ID), res = res)
r <- raster::focal(Ras, w = matrix(1, nrow=3, ncol=3), fun = mjrty, pad = TRUE)
crowns <- rasterToPolygons(x = r, na.rm = T, dissolve = T)

gc()

return(crowns)

# What I want
return(list(ttops = ttops, crowns = crowns))

}
Jean-Romain commented 3 years ago

What you want to do is not natively supported but the LAScatalog processing engine has enough flexibility to do it. First read this wiki page. And also ?lidR::lidR-LAScatalog-drivers.

Then, if you are still struggled please ask a question on gis.stackexchange with the tag lidr where I will be able to give you a detailed and comprehensive answer. You must also provide a minimal reproducible example. Your current example is not minimal (80% of the code is useless with regard to your question) and is not reproducible (no example data, no mention of the processing options, ...).

Jean-Romain commented 3 years ago

After thinking about it, you never mentioned that you want to write files. Thus

return(list(ttops = ttops, crowns = crowns))

is expected to work. Anyway please ask on gis.stackexchange being very specific on what you want to achieve.