Digital Elevation Model

In the previous lesson, you opened a digital elevation model. The digital elevation model (DEM), also known as a digital terrain model (DTM) represents the elevation of the Earth’s surface. The DEM represents the ground - and thus DOES NOT INCLUDE trees and buildings and other objects.

# load libraries
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.5-12, (SVN revision 1018)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/Rend_Aalamsuz/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/Rend_Aalamsuz/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
# set working directory to ensure R can find the file you wish to import
# setwd("working-dir-path-here")

First, let’s open and plot the digital elevation model.

# open raster data
lidar_dem <- raster(x = "D:/projects/R/DATA/earthanalyticswk3/BLDR_LeeHill/pre-flood/lidar/pre_DTM.tif")

# plot raster data
plot(lidar_dem,
     main = "Lidar Digital Elevation Model (DEM)")

Next, let’s open the digital SURFACE model (DSM). The DSM represents the top of the earth’s surface. Thus, it INCLUDES TREES, BUILDINGS and other objects that sit on the Earth.

# open raster data
lidar_dsm <- raster(x ="D:/projects/R/DATA/earthanalyticswk3/BLDR_LeeHill/pre-flood/lidar/pre_DSM.tif")

# plot raster data
plot(lidar_dsm,
     main = "Lidar Digital Surface Model (DSM)")

Raster Histograms - Distribution of Elevation Values

hist(lidar_dem,
     main = "Distribution of digital Elevation model Values",
     xlab = "Elevation (meters)", ylab = "Frequency",
     col = "lightgreen")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.

hist(lidar_dsm,
     main = "Distribution of digital SURFACE model Values",
     xlab = "Elevation (meters)", ylab = "Frequency",
     col = "lightgreen")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.

Canopy Height Model

The canopy height model (CHM) represents the HEIGHT of the trees. This is not an elevation value, rather it’s the height or distance between the ground and the top of the trees. Some canopy height models also include buildings so you need to look closely at your data to make sure it was properly cleaned before assuming it represents all trees!

Calculate Difference Between Two Rasters

There are different ways to calculate a CHM. One easy way is to subtract the DEM from the DSM.

DSM - DEM = CHM

This math gives you the residual value or difference between the top of the earth surface and the ground which should be the heights of the trees (and buildings if the data haven’t been “cleaned”).

# open raster data
lidar_chm <- lidar_dsm - lidar_dem

# plot raster data
plot(lidar_chm,
     main = "Lidar Canopy Height Model (CHM)")

Plots Using Breaks

#plot raster data
plot(lidar_chm,
     breaks = c(0, 2, 10, 20, 30),
     main = "Lidar Canopy Height Model",
     col = c("white", "brown", "springgreen", "darkgreen"))

Export a Raster We can export a raster file in R using the write.raster() function. Let’s export the canopy height model that we just created to your data folder.

#heck to see if an output directory exists
dir.exists("D:/projects/R/DATA/earthanalyticswk3/BLDR_LeeHill/pre-flood/lidar")
## [1] TRUE
## [1] FALSE

# if the output directory doesn't exist, create it
if (dir.exists("D:/projects/R/DATA/earthanalyticswk3/BLDR_LeeHill/pre-flood/lidar/outputs")) {
  print("the directory exists!")
  } else {
    # if the directory doesn't exist, create it
    # recursive tells R to create the entire directory path ("D:/projects/R/DATA/earthanalyticswk3/BLDR_LeeHill/pre-flood/lidar)
    dir.create("D:/projects/R/DATA/earthanalyticswk3/BLDR_LeeHill/pre-flood/lidar/outputs", recursive = TRUE)
  }
## [1] "the directory exists!"
# export CHM object to new GeotIFF
writeRaster(lidar_chm, "D:/projects/R/DATA/earthanalyticswk3/BLDR_LeeHill/pre-flood/lidar/outputslidar_chm.tiff",
            format = "GTiff",  # output format = GeoTIFF
            overwrite = TRUE) # CAUTION: if this is true, it will overwrite an existing file