Forest biomass mapping based on mean top-of-canopy height
The study aims to create a forest biomass map in Eberswalde, Germany based on mean top-of-canopy height (TCH) raster model. In the first step, a point cloud of Traunstein forest and inventory data were used to derive TCH-to-biomass relationship in Traunstein forest. We assume that the TCH-to-biomass relationship in Traunstein forest is similar with that of Eberswalde forest. With this assumption, in the second step, using the Traunstein TCH-to-biomass relationship to predict or map biomass in Eberswalde.
Viet Nguyen
Faculty of Forest and Environment Eberswalde University for Sustainable Development
Duc.Nguyen@hnee.de
September, 2022
🛈 Information: This project was hosted at GitHub.
💾 Data: Eberswalde forest point cloud
Traunstein forest point cloud
Traunstein forest inventory
Prepare workspace
📝 note: To ensure reproduce R environment, this project provide renv and R-4.02 folder.
# check working directory
getwd()## [1] "E:/OneDrive_HNEE/01_study/03_Master_FIT/03_semester_3/LiDAR/exercise"# Check path which R loads the packages from
.libPaths()## [1] "C:/Users/viet3/AppData/Local/R/win-library/4.2"
## [2] "C:/Program Files/R/R-4.2.2/library"# change path to renv folder which R should loads packages from
getwd()
.libPaths(new=paste0(getwd(), "/renv/library/R-4.0/x86_64-w64-mingw32"))
.libPaths()
# load packages
library(slidaRtools)library(data.table)## Warning: package 'data.table' was built under R version 4.2.3library(raster)## Warning: package 'raster' was built under R version 4.2.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.2.3library(rgl)## Warning: package 'rgl' was built under R version 4.2.3library(lidR)## Warning: package 'lidR' was built under R version 4.2.3##
## Attaching package: 'lidR'## The following objects are masked from 'package:raster':
##
## projection, projection<-## prepare parallel processing with lidR
# get how many CPU threads lidR using
get_lidr_threads()## [1] 6# set how many threads lidR should use
set_lidr_threads(0.8) # recommended to use maximum 80% of total threads
# check how many CPU threads lidR using again
get_lidr_threads()## [1] 9
1. Derive TCH-to-biomass relationship in Traunstein forest
# load Traunstein point cloud as dataframe
pc.df <- readRDS("data\\Traunstein\\Subplot_PointCloud_Transformed.rds")
head(pc.df)# see point cloud in 3D
display.point.cloud(pc.df, size=2)
# derive canopy height model (CHM) from point cloud
chm.ras <- raster.from.point.cloud(pc.df, res=1, func="max")plot(chm.ras)
chm.ras## class : RasterLayer
## dimensions : 200, 900, 180000 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 0, 900, 0, 200 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : layer
## values : -0.01, 40.49 (min, max)# Convert (melt) the CHM from a raster into to a XYZ-table
chm.df <- ras2xyzdf(chm.ras)head(chm.df)# Assign spatial grid indices (to each CHM pixel) with a plot resolution of 50 m
chm.df$SpatID <- calc.spatial.index(xcor=chm.df$X, ycor=chm.df$Y, res=50)# spatial grid indices columns (SpatID) was added
head(chm.df)tail(chm.df)display.point.cloud(chm.df, col.var="SpatID", col.lim=c(1, 72))
# Aggregation to get mean Z of each 50-m plot (which is TCH)
chm.dt <- data.table(chm.df)
head(chm.dt)agg.chm.dt <- chm.dt[, .(TCH = mean(Z, na.rm=T)), keyby=SpatID]
head(agg.chm.dt)# Load Traunstein inventory data
inv.df <- readRDS("Data\\Traunstein\\Subplot_Inventory.rds")# take a look at the inventory data
head(inv.df)# Assign spatial grid indices (to each tree) with a plot resolution of 50 m
inv.df$SpatID <- calc.spatial.index(xcor=inv.df$X, ycor=inv.df$Y, res=50)# take a look, SpatID column was added
head(inv.df)# Aggregation to get sum of biomass of single trees in each 50-m plot (which is AGB)
inv.dt <- data.table(inv.df)
head(inv.dt)agg.inv.dt <- inv.dt[, .(AGB = sum(AGB, na.rm=T)), keyby=SpatID]
head(agg.inv.dt)# Sort table by column SpatID
setorderv(agg.inv.dt, cols=c("SpatID"))
agg.inv.dt# Convert biomass from tons per quarter hectare to tons per hectare (which is the standard unit for AGB)
agg.inv.dt$AGB <- 4*agg.inv.dt$AGB
agg.inv.dt# Make scatterplot of biomass over height
plot(agg.inv.dt$AGB ~ agg.chm.dt$TCH, xlim=c(0, 30), ylim=c(0, 450))
# Combine the AGB column from the inventory table and the TCH column from the CHM table
metrics.dt <- cbind(agg.inv.dt, TCH=agg.chm.dt$TCH)
head(metrics.dt)# Use nls function to fit power law
nls.AGB.TCH <- nls(AGB ~ a*TCH^b, data=metrics.dt, start=expand.grid(a=0.5, b=2))
a <- coef(nls.AGB.TCH)[1]
b <- coef(nls.AGB.TCH)[2]
curve(a*x^b, add=T, col="red")
2. Preprocessing Eberswalde forest point cloud
# Load the Eberswalde lidar data
ew.las <- lidR::readLAS("data\\Eberswalde\\419500_5853000.laz")## Warning: Invalid data: ScanAngleRank greater than 90 degrees# Get only the lidar data.table from the las object
ew.dt <- ew.las@data
head(ew.dt)# Plot the point cloud
plot(ew.las)
The study area in Eberswalde includes ground and building area which are not our interest for biomass mapping. In the next steps, we classify the point cloud and set height values of points classified as ground or building to 0.
# Classify ground returns
ew.las <- classify_ground(ew.las, algorithm=csf())
# check classification
table(ew.las$Classification)##
## 1 2
## 4427232 1629389# Classify buildings and planar areas
ew.las <- segment_shapes(las=ew.las,
algorithm=shp_plane(k=10),
attribute="planar",
filter= ~Classification != 2L)## Eigenvalues computation: [============--------------------------------------] 25% (9 threads)Eigenvalues computation: [=============-------------------------------------] 26% (9 threads)Eigenvalues computation: [=============-------------------------------------] 27% (9 threads)Eigenvalues computation: [==============------------------------------------] 28% (9 threads)Eigenvalues computation: [==============------------------------------------] 29% (9 threads)Eigenvalues computation: [===============-----------------------------------] 30% (9 threads)Eigenvalues computation: [===============-----------------------------------] 31% (9 threads)Eigenvalues computation: [================----------------------------------] 32% (9 threads)Eigenvalues computation: [================----------------------------------] 33% (9 threads)Eigenvalues computation: [=================---------------------------------] 34% (9 threads)Eigenvalues computation: [=================---------------------------------] 35% (9 threads)Eigenvalues computation: [==================--------------------------------] 36% (9 threads)Eigenvalues computation: [==================--------------------------------] 37% (9 threads)Eigenvalues computation: [===================-------------------------------] 38% (9 threads)Eigenvalues computation: [===================-------------------------------] 39% (9 threads)Eigenvalues computation: [====================------------------------------] 40% (9 threads)Eigenvalues computation: [====================------------------------------] 41% (9 threads)Eigenvalues computation: [=====================-----------------------------] 42% (9 threads)Eigenvalues computation: [=====================-----------------------------] 43% (9 threads)Eigenvalues computation: [======================----------------------------] 44% (9 threads)Eigenvalues computation: [======================----------------------------] 45% (9 threads)Eigenvalues computation: [=======================---------------------------] 46% (9 threads)Eigenvalues computation: [=======================---------------------------] 47% (9 threads)Eigenvalues computation: [========================--------------------------] 48% (9 threads)Eigenvalues computation: [========================--------------------------] 49% (9 threads)Eigenvalues computation: [=========================-------------------------] 50% (9 threads)Eigenvalues computation: [=========================-------------------------] 51% (9 threads)Eigenvalues computation: [==========================------------------------] 52% (9 threads)Eigenvalues computation: [==========================------------------------] 53% (9 threads)Eigenvalues computation: [===========================-----------------------] 54% (9 threads)Eigenvalues computation: [===========================-----------------------] 55% (9 threads)Eigenvalues computation: [============================----------------------] 56% (9 threads)Eigenvalues computation: [============================----------------------] 57% (9 threads)Eigenvalues computation: [=============================---------------------] 58% (9 threads)Eigenvalues computation: [=============================---------------------] 59% (9 threads)Eigenvalues computation: [==============================--------------------] 60% (9 threads)Eigenvalues computation: [==============================--------------------] 61% (9 threads)Eigenvalues computation: [===============================-------------------] 62% (9 threads)Eigenvalues computation: [===============================-------------------] 63% (9 threads)Eigenvalues computation: [================================------------------] 64% (9 threads)Eigenvalues computation: [================================------------------] 65% (9 threads)Eigenvalues computation: [=================================-----------------] 66% (9 threads)Eigenvalues computation: [=================================-----------------] 67% (9 threads)Eigenvalues computation: [==================================----------------] 68% (9 threads)Eigenvalues computation: [==================================----------------] 69% (9 threads)Eigenvalues computation: [===================================---------------] 70% (9 threads)Eigenvalues computation: [===================================---------------] 71% (9 threads)Eigenvalues computation: [====================================--------------] 72% (9 threads)Eigenvalues computation: [====================================--------------] 73% (9 threads)Eigenvalues computation: [=====================================-------------] 74% (9 threads)Eigenvalues computation: [=====================================-------------] 75% (9 threads)Eigenvalues computation: [======================================------------] 76% (9 threads)Eigenvalues computation: [======================================------------] 77% (9 threads)Eigenvalues computation: [=======================================-----------] 78% (9 threads)Eigenvalues computation: [=======================================-----------] 79% (9 threads)Eigenvalues computation: [========================================----------] 80% (9 threads)Eigenvalues computation: [========================================----------] 81% (9 threads)Eigenvalues computation: [=========================================---------] 82% (9 threads)Eigenvalues computation: [=========================================---------] 83% (9 threads)Eigenvalues computation: [==========================================--------] 84% (9 threads)Eigenvalues computation: [==========================================--------] 85% (9 threads)Eigenvalues computation: [===========================================-------] 86% (9 threads)Eigenvalues computation: [===========================================-------] 87% (9 threads)Eigenvalues computation: [============================================------] 88% (9 threads)Eigenvalues computation: [============================================------] 89% (9 threads)Eigenvalues computation: [=============================================-----] 90% (9 threads)Eigenvalues computation: [=============================================-----] 91% (9 threads)Eigenvalues computation: [==============================================----] 92% (9 threads)Eigenvalues computation: [==============================================----] 93% (9 threads)Eigenvalues computation: [===============================================---] 94% (9 threads)Eigenvalues computation: [===============================================---] 95% (9 threads)Eigenvalues computation: [================================================--] 96% (9 threads)Eigenvalues computation: [================================================--] 97% (9 threads)Eigenvalues computation: [=================================================-] 98% (9 threads)Eigenvalues computation: [=================================================-] 99% (9 threads)Eigenvalues computation: [==================================================] 100% (9 threads)# Label all points as buildings, for which 20% or more of their neighbors are planar
metrics <- point_metrics(ew.las, ~list(PlanarNeighborhood=mean(planar)), k=10) ## Metrics computation: [=-------------------------------------------------] 2% (1 threads)Metrics computation: [=-------------------------------------------------] 3% (1 threads)Metrics computation: [==------------------------------------------------] 4% (1 threads)Metrics computation: [==------------------------------------------------] 5% (1 threads)Metrics computation: [===-----------------------------------------------] 6% (1 threads)Metrics computation: [===-----------------------------------------------] 7% (1 threads)Metrics computation: [====----------------------------------------------] 8% (1 threads)Metrics computation: [====----------------------------------------------] 9% (1 threads)Metrics computation: [=====---------------------------------------------] 10% (1 threads)Metrics computation: [=====---------------------------------------------] 11% (1 threads)Metrics computation: [======--------------------------------------------] 12% (1 threads)Metrics computation: [======--------------------------------------------] 13% (1 threads)Metrics computation: [=======-------------------------------------------] 14% (1 threads)Metrics computation: [=======-------------------------------------------] 15% (1 threads)Metrics computation: [========------------------------------------------] 16% (1 threads)Metrics computation: [========------------------------------------------] 17% (1 threads)Metrics computation: [=========-----------------------------------------] 18% (1 threads)Metrics computation: [=========-----------------------------------------] 19% (1 threads)Metrics computation: [==========----------------------------------------] 20% (1 threads)Metrics computation: [==========----------------------------------------] 21% (1 threads)Metrics computation: [===========---------------------------------------] 22% (1 threads)Metrics computation: [===========---------------------------------------] 23% (1 threads)Metrics computation: [============--------------------------------------] 24% (1 threads)Metrics computation: [============--------------------------------------] 25% (1 threads)Metrics computation: [=============-------------------------------------] 26% (1 threads)Metrics computation: [=============-------------------------------------] 27% (1 threads)Metrics computation: [==============------------------------------------] 28% (1 threads)Metrics computation: [==============------------------------------------] 29% (1 threads)Metrics computation: [===============-----------------------------------] 30% (1 threads)Metrics computation: [===============-----------------------------------] 31% (1 threads)Metrics computation: [================----------------------------------] 32% (1 threads)Metrics computation: [================----------------------------------] 33% (1 threads)Metrics computation: [=================---------------------------------] 34% (1 threads)Metrics computation: [=================---------------------------------] 35% (1 threads)Metrics computation: [==================--------------------------------] 36% (1 threads)Metrics computation: [==================--------------------------------] 37% (1 threads)Metrics computation: [===================-------------------------------] 38% (1 threads)Metrics computation: [===================-------------------------------] 39% (1 threads)Metrics computation: [====================------------------------------] 40% (1 threads)Metrics computation: [====================------------------------------] 41% (1 threads)Metrics computation: [=====================-----------------------------] 42% (1 threads)Metrics computation: [=====================-----------------------------] 43% (1 threads)Metrics computation: [======================----------------------------] 44% (1 threads)Metrics computation: [======================----------------------------] 45% (1 threads)Metrics computation: [=======================---------------------------] 46% (1 threads)Metrics computation: [=======================---------------------------] 47% (1 threads)Metrics computation: [========================--------------------------] 48% (1 threads)Metrics computation: [========================--------------------------] 49% (1 threads)Metrics computation: [=========================-------------------------] 50% (1 threads)Metrics computation: [=========================-------------------------] 51% (1 threads)Metrics computation: [==========================------------------------] 52% (1 threads)Metrics computation: [==========================------------------------] 53% (1 threads)Metrics computation: [===========================-----------------------] 54% (1 threads)Metrics computation: [===========================-----------------------] 55% (1 threads)Metrics computation: [============================----------------------] 56% (1 threads)Metrics computation: [============================----------------------] 57% (1 threads)Metrics computation: [=============================---------------------] 58% (1 threads)Metrics computation: [=============================---------------------] 59% (1 threads)Metrics computation: [==============================--------------------] 60% (1 threads)Metrics computation: [==============================--------------------] 61% (1 threads)Metrics computation: [===============================-------------------] 62% (1 threads)Metrics computation: [===============================-------------------] 63% (1 threads)Metrics computation: [================================------------------] 64% (1 threads)Metrics computation: [================================------------------] 65% (1 threads)Metrics computation: [=================================-----------------] 66% (1 threads)Metrics computation: [=================================-----------------] 67% (1 threads)Metrics computation: [==================================----------------] 68% (1 threads)Metrics computation: [==================================----------------] 69% (1 threads)Metrics computation: [===================================---------------] 70% (1 threads)Metrics computation: [===================================---------------] 71% (1 threads)Metrics computation: [====================================--------------] 72% (1 threads)Metrics computation: [====================================--------------] 73% (1 threads)Metrics computation: [=====================================-------------] 74% (1 threads)Metrics computation: [=====================================-------------] 75% (1 threads)Metrics computation: [======================================------------] 76% (1 threads)Metrics computation: [======================================------------] 77% (1 threads)Metrics computation: [=======================================-----------] 78% (1 threads)Metrics computation: [=======================================-----------] 79% (1 threads)Metrics computation: [========================================----------] 80% (1 threads)Metrics computation: [========================================----------] 81% (1 threads)Metrics computation: [=========================================---------] 82% (1 threads)Metrics computation: [=========================================---------] 83% (1 threads)Metrics computation: [==========================================--------] 84% (1 threads)Metrics computation: [==========================================--------] 85% (1 threads)Metrics computation: [===========================================-------] 86% (1 threads)Metrics computation: [===========================================-------] 87% (1 threads)Metrics computation: [============================================------] 88% (1 threads)Metrics computation: [============================================------] 89% (1 threads)Metrics computation: [=============================================-----] 90% (1 threads)Metrics computation: [=============================================-----] 91% (1 threads)Metrics computation: [==============================================----] 92% (1 threads)Metrics computation: [==============================================----] 93% (1 threads)Metrics computation: [===============================================---] 94% (1 threads)Metrics computation: [===============================================---] 95% (1 threads)Metrics computation: [================================================--] 96% (1 threads)Metrics computation: [================================================--] 97% (1 threads)Metrics computation: [=================================================-] 98% (1 threads)Metrics computation: [=================================================-] 99% (1 threads)Metrics computation: [==================================================] 100% (1 threads)ew.las <- add_attribute(ew.las, x=ifelse(metrics$PlanarNeighborhood < 0.2, 0, 1),name="Building")
# check classification
head(ew.las)## class : LAS (v1.2 format 1)
## memory : 4.8 Kb
## extent : 419500, 419500, 5853424, 5853498 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 33N
## area : 0 m²
## points : 6 points
## density : 0 points/m²table(ew.las@data$Building)##
## 0 1
## 5685948 370673table(ew.las@data$planar)##
## FALSE TRUE
## 5854545 202076# plot Eberswalde point cloud with building coloring
plot(ew.las, color = "Building")
# plot Eberswalde point cloud with planar coloring
plot(ew.las, color = "planar")
# Terrain normalize the point cloud
norm.ew.las <- lidR::normalize_height(ew.las, knnidw())## Inverse distance weighting: [=============================================-----] 91% (9 threads)Inverse distance weighting: [==============================================----] 92% (9 threads)Inverse distance weighting: [==============================================----] 93% (9 threads)Inverse distance weighting: [===============================================---] 94% (9 threads)Inverse distance weighting: [===============================================---] 95% (9 threads)Inverse distance weighting: [================================================--] 96% (9 threads)Inverse distance weighting: [================================================--] 97% (9 threads)Inverse distance weighting: [=================================================-] 98% (9 threads)Inverse distance weighting: [=================================================-] 99% (9 threads)Inverse distance weighting: [==================================================] 100% (9 threads)# plot terrain normalized point cloud
plot(norm.ew.las)
## Set all building points to 0 height
# Get the data.table from the las object
norm.ew.dt <- norm.ew.las@data
# Convert to data.frame
norm.ew.df <- data.frame(norm.ew.dt)
# Set all building points to height 0
norm.ew.df[norm.ew.df$Building == 1, "Z"] <- 0
# Set all points planar area to height 0
norm.ew.df[norm.ew.df$planar, 'Z'] =0# Check if the buildings have disappeared
display.point.cloud(norm.ew.df)
# create CHM for Eberswalde forest
ew.chm.ras <- raster.from.point.cloud(norm.ew.dt, res=1, func="max")ew.chm.ras## class : RasterLayer
## dimensions : 500, 500, 250000 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 419500, 420000, 5853000, 5853500 (xmin, xmax, ymin, ymax)
## crs : NA
## source : memory
## names : layer
## values : -0.135, 42.821 (min, max)plot(ew.chm.ras)
3. Using the Traunstein TCH-to-biomass relationship to predict or map biomass in Eberswalde
There are 2 approaches:
3.1. Raster aggregate: the fast solution
This is a very fast solution, because TCH is a raster metric, which can be calculated directly using raster aggregation. For a more general solution using XYZ-table conversion and spatial indexing (second approach). The more general solution can be applied to raster metrics and point cloud metrics (e.g., MCH).
# Aggregate CHM raster for mean top-of-canopy height (TCH)
tch.50m.ras <- raster::aggregate(ew.chm.ras, fact=50, fun=mean)
plot(tch.50m.ras)
# Predict Eberswalde biomass (AGB) from TCH using Traunstein TCH-to-biomass relationship
agb.50m.ras <- a*tch.50m.ras^b
plot(agb.50m.ras)
hist(agb.50m.ras)
# Assign CRS to raster
crs(agb.50m.ras) <- CRS("+init=epsg:32633")## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.3.2. Spatial indexing: the versatile preference
This approach does not use raster aggregation. Instead conversion to XYZ-table and spatial indexing is used. This solution is not restricted to raster metrics, i.e., it also works for making maps from point cloud metrics.
# Convert the CHM to a XYZ-table
ew.chm.df <- ras2xyzdf(ew.chm.ras)head(ew.chm.df)# Calculate spatial ID for the CHM data.frame (i.e. each raster pixel)
ew.chm.df$PlotID <- calc.spatial.index(xcor=ew.chm.df$X, ycor=ew.chm.df$Y, res=50)head(ew.chm.df)# Aggregate the CHM for top-of-canopy height (TCH)
ew.chm.dt <- data.table(ew.chm.df)
agg.ew.tch.dt <- ew.chm.dt[,.(TCH=mean(Z, na.rm=T)), by='PlotID']
agg.ew.tch.dt# Predict biomass from TCH by applying the power law coefficients
# a and b, which we have fitted earlier on the Traunstein data
agg.ew.tch.dt$AGB <- a*agg.ew.tch.dt$TCH^b
head(agg.ew.tch.dt)## Make a map of biomass with each pixel representing 50 m x 50 m
# First make a matrix with the AGB values
agb.mx <- matrix(agg.ew.tch.dt$AGB, nrow=10, ncol=10)
# Convert to raster
agb.50m.ras2 <- raster(agb.mx)
plot(agb.50m.ras2)
# The raster shows the typical 90 degree rotation of matrix-to-raster conversion. Rotate it back by 90 degrees by transposing it (diagonal reflection) and then flipping it upside down (horizontal reflection)
agb.50m.ras2 <- t(agb.50m.ras2)
plot(agb.50m.ras2)
agb.50m.ras2 <- flip(agb.50m.ras2, direction="y")
plot(agb.50m.ras2)
# Set the extent of the raster to correct UTM coordinates
agb.50m.ras2 <- setExtent(agb.50m.ras2, extent(ew.chm.ras))
plot(agb.50m.ras2)
# Assign the coordinate system to the raster
crs(agb.50m.ras2) <- CRS("+init=epsg:32633")
# let compare the biomass map from two approach
par(mfrow = c(1,2)) # create a 1 row x 2 column plotting matrix
plot(agb.50m.ras, main='approach 1')
plot(agb.50m.ras2, main='approach 2')