Movatterモバイル変換


[0]ホーム

URL:


Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

Sign up
Appearance settings

Segment individual trees and compute metrics

Jean-Romain edited this pageJun 11, 2020 ·7 revisions

Document updated on February 13th 2020 and up-to-date withlidR v3.0.0

The following example shows how to process your file entirely in the R environment using thelidR package. The file used in this work can be downloadedhere. This work is the "lidR version" of the blog article published onquantitative ecology.org. We will use the same file to do (almost) the same things but using fewer lines of code.

Preprocessing

Read the file

First we need to read the file. For this we only need the spatial coordinates.ScanAngle,Intensity and so on are not useful. WhileClassification of the points is important for computing the digital terrain model, this task was not already done so we must perform that task ourselves.

library(lidR)las= readLAS("Example.las")plot(las)

Classify ground points

classify_ground provides a several algorithm to classify ground points. This function is convenient for small-to-medium plots like the one we are processing but I do not recommend it for large catalog of file files. Here we use thecsf algorithm because it works well without need to tune the parameters.

las= classify_ground(las, csf())plot(las,color="Classification")

Height normalize the dataset

We need to set the ground at 0. We could subtract the DTM to obtain ground points at 0 but here we won't use a DTM but we will rather interpolate each point exactly.

It is important to notice here that neither the ground classification nor the DTM interpolation where performed using a buffer around the region of interest. Thus we can expect a poor height normalization at the egdes of the dataset. But we will consider that fair enough for the purpose of this tutorial.

las= normalize_height(las, tin())plot(las)

Tree segmentation

There are several methods to segment the tree inlidR the following will use a watershed, that is far to be the best, but is good enough for this easy to segment example.

Compute a canopy height model

In the next steps we will use an algorithm that requires an canopy height model. This step can be skipped if you chose an algorithm that performed the segmentation at the point cloud level. So, let's compute a digital surface model with the pit-free algorithm (see alsocanopy height models in lidR).

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))

Optionally we can smooth this CHM using theraster package

# smoothing post-process (e.g. two pass, 3x3 median convolution)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

Segment the trees

The segmentation can be achieved withsegment_trees. Here I chose thewatershed algorithm with a threshold of 4 meters. The point cloud has been updated and each point now has a number that refers to an individual tree (treeID). Points that not trees are assigned the id valueNA.

algo= watershed(chm,th=4)las= segment_trees(las,algo)# remove points that are not assigned to a treetrees= filter_poi(las,!is.na(treeID))plot(trees,color="treeID",colorPalette= pastel.colors(100))

Compute some metrics and hulls

hulls= delineate_crowns(las,func=.stdmetrics)spplot(hulls,"Z")

Deal with a raster

In the previous example, even if the segmentation is done using a canopy height model, the classification has been made on the point cloud. This is because lidR is a point cloud oriented library. But one may want to get the raster to work with rasters. In that case the functionwatershed can be used standalone:

crowns= watershed(chm,th=4)()plot(crowns,col= pastel.colors(100))

Once you are working with rasters the lidR package is not implied anymore. User can rely on theraster package for further analysis. For example:

contour= rasterToPolygons(crowns,dissolve=TRUE)plot(chm,col= height.colors(50))plot(contour,add=T)

Clone this wiki locally


[8]ページ先頭

©2009-2025 Movatter.jp