2

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?

1 Answer 1

3

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")
1
  • looking forward to classify_poi() thanks for the help JRR!! :) Sep 23, 2021 at 19:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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