Vuorijärvensuo Individual Tree Detection Demo
Henri Riihimäki
24 May 2018
Introduction
For running the demo, please download the materials from this link. The project folder should contain itd_demo.las and reference.csv. In addition it should have a bucking folder which contains the MARV42Polkyttaja.exe, made by Ilkka Korpela. For matters related to the script / workflow, please contact henri.riihimaki@helsinki.fi. For matters related to bucking and field reference contact ilkka.korpela@helsinki.fi. The following site description was provided by Ilkka Korpela (edited):
The study site is a drained peatland (cottongrass-bog). Main tree species in the area is Scots pine (Pinus sylvestris). First ditches were done in the 1920s. Clear-cut in 1946, planted in about 1950. Current site type Puolukkaturvekangas in the Finnish classification system.
Reference data collection was done by the students of field course FOR110B in 12 teams of two students, during 2016. The study area size is 132 x 164 meters (.las data clipped approximately). Trees were first mapped in LiDAR+aerial photos - then the tree map was verified and filled with omission trees of sufficient tree dbh (stem diameter). The so called photo-trees have thus an estimate of top XYZ, height (thru the DEM), and species. The species classification in image interpretation was pine, spruce, birch and dead trees. Students visited all trees, observed species, crown status, stem defects and dbh. A portion of trees was sampled and measured for d6 (stem diameter 6 m above ground) and height and height of the lowest living branch.
Airborne laser scanning data was collected with Optech Titan multispectral lidar in 2016: see documentation from http://www.helsinki.fi/%7Ekorpela/HYDE_REF/English_documentation.pdf (p. 74, accessed 23.05.2018).
ITD workflow follows the Andrew Plowright’s ForestTools vignette: https://cran.r-project.org/web/packages/ForestTools/vignettes/treetopAnalysis.html
Disclaimer regarding the use of LAStools. User must have a proper LAStools license to use the program; read LAStools LICENSE terms before running script: http://www.cs.unc.edu/~isenburg/lastools/LICENSE.txt This demo is strictly for educational purposes!
Part 1. Create a canopy height model
First: clear workspace and check that you have everything you need in your workspace:
rm(list=ls(all=T))
list.files()## [1] "itd_demo.las" "itd_demo.Rmd" "ITD_Demo.Rproj" "MARV42polkky"
## [5] "reference.csv"For the first time install the needed packages (e.g. install.packages(“raster”) : and load the libraries
library(raster); library(ForestTools); library(sp); library(rgdal); library(tidyverse); library(sf); library(PerformanceAnalytics); library(gridExtra); library(Metrics);library(RColorBrewer)Next, setup LAStools directory and a function that runs LAStools through R
See the info from the input .las-file
f("lasinfo -i itd_demo.las")## [1] "C:/HYapp/LAStools/bin/lasinfo -i itd_demo.las"We see that the itd_demo.las has already been normalized. Next we create a pit free canopy height model (CHM). Here, we follow this workflow: https://rapidlasso.com/2014/11/04/rasterizing-perfect-canopy-height-models-from-lidar/
f("lasthin -i itd_demo.las -subcircle 0.1 -step 0.25 -highest -o tmp.laz")
f("las2dem -i tmp.laz -elevation -step 0.5 -kill 1.0 -o chm_00.asc")
f("las2dem -i tmp.laz -elevation -drop_z_below 2 -step 0.5 -kill 1.0 -o chm_02.asc")
f("las2dem -i tmp.laz -elevation -drop_z_below 5 -step 0.5 -kill 1.0 -o chm_05.asc")
f("las2dem -i tmp.laz -elevation -drop_z_below 10 -step 0.5 -kill 1.0 -o chm_10.asc")
f("las2dem -i tmp.laz -elevation -drop_z_below 15 -step 0.5 -kill 1.0 -o chm_15.asc")
f("las2dem -i tmp.laz -elevation -drop_z_below 20 -step 0.5 -kill 1.0 -o chm_20.asc")
# f("las2dem -i tmp.laz -elevation -drop_z_below 25 -step 0.5 -kill 1.0 -o chm_25.asc") # use more if needed
f("lasgrid -i chm*.asc -merged -step 0.5 -elevation -highest -o chm_pf.asc")Part 2. Individual Tree Detection with ForestTools
Load the result to R. (Set small marginals)
par(mar = c(0.1, 0.1, 0.1, 0.1))
(chm <- raster("chm_pf.asc"))## class : RasterLayer
## dimensions : 416, 399, 165984 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : 2516759, 2516959, 6857749, 6857957 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : F:\FOR254\ITD_Demo\chm_pf.asc
## names : chm_pf# Check
plot(chm)
Define a function which controls local maxima search window size Here, I define that the search window size (in meters) is one tenth of the tree height
diam <- function(tree_height){ tree_height * 0.1 }
h <- seq(2,36,2)
par(mai = rep(1,4))
plot(diam(h), h, type="l", xlab="window radius, m", main = "Window size for local maxima (treetop) detection")
Next we use ForestTools to locate the treetops and segment the crowns The ForestTools::mcws function uses the local maxima (treetops) as markers which controls the crown segmentation:
ttops <- vwf(chm, winFun = diam, minHeight = 5) # drop trees under 5 m
crowns <- mcws(ttops, chm, minHeight = 5, verbose = FALSE, format = "polygons")See summary for the treetops and the crowns and write the result files:
sp_summarise(ttops)## Value
## TreeCount 1129sp_summarise(crowns, variables = c("height", "crownArea"))## Value
## TreeCount 1129.000000
## heightMean 15.696422
## heightMedian 15.960000
## heightSD 2.737337
## heightMin 5.040000
## heightMax 22.309999
## crownAreaMean 11.275465
## crownAreaMedian 10.500000
## crownAreaSD 5.613104
## crownAreaMin 0.250000
## crownAreaMax 42.750000writeOGR(ttops, getwd(), "treetops", driver = "ESRI Shapefile", overwrite_layer = T)
writeOGR(crowns, getwd(), "crowns", driver = "ESRI Shapefile", overwrite_layer = T)Part 3. Result evaluation
Reference data consist of 1226 trees measured on-site. The main species in the area is pine.
Read the reference data and remove fallen trees, unfound and stump observations
(d <- read_csv2("reference.csv"))## # A tibble: 1,226 x 13
## ID plaji plk elava x y z_foto h_foto d13 d6 h
## <chr> <int> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 11 1 2.52e6 6.86e6 170. 16.2 212 17.3 NA
## 2 2 1 11 1 2.52e6 6.86e6 171. 17.8 187 15.2 NA
## 3 3 1 11 1 2.52e6 6.86e6 173. 19.6 257 21.1 NA
## 4 4 4 11 1 2.52e6 6.86e6 168. 15.0 125 10 14.5
## 5 5 1 11 1 2.52e6 6.86e6 168. 14.1 126 10.0 NA
## 6 6 2 11 1 2.52e6 6.86e6 164. 10.0 120 9.54 NA
## 7 7 4 11 1 2.52e6 6.86e6 167. 13.9 144 11.6 NA
## 8 8 2 11 1 2.52e6 6.86e6 166. 12.2 129 10.3 NA
## 9 9 1 12a1 1 2.52e6 6.86e6 171. 17.2 214 17.5 NA
## 10 10 1 11 1 2.52e6 6.86e6 174. 20.4 252 20.7 NA
## # ... with 1,216 more rows, and 2 more variables: hc <dbl>, h_est <dbl>d$plaji <- as.factor(d$plaji) # convert factor variables to factors
d$plk <- as.factor(d$plk)
d <- subset(d, !(plk %in% c("23", "31", "41")))Plot scatterplots & see correlations
chart.Correlation(d[,7:13])
Next we convert the field reference to a simple feature (spatial) format by using the sf package.
Note: The coordinate system is Finland Zone 2 (EPSG: 2392). The current Finnish coordinate system is EUREF FIN TM35 (EPSG: 3067), the user could transform the data with sf::st_transform function (not used here)
sd <- st_as_sf(d, coords = c("x","y")) %>% st_set_crs(2392) # %>% st_transform(3067))Convert the treetops and crowns to simple feature format, and add an ID column for the crowns
ttops <- st_as_sf(ttops) %>% st_set_crs(2392) # add transform if wanted
crowns <- st_as_sf(crowns) %>% st_set_crs(2392)
crowns$crownID <- row.names(crowns) # id columnCompare the data sets
summary(crowns$height) # ITD## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.04 14.54 15.96 15.70 17.42 22.31summary(d$h_est) # Field## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.28 15.05 16.36 16.27 17.70 22.64Finally, join the two data sets ( see https://en.wikipedia.org/wiki/DE-9IM for further information)
st_join(crowns, sd, join = st_intersects) -> cr_to_ref # crowns intersecting with reference
st_join(sd, crowns, join = st_covered_by) -> ref_to_cr # reference trees covered by a crownCalculate omission and comission statistics Compare the data sets
n_comission <- sum(is.na(cr_to_ref$d13)) # extra crowns that don't have any attributes
n_omission <- sum(is.na(ref_to_cr$height)) # reference trees that are not covered by a crownFind out missed observations because the tree is under another canopy or the canopies are merged i.e. a segmented crown is containing more than one reference observation:
- nrow(crowns) = number of segmented crowns
- nrow(cr_to_ref2) = number of reference observations within the segmented crowns
- difference is the number of extra trees
- add the value to the trees that didn’t weren’t covered by ITD-segments
subset(cr_to_ref, !(is.na(d13))) -> cr_to_ref2 # drop comission trees from this
n_omission <- n_omission + nrow(cr_to_ref2)-nrow(crowns)Outputs: 1. unique crowns with reference data 2. all crowns
itd_res <- subset(ref_to_cr, !(is.na(ref_to_cr$height)))See results:
n_omission/nrow(ref_to_cr)*100## [1] 7.610475n_comission/nrow(ref_to_cr)*100## [1] 10.72013nrow(itd_res)/nrow(ref_to_cr)*100 # crowns which have at least one reference tree## [1] 99.509Part 4 - from height to volume
First estimate the diameter at breast height with a model. Model input is the height of the treetop. The model used is from Kalliovirta & Tokola 2005, http://www.metla.fi/silvafennica/full/sf39/sf392227.pdf
Estimate d13 with model nr. 2.1.1. The study reports a 0.64 correlation for tree-level data. Compare the results to reference data (itd_res)
mod211 <- function(height) {
d13 <- (-0.785 + 1.177*sqrt(height))^2
return(d13)
}
mod211(itd_res$height*10) -> itd_res$d13_est_211
mod211(ttops$height*10) -> ttops$d13_est_211
plot( x = itd_res$d13, y = itd_res$d13_est_211, xlab="Actual", ylab="Predicted", main = "Model 2.1.1")
cor(itd_res$d13, itd_res$d13_est_211, use = "complete.obs")## [1] 0.5966722We can see that the correlation is close to the correlation of the original study.
Next, prepare ITD tree data (ttops) for bucking. Note that the system separator must be “.” The program needs PlotID, TreeID, tree species, d13 and height in an ascii file This version of the script doesn’t have tree species identification included, but the data consist mainly of pine, thus we assume that all the trees are pines (species code = 1). R will open the program as a pop-up window. Close the program after executing it.
#ITD
out_data_itd <- ttops
st_geometry(out_data_itd) <- NULL
out_data_itd$plotID <- 1
out_data_itd$species <- 1
out_data_itd$d13_est_211 <- out_data_itd$d13_est_211/10 # convert to cm
out_data_itd <- dplyr::select(out_data_itd, c(plotID,treeID, species,d13_est_211, height))
head(out_data_itd)## plotID treeID species d13_est_211 height
## 1 1 1 1 23.99394 19.12
## 2 1 2 1 22.71583 18.15
## 3 1 3 1 25.78867 20.48
## 4 1 4 1 24.33681 19.38
## 5 1 5 1 24.12580 19.22
## 6 1 6 1 24.71938 19.67# Reference
out_data_ref <- d
out_data_ref$plotID <- 1
out_data_ref$plaji <- as.numeric(out_data_ref$plaji)
out_data_ref$plaji[is.na(out_data_ref$plaji)] <- 1 # pine by default
out_data_ref$ID <- seq(1,nrow(out_data_ref),1) # orig_id has characters, program requires numeric ID
out_data_ref$d13 <- out_data_ref$d13/10 # convert to cm
out_data_ref <- dplyr::select(out_data_ref, c(plotID, ID, plaji, d13, h_est))
head(out_data_ref) # names are different but variables are the same## # A tibble: 6 x 5
## plotID ID plaji d13 h_est
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 21.2 16.2
## 2 1 2 1 18.7 17.8
## 3 1 3 1 25.7 19.6
## 4 1 4 4 12.5 14.5
## 5 1 5 1 12.6 14.1
## 6 1 6 2 12 10.0Create a directory, note that it must be in C:/data, Then write the output data and execute the bucking program for it, and read the results back to R. Repeat the bucking for both, ITD and reference files: [File -> open -> *executes program -> close program]
if(dir.exists("C:/data/") == F) { dir.create("C:/data/")}
# ITD
write.table(out_data_itd, "C:/data/in_trees.txt", dec=".", sep = ",", row.names = F, col.names = F)
# Run bucking
system("MARV42polkky/MARV42Polkyttaja.exe",wait=T)
results_itd <- read_csv("C:/data/out_trees.txt")
rename(results_itd, treeID = TreeID) %>%
left_join(ttops, ., by ="treeID") -> results_itd
# Reference
write.table(out_data_ref, "C:/data/in_trees.txt", dec=".", sep = ",", row.names = F, col.names = F)
# Run bucking
system("MARV42polkky/MARV42Polkyttaja.exe",wait=T)
results_ref <- read_csv("C:/data/out_trees.txt") Join the itd-data to crown segments
st_join(crowns,results_itd, join = st_intersects) -> itd_res2Finally, plot some of the results
Part 5. Results:
ggplot() +
geom_sf(data = itd_res2, aes(fill = SawVol)) +
scale_fill_distiller(palette = "Greens", direction = 1, breaks = seq(0,700,100))+
ggtitle("Sawlog volume")+
labs(fill="litres")+
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank()) 
ggplot() +
geom_sf(data = itd_res2, aes(fill = PulpVol)) +
scale_fill_distiller(palette = "Greens", direction = 1, limits = c(0,250), breaks = seq(0,250,50))+
ggtitle("Pulpwood volume")+
labs(fill="litres")+
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank()) 
ggplot() +
geom_sf(data = itd_res2, aes(fill = ToTVol)) +
scale_fill_distiller(palette = "Greens", direction = 1, limits = c(0,700), breaks = seq(0,700,100))+
ggtitle("Total volume")+
labs(fill="litres")+
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank())
Lastly, compare bucking results
res_mat <- as.data.frame(matrix(NA,3,4))
colnames(res_mat) <- c("Sawlog volume", "Pulp volume", "Residue volume", "Total volume")
rownames(res_mat) <- c("Reference", "ITD", "Difference")
res_mat[1,1]<- sum(results_ref$SawVol)/1000 #m3
res_mat[1,2]<- sum(results_ref$PulpVol)/1000 #m3
res_mat[1,3]<- sum(results_ref$TopVol)/1000 #m3
res_mat[2,1]<- sum(results_itd$SawVol)/1000 #m3
res_mat[2,2]<- sum(results_itd$PulpVol)/1000 #m3
res_mat[2,3]<- sum(results_itd$TopVol)/1000 #m3
res_mat[1,4] <- sum(res_mat[1,1:3])
res_mat[2,4] <- sum(res_mat[2,1:3])
res_mat[3,1] <- res_mat[1,1]-res_mat[2,1]
res_mat[3,2] <- res_mat[1,2]-res_mat[2,2]
res_mat[3,3] <- res_mat[1,3]-res_mat[2,3]
res_mat[3,4] <- res_mat[1,4]-res_mat[2,4]
area_size <- 2.165 # ha (field plot size 132 m x 164 m, dem cropped approximately to same size, might be small differences)
res_mat # absolute m3## Sawlog volume Pulp volume Residue volume Total volume
## Reference 131.4683 159.8456 5.8000 297.1139
## ITD 140.5440 136.6579 5.0555 282.2574
## Difference -9.0757 23.1877 0.7445 14.8565res_mat/area_size # m3 / ha## Sawlog volume Pulp volume Residue volume Total volume
## Reference 60.724388 73.83169 2.6789838 137.235058
## ITD 64.916397 63.12143 2.3351039 130.372933
## Difference -4.192009 10.71025 0.3438799 6.862125We can see that ITD gives a slight overestimate of sawlog volume, but underestimates the pulpwood and total volume. This can be explained by small omission trees, whereas larger trees are well detected.