4

I am using lidR package to segment my trees using lastrees directly from the points clouds.

The function works fine, however I would like to export to vector file the canopy boundaries from tree_segmentation that I did directly from my point cloud.

Is it possible to extract canopy boundaries from segmentation lastrees in order to analyze my results I got?

las = lastrees(lasn,li2012(R = 3, speed_up = 5))
plot(las, color = "treeID", colorPalette = col2)
Andre Silva
  • 10,259
  • 12
  • 54
  • 106
Juan Guerra
  • 147
  • 1
  • 8
  • 1
    Do you mean you want vectors defined by the points of each colour blob obtained when you do a plot like plot(las$X, las$Y, col=las$treeID)? Then you want to maybe look at the convex hull or the convex "alpha shape" defined by each group of points with the same treeID... – Spacedman Feb 14 '19 at 16:01
  • 2
    It is unclear what you are asking. What do you mean by export to vector file the canopy boundaries? Please edit your question to clarify. Does this document help you? – JRR Feb 14 '19 at 16:07
  • 1
    Looks like tree_hulls is what I was thinking of. That can be written to a shapefile or geopackage using writeOGR from the rgdal package. – Spacedman Feb 14 '19 at 16:10

1 Answers1

5

This example is directly from the documentation @JRR provided. You can use writeOGR() from the rgdal package to write the convex hull polygons representing tree canopies to shapefile.

library(lidR)
library(rgdal)

las = readLAS("/path/to/your/points.las")
plot(las)

enter image description here

# Classify ground points
las = lasground(las, csf())
plot(las, color = "Classification")

enter image description here

# Normalize points
las = lasnormalize(las, tin())
plot(las)

enter image description here

# Calculate the canopy height model (CHM)
algo = pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
chm  = grid_canopy(las, 0.5, algo)
plot(chm, col = height.colors(50))

enter image description here

# Optionally smooth the CHM with a 3x3 pixel moving window with a median statistic
ker = matrix(1,3,3)
chm = focal(chm, w = ker, fun = median)
chm = focal(chm, w = ker, fun = median)

plot(chm, col = height.colors(50)) # check the image

enter image description here

algo = watershed(chm, th = 4)
las  = lastrees(las, algo)

# remove points that are not assigned to a tree
trees = lasfilter(las, !is.na(treeID))

plot(trees, color = "treeID", colorPalette = pastel.colors(100))

enter image description here

# Calculate tree metrics and convex hulls
metric = tree_metrics(las, .stdtreemetrics)
hulls  = tree_hulls(las)
hulls@data = dplyr::left_join(hulls@data, metric@data)

spplot(hulls, "Z")

enter image description here

# Write to shapefile
writeOGR(obj=hulls, dsn="/Users/aaron/Desktop/", layer="hulls", driver="ESRI Shapefile") # this is in equal area projection
Aaron
  • 51,658
  • 28
  • 154
  • 317