LiDAR data processing
What is LiDAR?
LiDAR (Light detection and ranging) is a data acquisition technology that uses light emitted from a laser pulse (up to 400000 pulses of light per second) to measure the travel time of the pulse from the source to the target and back (returns). Using this approach the instrument collects three dimensional data in large volumes, high density and at unprecedented accuracy (Figure 1). The resolution of the data collected depends on the pulse frequency, that is, resolution increases with increasing pulse frequency.
What types of variables can be extracted from LiDAR data?
LiDAR data is provided in LAS format and classified into ground and non-ground (vegetation) returns. This type of classification is usually performed by the data provider.
Canopy density map -. Forest canopy density is the ratio of the points classified as vegetation to the point classified as ground. This variable, together with canopy height, is commonly used for biomass estimation and the assessment of vegetation coverage.
Digital Terrain Model (DTM) -. The DTM represents the bare ground, it can be calculated from the average value of points that have been classified as ground in each cell.
- Canopy height metrics -. The forest canopy height is calculated by subtracting the DTM from the first return. Other canopy height metrics such as maximum, minimum, mean, standard deviation, variance and percentiles can also be calculated.
Setting the R environment
The package rgrass7 provides the interface between GRASS GIS 7 and the R software that allows the use of all the GRASS GIS commands from the R command line. To load the rgrass7 package into your R session, it is necessary to first install the package using the command install.package(“rgrass7”) and then:
library(rgrass7) #Load the rgrass7 packageIt is recommended to first create a folder where the datasets (LAS files, download a sample dataset from HERE) are located and the outputs and results will be saved. To set the working directory in R to the folder created:
setwd("/home/garzonc/Desktop/DIARS/") #Set the working directoryNow you need to create a directory that will contain the GRASS GIS database for the working session:
dir.create("grassdata") #Create a grassdata directoryTo initiate GRASS into your R session you need to set the gisBase argument to GRASS binaries and libraries directory path, this should be something like C:\Program Files\GRASS GIS 7.00beta4 if your are using Windows (check HERE for other operating systems). This will create a new GRASS location called Sylt_island inside the GRASS GISDBASE.
myGRASS <- "/usr/local/grass-7.0.4svn" #Set the location of the grass binaries and libraries
myGISDbase <- "/home/garzonc/Documents/grassData" #Set the directory path to your GISDbase
myLocation <- "DIARS" #Name the Location
myMapset <- "PERMANENT" #Name the mapset
# Initiate the GRASS GIS session
initGRASS(myGRASS, home = tempdir(), gisDbase = myGISDbase,
location = myLocation, mapset = myMapset,
override = TRUE)## gisdbase /home/garzonc/Documents/grassData
## location DIARS
## mapset PERMANENT
## rows 170
## columns 185
## north 6080248
## south 6079942
## west 463974.6
## east 464307.6
## nsres 1.8
## ewres 1.8
## projection +proj=utm +no_defs +zone=32
## +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000
## +to_meter=1And now set the projection system, for this example we need WGS84 UTM zone 32 (epgs: 32632).
execGRASS("g.proj", flags = c("c"), parameters = list(epsg = 32632)) #Set the projection system## Projection information updated
You can search for the epgs code of your projection system here
To set the geographic region in GRASS including the extent of the bounding box (i.e. the region where all the calculations will be performed) and the resolution:
# Set the spatial region where all
# calculations will be performed
execGRASS("g.region", parameters = list(n = "6075220.40",
s = "6074500.40", e = "454719.89", w = "453999.89",
res = "0.5"))
Importing the LiDAR data into R via GRASS GIS
Once the GRASS in R environment is set you can import the LiDAR data. The following for loop will create 4 vector point files from the binary LAS files using the v.in.lidar GRASS module. Note that this will take a couple of hours to run depending on your computer configuration because this files contain altogether 16568266 points.
# Select the LAS files
for (i in c("GLOBAL-H_454000_6074500.las",
"GLOBAL-H_454000_6075000.las", "GLOBAL-H_454500_6074500.las",
"GLOBAL-H_454500_6075000.las"))
# Create vector files
{
execGRASS("v.in.lidar", flags = c("r",
"o", "overwrite", "quiet"), parameters = list(input = i,
output = substr(i, 8, 23)))
}Now combine the 4 imported vector map layers into one single vector map layer using the v.patch GRASS module (this step is also a time consuming task):
execGRASS("v.patch", flags = c("e", "overwrite",
"quiet"), parameters = list(input = c("H_454000_6074500",
"H_454000_6075000", "H_454500_6074500",
"H_454500_6075000"), output = "las_file"))The points in the LAS files are already classified as vegetation, non-vegetation and, ground. The next step is to extract in a single vector file each of these categories. Start with the points classified as “vegetation” (categories high, medium and low vegetation) and create a new vector map (v.extract GRASS module):
execGRASS("v.extract", flags = c("overwrite",
"quiet"), parameters = list(input = "las_file",
where = "class IN ('Low Vegetation','High Vegetation','Medium Vegetation')",
output = "las_file_vegetation"))Select the points classified as “ground” and create a new vector map:
execGRASS("v.extract", flags = c("overwrite",
"quiet"), parameters = list(input = "las_file",
where = "class IN ('Ground')", output = "las_file_ground"))Select the “non-vegetation” points (categories Unclassified, ground, building and water) and create a new vector map:
execGRASS("v.extract", flags = c("overwrite",
"quiet"), parameters = list(input = "las_file",
where = "class IN ('Unclassified','Ground','Building' ,'Water')",
output = "las_file_non_vegetation"))Export each of the 3 GRASS vector maps containing the vegetation, ground and non-vegetation points to GRASS ASCII map:
execGRASS("v.out.ascii", flags = "overwrite",
parameters = list(input = "las_file_vegetation",
output = "las_file_vegetation", format = "point"))
execGRASS("v.out.ascii", flags = "overwrite",
parameters = list(input = "las_file_non_vegetation",
output = "las_file_non_vegetation",
format = "point"))
execGRASS("v.out.ascii", flags = "overwrite",
parameters = list(input = "las_file_ground",
output = "las_file_ground", format = "point"))You can visualize a subset of the point cloud with the RGL package you need to first install the rgl package using the install.package("rgl") command:
library(rgl) #Load the rgl package
# Import the vector files into R
ground <- read.table("las_file_ground", sep = "|",
dec = ".")
vegetation <- read.table("las_file_vegetation",
sep = "|", dec = ".")
# Select the subset of points to plot
ground <- ground[which(ground[, 1] < 454100 &
ground[, 2] > 6074800 & ground[, 2] <
6074900), ]
vegetation <- vegetation[which(vegetation[,
1] < 454100 & vegetation[, 2] > 6074800 &
vegetation[, 2] < 6074900), ]
# Create the 3d plot
plot3d(ground[, 1], ground[, 2], ground[,
3], col = "coral", aspect = c(1, 1, 0.1),
xlab = "", ylab = "", zlab = "", box = F,
axes = F)
plot3d(vegetation[, 1], vegetation[, 2],
vegetation[, 3], col = "green", aspect = c(1,
1, 0.1), add = T)