Crown Segementation on Canopy Height Models from low budget UAV generated digital terrain models
The reliable and reproducible segmentation of tree entities also called tree crowns is the basic need for all individual based determination of structure, competition ecological traits and so on. While this kind of segmentation is even hard in real live there are tons of approaches to perform it on LiDAR data. It is crucial to understand that only high end terrestrial 3D scanning together with terrestrial LiDAR will provide reliable data for doing so. Nevertheless let's try. We have to keep in mind that UAV data is not really 3D because you will find only a small vertical variation in the height of retrieved point in a column. Due to this we are not able to use the well known LiDAR Segmentation approaches on 3D data. We have to rely on a canopy height model (CHM)
Canopy Height Models (CHM) are usually derived from airborne LiDAR data. Since they are rasterized, they can be comprehended as an aggregated simplification of the 3D point cloud model.They are typically representing the tree canopy. A reliable CHM has removed all other above-ground features such as buildings etc.
Even using LiDAR data there are some massive limitations with the most CHMs. Same and even worse with UAV 3D point clouds. Nevertheless CHMs are very useful for all kind of forest and ecological management, for conservation biology and a lot of other topics.
The aim of this short tutorial is to show some possibilities to do so.
To derive Orthoimages or point clouds one can use between different tools the tutorial data is produced using Agisoft Photoscan which is a great tool for deriving point clouds and all kind of surface models. In the end of the processing chain you will have an orthorectified image and a dense point cloud. We will use some UAV data from the Marburg Open Forest project.
First we have to set up the project. Because a bunch of software is needed we create a fix structure and link them against R.
devtools::install_github("gisma/uavRst", ref = "develop")
require(raster)
require(mapview)
require(link2GI)
# proj subfolders
projRootDir<-tempdir()
paths<-link2GI::initProj(projRootDir = projRootDir,
projFolders = c("data/","data/ref/","output/","run/","las/"),
global = TRUE,
path_prefix = "path_")
# get some colors
pal = mapview::mapviewPalette("mapviewTopoColors")
# get the data
# NOTE file size is about 80MB
utils::download.file(url="https://github.com/gisma/gismaData/raw/master/uavRst/data/477369_800_5631924_000_477469_800_5632024_000.las",
destfile=paste0(path_data,"lasdata.las"))
# make the folders and linkages
# NOTE this will take some time - especially running Windows
giLinks<-uavRst::linkAll() We just use or findings from the dsm/dtm tutorial. First we calculate a DSM, second a DTM.
# set the gridsize
actual_grid_size <- 0.5
# n of aspline iterations
lev_max <- 3
# calculate DSM uav point cloud
dsm1 <- uavRst::pc2D_dsm(laspcFile = paste0(path_data,"lasdata.las"),
gisdbasePath = projRootDir,
sampleMethod = "max",
targetGridSize = actual_grid_size,
giLinks = giLinks)## :: create copy of the las file at the working directory...
## :: get extent of the point cloud
## :: link to GRASS
## :: sampling max altitudes using : 0.5 meter grid size
## :: filling no data values if so
raster::plot(dsm1, col = pal(100), main =paste(actual_grid_size,"m DSM "))
# calculate DTM from uav point cloud
dtm1 <- uavRst::pc2D_dtm(laspcFile = paste0(path_data,"lasdata.las"),
gisdbasePath = tempdir(),
tension = 20 ,
sampleGridSize = 25,
targetGridSize = 0.5,
giLinks = giLinks)## :: create copy of the las file at the working directory...
## :: create DTM by interpolation to a raster size of: 0.5
raster::plot(dtm1,col = pal(100), main =paste(actual_grid_size,"m DTM "))
Then we need to define first some arguments. Basically what is the minimum , maximum height of an tree and whats the maximum area a crown can cover.
Second we take the results and simply subtract them from each other.
# setting up the paramters
proj4 = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 "
maxCrownArea = 150
minTreeAlt = 5
# due to the different smoothing approaches we should resample
dsm1 <- raster::resample(dsm1, dtm1 , method = 'bilinear')
# calculate CHM
chmR <- dsm1 - dtm1
# reset negative values to 0
chmR[chmR<0]<-0
raster::plot(chmR,col = pal(100), main =paste(actual_grid_size,"m CHM"))
## Crown Segmentation
Now we can start with the segmentation of the derived CHM. It is twofold again. First we have to identify potential treepositions, second we use this positions to start a seeded region growing tree crown segmentation algorithm.
# identify potential treepositions
tPos <- uavRst::treepos_GWS(chm = chmR,
minTreeAlt = minTreeAlt,
maxCrownArea = maxCrownArea,
join = 1,
thresh = 0.4,
split=TRUE,
cores=4,
giLinks = giLinks )##
## :: start crown identification...
## :: run pre-segmentation...
## :: filter results...
## :: find max height position...
## :: calculate chm statistics
# again we fix spurious location errors
tPos<-raster::resample(tPos, chmR , method = 'bilinear')
# and reset negative tree heights to zero
tPos[tPos<=0]<-0
# writing of the mandantory raster for the seeded segmentation
raster::writeRaster(tPos,paste0(path_run,"treePos.tif"),overwrite = TRUE,NAflag = 0)
raster::writeRaster(chmR,paste0(path_run,"chm.sdat"),overwrite = TRUE,NAflag = 0)# call tree crown segmentation
crownsGWS <- uavRst::chmseg_GWS( treepos = tPos,
chm = chmR,
minTreeAlt = minTreeAlt,
neighbour = 1,
thVarFeature = 1.5,
thVarSpatial = 1.5,
thSimilarity = 0.004,
majorityRadius = 1.0,
minTreeAltParam = "chmQ20",
giLinks = giLinks )## ::: run main segmentation...
## :: calculate chm statistics
## segmentation finsihed...
crownsITC<- uavRst::chmseg_ITC(chm = chmR,
EPSG =3064,
movingWin = 7,
TRESHSeed = 0.45,
TRESHCrown = 0.55,
minTreeAlt = minTreeAlt,
maxCrownArea = maxCrownArea)
raster::plot(chmR, col = pal(50), main =paste(actual_grid_size,"m CHM GWS Algorithm"))
raster::plot(crownsGWS,col=grey(0:100/100,alpha=0.25), add = T)
raster::plot(chmR, col = pal(50), main =paste(actual_grid_size,"m CHM ITC Segment Algorithm"))
raster::plot(crownsITC,col=grey(0:100/100,alpha=0.25), add = T)
The derived CHM is not that bad at all. It shows a pretty fair normalisation of the canopy height even in pretty structured terrain. As you may notice there is an altitude difference of almost 60 m in an area covering 100 by 100 meter. The tree canopy is also pretty structured for a roughly about 120 years old beech wood. Especially low and pretty large poor structured areas are detected poorly. Nevertheless this is just a question of parameterisation. To summarize the result can be compete to the LiDAR derived results.
We want to have a result that keeps the canopy structures as much as possible. This is even more important because of the relative homogeneous UAV point cloud.
In the end it is like most applied science, try it and write it down...