2

I'm using the R packages lidR and TreeLS to process lidar data and I made a function that takes in a single LAS file and splits out a shapefile, it looks like this:

lasfile = 'tile_1.las'
las = readLAS(lasfile)

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

stemdetection.LAS = function(las, vox, dens, n, rad)
{
  thin = tlsSample(las, voxelize(vox))
  map = treeMap(thin, map.hough(hmin = 0.25, hmax = 16, hstep = 0.5, pixel_size = 0.05, max_radius = rad, min_density = dens))
  tls = stemPoints(las, map)#classify stem points
  seg = stemSegmentation(tls, sgmt.ransac.circle(n))#extract measures
  cs <- crs(tls)
  thin <- NULL
  map <- NULL
  coordinates(seg) = cbind(seg$X, seg$Y)
  output <- writeOGR(seg, td, paste(tools::file_path_sans_ext(basename(lasfile)), "_stems", sep=""), driver="ESRI Shapefile", overwrite_layer=TRUE)
  return(output)
}

I'm trying to change the function to work on a LAScatalog but having some issues with getting the name of each las file (chunk). Here is one example of what I tried:

stemdetection.LAScluster = function(las, vox, dens, n, rad)
{
  las <- readLAS(las, select = 'xyzi')# Read the LAScluster
  if (is.empty(las)) return(NULL)              # Exit early (see documentation)

  thin = tlsSample(las, voxelize(vox))
  map = treeMap(thin, map.hough(hmin = 0.25, hmax = 16, hstep = 0.5, pixel_size = 0.05, max_radius = rad, min_density = dens))
  tls = stemPoints(las, map)#classify stem points
  seg = stemSegmentation(tls, sgmt.ransac.circle(n))#extract measures
  cs <- crs(tls)
  thin <- NULL; map <- NULL ; tls = NULL
  coordinates(seg) = cbind(seg$X, seg$Y)
  output <- writeOGR(seg, td, paste(tools::file_path_sans_ext(basename(ctg@data$File)), "_stems", sep=""), driver="ESRI Shapefile", overwrite_layer=TRUE)
  return(output)
}

stemdetection.LAScatalog = function(las, vox, dens, n, rad)
{
  catalog_apply(las, stemdetection, vox = vox, dens = dens, n = n, rad = rad)
}

Is there a way to append the filename to the LAS object while (or after) performing the readLAS function?

Thanks in advance!

2

There is no (documented) way to get the filename because it does not make sense in the general case. A chunk can be a single file or can be built out of several files and in this case the filename is meaningless. Sometime some chunks may be smaller than a file and thus several chunks share the same filename.

You don't need to call any write*() function. the LAScatalog processing engine already manage that for you. Your code should look like

stemdetection.LAScluster = function(las, vox, dens, n, rad)
{
  bbox <- raster::extent(las)
  las <- readLAS(las) # no select = 'xyzi' it is not considered
  [...]
  coordinates(seg) = cbind(seg$X, seg$Y)

  seg <- raster::crop(seg, bbox) # code to remove the buffer
  return(seg)
}

stemdetection.LAScatalog = function(las, vox, dens, n, rad)
{
  options = list(need_buffer = TRUE, automerge = TRUE)
  catalog_apply(las, stemdetection, vox = vox, dens = dens, n = n, rad = rad, .options = options)
}

Then you can use the engine options

opt_select(ctg) <- "xyzi"
opt_chunk_buffer(ctg) <- 10
opt_chunk_size(ctg) <- 0
opt_output_files(ctg)<- "/path/to/folder/{ORIGINALFILENAME}_stem"

out <- stemdetection(ctg, ...)
| improve this answer | |
  • Thank you once again @JRR. You're not looking for a postdoc by any chance, are you? : - ) – Kartograaf Mar 12 at 19:43

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.