lidR – Resolving Issues with Classifying by Shapefile Using merge_spatial in lidR

functionlidrmerger

I have successfully classified buildings via a buildings shapefile input for a single las file but I am having trouble with making a function to use this feature in the catalog format.

Here is the code for a single las file that works well. I exported the mixed conifer las data built into the package and used it for this example. Then in QGIS I drew a heart shaped polygon to use as an example.

# Open R Project, then load packages
#---------------------------------
packages <- c("tidyverse","here","lidR","raster","rgdal","sf","sp", "future")
lapply(packages, library, character.only = TRUE)

# load las
file <- here::here("coconutgrove/mixed_conifer/mixedconifer.las")
ctg <- lidR::catalog(file)
lascrs <- crs(ctg)
# load shp
shp_file <- here("coconutgrove/mixed_conifer/mixed_conifer_heart.shp")
shp <- st_read(shp_file)
shp <- sf::st_transform(shp, crs = st_crs(lascrs))

# =============== single las file - classify buildings by shapefile input 
start.time <- Sys.time() #start time buildings clip

las <- merge_spatial(las, shp, attribute = "in_poly")
las$Classification[las$in_poly] <- LASBUILDING
las$in_poly <- NULL
# end.time <- Sys.time() #end time
time.taken <- end.time - start.time
print(time.taken) #merge spatial

w <- writeLAS(las, here("coconutgrove/las_out/bldgcls_mixed_conifer.las"))

enter image description here

This code successfully executed merge_spatial and classified the points within the polygon as buildings. Below is the code based on (https://cran.r-project.org/web/packages/lidR/vignettes/lidR-catalog-apply-examples.html) to create a user defined function to execute the merge_spatial feature on a las catalog.

# Open R Project, then load packages
#---------------------------------
packages <- c("tidyverse","here","lidR","raster","rgdal","sf","sp", "future")
lapply(packages, library, character.only = TRUE)

# load las
file <- here::here("coconutgrove/mixed_conifer/mixedconifer.las")
ctg <- lidR::catalog(file)
lascrs <- crs(ctg)
# load shp
shp_file <- here("coconutgrove/mixed_conifer/mixed_conifer_heart.shp")
shp <- st_read(shp_file) 
shp <- sf::st_transform(shp, crs = st_crs(lascrs))


# Build the user-defined function that applies merge spatial.
# ================= generic function creation ===============================================

shp_merge_spatial = function(las, ...)
{
  UseMethod("shp_merge_spatial", las)
}


shp_merge_spatial.LAScluster = function(las)
{
  # The function is automatically fed with LAScluster objects
  # Here the input 'las' will a LAScluster
  
  las <- readLAS(las)                          # Read the LAScluster
  if (is.empty(las)) return(NULL)              # Exit early (see documentation)
  
  las <- merge_spatial(las,shp, attribute = "in_poly")      # merge_spatial
  las <- las$Classification[las$in_poly] <- LASBUILDING           # classify
  las$in_poly <- NULL
  return(las)                                  # Return the point cloud
  
}

shp_merge_spatial.LAScatalog = function(las)
{
  catalog_apply(las, shp_merge_spatial.LAScluster)
}

# ======= CATALOG OPTIONS

#opt_filter(ctg)       <- "-drop_z_below 0"
opt_chunk_buffer(ctg) <- 10
opt_chunk_size(ctg)   <- 0
#opt_cores(ctg)        <- 2
opt_output_files(ctg) <- here("coconutgrove/mixed_conifer/ctg_out/{*}_bldg_cls")

start.time <- Sys.time() #start time buildings clip

output <- shp_merge_spatial.LAScluster(ctg)

end.time <- Sys.time() #end time
time.taken <- end.time - start.time
print(time.taken) #merge spatial function

Then the only error message that occurs is:

Warning message:
In las$in_poly <- NULL : Coercing LHS to a list

If I comment out that line of code then there are no error messages yet nothing is written to the output folder I designated in opt_output_files.

Can you advise?

Best Answer

You were almost there but you made some R errors. Below a full reproducible example

shp_merge_spatial = function(las, shp) UseMethod("shp_merge_spatial", las)

shp_merge_spatial.LAScluster = function(las, shp) {
  # The function is automatically fed with LAScluster objects
  # Here the input 'las' will a LAScluster
  
  las <- readLAS(las)                          # Read the LAScluster
  if (is.empty(las)) return(NULL)              # Exit early (see documentation)
  
  las <- merge_spatial(las, shp, attribute = "in_poly")      # merge_spatial
  las$Classification[las$in_poly] <- LASBUILDING           # classify
  return(las)                             # Return the point cloud
}

shp_merge_spatial.LAScatalog = function(las, shp){
  opt_chunk_buffer(las) <- 0
  catalog_sapply(las, shp_merge_spatial, shp = shp)
}

library(lidR)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile, chunk_size = 150)

shp   <- system.file("extdata", "lake_polygons_UTM17.shp", package = "lidR")
lakes <- sf::st_read(shp)

opt_output_files(ctg) <-"{tempdir()}/{ID}_bldg_cls"
output <- shp_merge_spatial(ctg, lakes)

plot(output)
las = readLAS(output)
plot(las, color = "Classification")

This is a bit complex because merge_spatial does not have LAScatalog variant. But in lidR 3.2.0 that will be released probably next week there is classify_poi()

library(lidR)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile, chunk_size = 150)

shp   <- system.file("extdata", "lake_polygons_UTM17.shp", package = "lidR")
lakes <- sf::st_read(shp)

opt_chunk_size(ctg) <- 150
opt_output_files(ctg) <-"{tempdir()}/{ID}_bldg_cls"

output = classify_poi(ctg, LASBUILDING, roi = lakes)

plot(output)
las = readLAS(output)
plot(las, color = "Classification")
Related Question