Finding the tree tops
Now to build the Canopy Height Model CHM from lidar .tif files using R & QGIS. Starting in R at the last line of previous post;
setwd(file.path("C:/Users/myname/Documents/RWorking"))
R needs libraries loading to manipulate data. R loads missing libraries automatically, but it is best to deliberately load them to make certain). For this process the following are needed;
library(ForestTools)
# ForestTools is the core library for this operation, and its loading initiates loading of the other libraries. .
library(terra)
library(sf)
library(sp)
library(raster)
# R may add library information, but it is usually OK to proceed from here.
My data (from linz.data) comprises two files; a Digital Elevation Model or DEM, and a Digital Surface Model - DSM. DEMs are bare earth/ground height points, while DSMs are points of all items detected above ground. These are height-encoded .tif files, which look like greyscale photos but each pixel contains a height value, instead of a colour (the grayscale effect is a height to grey palette). Here is a DSM (bright dots lower left are an overhead power line), red line is a buffer of site boundary to include all trees that cross the boundary.
GIS files are often comprise many sub-images (tiles) so it is useful to cut out (clip) to the smallest area. This simplifies the interface, and streamlines processing.
To make above image, import DSM files into QGIS, then merge them: Menu > Raster > Miscellaneous > Merge, which brings up the Merge dialog. For most occasions just select files via dialog and go OK, and files merge. The new merged file appears at top of layers panel named Merged. Save with a sensible name (click layer name and press F2 to enter line for renaming), and save by clicking file name in the layers panel, choose Export, which brings up the dialog Save Raster Layer as.... Dialog Settings: Output mode set to Raw data Format should be on Geotiff Save file to your working directory and click OK.
Now to make the CHM. This is made by subtracting the DEM from the DSM using the Raster Calculator* where the lowest point is the ground for each pixel and the upper height is the top of the vegetation canopy.
*QGIS has two dialog calculators and users often try to use the wrong calculator. The correct one for this post is via Menu > Raster > Raster Calculator. Expressions are built by selecting layer names and entering functions either from dialog buttons or keyboard entry.
Sometimes there is a glitch and the minus sign has no space either side, or more than one space. Check before you run the calculation as it's a source of errors.
)
These operations are easier if both images are clipped to the same size. The red line is the land title plus a twenty metre buffer {BOLD} shown by the green line which is used as the clip boundary.
Select polygon layer to make a buffer from Menu > Vector > Geoprocessing Tools > Buffer
and in dialog Check object you want to buffer is selected: Select layer in QGIS layers panel, and click Enter a distance. I used 20m to capture the full width of any tree canopy crossing boundary - check aerial visually (and ground check) to estimate this. End cap and join styles determine whether buffer is rounded or angular. Keep other defaults, and click Okay result should be a larger polygon Save this to your working directory by clicking on layer name in the layer Panel, and choose Export > Save features as.
Now clip image and discard all data outside buffer.
Menu: Raster > Extraction > clip raster by mask layer. Then in dialog under input layer > Select the file you want to clip, and then select both of these:
Match the extent of the dipped raster to the extent of the mask layer
Keep resolution of input raster
Accept all other defaults, click OK and the result will pop up.
In dialog: Clipping extent, Right drop down arrow > Clip raster by layer, and select your mask layer, in this case the title polygon. Ensure 'Resolution of input raster' is selected and press OK.
The next post will be on using the CHM to get a good idea of total number of trees, delineate outline of each tree canopy, and then identify centres of tree tops as in this image
