The code below presents a tree segmentation workflow from Airborne Laser Scanning (lidar remote sensing) data. The workflow is based on functions from R packages lidaRtRee and lidR, and it includes the following steps:

Steps 1 and 3 are documented in (Monnet et al. 2010; Monnet 2011). The detection performance of this algorithm was evaluated in a benchmark (Eysn et al. 2015).

Licence: CC-BY / source page

Material

Field inventory

The field inventory corresponds to a 50m x 50m plot located in the Chablais mountain (France) (Monnet 2011, 34). All trees with a diameter at breast height above 7.5 cm are inventoried. The data is available in package lidaRtRee.

# load dataset from package (default)
data(treeinventorychablais3, package = "lidaRtRee")

Otherwise you can load your own data provided positions and heights are measured.

# import field inventory
fichier <- "chablais3_listeR.csv"
tree.inventory <- read.csv(file = fichier, sep = ";", header = F)
names(tree.inventory) <- c("x", "y", "d", "h", "n", "s", "e", "t")
# save as rda for later access
# save(tree.inventory,file="tree.inventory.rda")

Attributes are:

  • x easting coordinate in Lambert 93
  • y northing coordinate in Lambert 93
  • d diameter at breast height (cm)
  • h tree height (m)
  • n tree number
  • s species abreviated as GESP (GEnus SPecies)
  • e appearance (0: missing or lying, 1: normal, 2: broken treetop, 3: dead with branches, 4: snag)
  • t tilted (0: no, 1: yes)
##          x       y    d    h n    s e t
## 1 974353.3 6581643 37.6 23.6 1 PIAB 1 0
## 2 974351.0 6581648 15.7 13.9 2 PIAB 1 0
## 3 974348.5 6581650 34.2 23.0 3 PIAB 1 0

Function plotTreeInventory is designed to plot forest inventory data.

# display inventoried trees
lidaRtRee::plotTreeInventory(
  treeinventorychablais3[, c("x", "y")],
  treeinventorychablais3$h,
  species = as.character(treeinventorychablais3$s)
)

We define the region of interest (ROI) to crop ALS data to corresponding extent before further processing. ROI is set on the extent of tree inventory, plus a 10 meter buffer.

# buffer to apply around ROI (meters)
ROI.buffer <- 10
# ROI limits
ROI.range <- data.frame(round(apply(treeinventorychablais3[,c("x","y")],2,range)))

Airborne Laser Scanning data

ALS data is loaded with functions of package lidR. First a catalog of files containing ALS data is built. Then points located inside our ROI are loaded.

# directory for laz files
lazdir <- "/directory_of_las_files/"
# build catalog of files
cata <- lidR::readLAScatalog(lazdir)
# set coordinate system
sp::proj4string(cata) <- sp::CRS(SRS_string = "EPSG:2154")
# set sensor type
lidR::sensor(cata) <- "ALS"
# extract points in ROI plus additional buffer
laschablais3 <- lidR::clip_rectangle(cata,
                                       ROI.range$x[1]-ROI.buffer-5,
                                       ROI.range$y[1]-ROI.buffer-5,
                                       ROI.range$x[2]+ROI.buffer+5,
                                       ROI.range$y[2]+ROI.buffer+5)
# save as rda for easier access:
# save(laschablais3, file="laschablais3.rda", compress = "bzip2")
## class        : LAS (v1.2 format 1)
## memory       : 7 Mb 
## extent       : 974326, 974408, 6581619, 6581702 (xmin, xmax, ymin, ymax)
## coord. ref.  : RGF93 / Lambert-93 
## area         : 6803.003 m²
## points       : 92.1 thousand points
## density      : 13.54 points/m²

Data preparation

Digital Elevation Models

From the ALS point cloud, digital elevation models are computed (Monnet 2011, 43–46). The Digital Terrain Model (DTM) represents the ground surface, it is computed by bilinear interpolation of points classified as ground. The Digital Surface Model (DSM) represents the canopy surface, it is computed by retaining in each raster’s pixel the value of the highest point included in that pixel. The Canopy Height Model (CHM) is the normalized height of the canopy. It is computed by subtracting the DTM to the DSM.

# terrain model computed from points classified as ground
dtm <- lidaRtRee::points2DTM(lidR::filter_ground(laschablais3),
                             res=0.5, ROI.range$x[1]-ROI.buffer,
                             ROI.range$x[2]+ROI.buffer, ROI.range$y[1]-ROI.buffer,
                             ROI.range$y[2]+ROI.buffer)
# surface model
dsm <- lidaRtRee::points2DSM(laschablais3,res=0.5, dtm@extent@xmin,
                             dtm@extent@xmax, dtm@extent@ymin, dtm@extent@ymax)
# canopy height model
chm <- dsm - dtm