Data Activity: Visualize Elevation Change using LiDAR in R to Better Understand

This tutorial focuses on how to visualize digital elevation models derived from LiDAR data in R. The tutorial is part of the Data Activities that can be used with the Quantifying The Drivers and Impacts of Natural Disturbance Events Teaching Module.

Research Question: How do We Measure Changes in Terrain?

Questions

1. How can LiDAR data be collected?
2. How might we use LiDAR to help study the 2013 Colorado Floods?

Additional LiDAR Background Materials

This data activity below assumes basic understanding of remote sensing and associated landscape models and the use of raster data and plotting raster objects in R. Consider using these other resources if you wish to gain more background in these areas.

Using LiDAR Data

LiDAR data can be used to create many different models of a landscape. This brief lesson on "What is a CHM, DSM and DTM? About Gridded, Raster LiDAR Data" explores three important landscape models that are commonly used.

1. How might we use a CHM, DSM or DTM model to better understand what happened in the 2013 Colorado Flood?
2. Would you use only one of the models or could you use two or more of them together?

In this Data Activity, we will be using Digital Terrain Models (DTMs).

More Details on LiDAR

If you are particularly interested in how LiDAR works consider taking a closer look at how the data are collected and represented by going through this tutorial on "Light Detection and Ranging."

Digital Terrain Models

We can use a digital terrain model (DTM) to view the surface of the earth. By comparing a DTM from before a disturbance event with one from after a disturbance event, we can get measurements of where the landscape changed.

First, we need to load the necessary R packages to work with raster files and set our working directory to the location of our data.

# load libraries
library(raster)   # work with raster files
library(rgdal)    # work with raster files

## Set your working directory to ensure R can find the file we wish to import and where we want to save our files. Be sure to move the downloaded files into your working directory!
wd <- "~/Git/data/" # This will depend on your local environment
setwd(wd)

Then we can read in two DTMs. The first DTM preDTM3.tif is a model from data collected 26-27 June 2013 and the postDTM3.tif is a model from data collected on 8 October 2013.

# Load DTMs into R
DTM_pre <- raster(paste0(wd,"disturb-events-co13/lidar/pre-flood/preDTM3.tif"))
DTM_post <- raster(paste0(wd,"disturb-events-co13/lidar/post-flood/postDTM3.tif"))

# View raster structure
DTM_pre

## class      : RasterLayer 
## dimensions : 2000, 2000, 4e+06  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 473000, 475000, 4434000, 4436000  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs 
## source     : /Users/olearyd/Git/data/disturb-events-co13/lidar/pre-flood/preDTM3.tif 
## names      : preDTM3

DTM_post

## class      : RasterLayer 
## dimensions : 2000, 2000, 4e+06  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 473000, 475000, 4434000, 4436000  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs 
## source     : /Users/olearyd/Git/data/disturb-events-co13/lidar/post-flood/postDTM3.tif 
## names      : postDTM3

Among the information we now about our data from looking at the raster structure, is that the units are in meters for both rasters.

Hillshade layers are models created to add visual depth to maps. It represents what the terrain would look like in shadow with the sun at a specific azimuth. The default azimuth for many hillshades is 315 degrees -- to the NW.

# Creating hillshade for DTM_pre & DTM_post
# In order to generate the hillshde, we need both the slope and the aspect of
# the extent we are working on. 

DTM_pre_slope <- terrain(DTM_pre, opt="slope",units = "radians")
DTM_pre_aspect <- terrain(DTM_pre, opt ="aspect", units = "radians")
DTM_pre_hillshade <- hillShade(DTM_pre_slope, DTM_pre_aspect)

DTM_post_slope <- terrain(DTM_post, opt="slope",units = "radians")
DTM_post_aspect <- terrain(DTM_post, opt ="aspect", units = "radians")
DTM_post_hillshade <- hillShade(DTM_post_slope, DTM_post_aspect)

Now we can plot the raster objects (DTM & hillshade) together by using add=TRUE when plotting the second plot. To be able to see the first (hillshade) plot, through the second (DTM) plot, we also set a value between 0 (transparent) and 1 (not transparent) for the alpha= argument.

# plot Pre-flood w/ hillshade
plot(DTM_pre_hillshade,
        col=grey(1:90/100),  # create a color ramp of grey colors for hillshade
        legend=FALSE,         # no legend, we don't care about the grey of the hillshade
        main="Four Mile Canyon Creek, Boulder County\nPre-Flood",
        axes=FALSE)           # makes for a cleaner plot, if the coordinates aren't necessary

plot(DTM_pre, 
        axes=FALSE,
        alpha=0.3,   # sets how transparent the object will be (0=transparent, 1=not transparent)
        add=TRUE)  # add=TRUE (or T), add plot to the previous plotting frame

# plot Post-flood w/ hillshade
# note, no add=T in this code, so new plotting frame. 
plot(DTM_post_hillshade,
        col=grey(1:90/100),  
        legend=FALSE,
        main="Four Mile Canyon Creek, Boulder County\nPost-Flood",
        axes=FALSE)

plot(DTM_post, 
        axes=FALSE,
        alpha=0.3,
        add=T)

Questions?

1. What does the color scale represent?
2. Can you see changes in these two plots?
3. Zoom in on the main stream bed. Are changes more visible? Can you tell where erosion has occurred? Where soil deposition has occurred?

Digital Elevation Model of Difference (DoD)

A Digital Elevation Model of Difference (DoD) is a model of the change (or difference) between two other digital elevation models - in our case DTMs.

# DoD: erosion to be neg, deposition to be positive, therefore post - pre
DoD <- DTM_post-DTM_pre

plot(DoD,
        main="Digital Elevation Model of Difference (DoD)",
        axes=FALSE)

Here we have our DoD, but it is a bit hard to read. What does the scale bar tell us?

Everything in the yellow shades are close to 0m of elevation change, those areas toward green are up to 10m increase of elevation, and those areas to red and white are a 5m or more decrease in elevation.

We can see a distribution of the values better by viewing a histogram of all the values in the DoD raster object.

# histogram of values in DoD
hist(DoD)

Most of the areas have a very small elevation change. To make the map easier to read, we can do two things.

1. Set breaks for where we want the color to represent: The plot of the DoD above uses a continuous scale to show the gradation between the loss of elevation and the gain in elevation. While this provides a great deal of information, in this case with much of the change around 0 and only a few outlying values near -5m or 10m a categorical scale could help us visualize the data better. In the plotting code we can set this with the breaks= argument in the plot() function. Let's use breaks of -5, -1, -0.5, 0.5, 1, 10 -- which will give use 5 categories.

2. Change the color scheme: We can specify a color for each of elevation categories we just specified with the breaks.

ColorBrewer 2.0 is a great reference for choosing mapping color palettes and provide the hex codes we need for specifying the colors of the map. Once we've chosen appropriate colors, we can create a vector of those colors and then use that vector with the `col=` argument in the `plot()` function to specify these.

Let's now implement these two changes in our code.

# Color palette for 5 categories
difCol5 = c("#d7191c","#fdae61","#ffffbf","#abd9e9","#2c7bb6")

# Alternate palette for 7 categories - try it out!
#difCol7 = c("#d73027","#fc8d59","#fee090","#ffffbf","#e0f3f8","#91bfdb","#4575b4")

# plot hillshade first
plot(DTM_post_hillshade,
        col=grey(1:90/100),  # create a color ramp of grey colors
        legend=FALSE,
        main="Elevation Change Post Flood\nFour Mile Canyon Creek, Boulder County",
        axes=FALSE)

# add the DoD to it with specified breaks & colors
plot(DoD,
        breaks = c(-5,-1,-0.5,0.5,1,10),
        col= difCol5,
        axes=FALSE,
        alpha=0.3,
        add =T)

Question

Do you think this is the best color scheme or set point for the breaks? Create a new map that uses different colors and/or breaks. Does it more clearly show patterns than this plot?

Optional Extension: Crop to Defined Area

If we want to crop the map to a smaller area, say the mouth of the canyon where most of the erosion and deposition appears to have occurred, we can crop by using known geographic locations (in the same CRS as the raster object) or by manually drawing a box.

Method 1: Manually draw cropbox

# plot the rasters you want to crop from 
plot(DTM_post_hillshade,
        col=grey(1:90/100),  # create a color ramp of grey colors
        legend=FALSE,
        main="Four Mile Canyon Creek, Boulder County\nPre-Flood",
        axes=FALSE)

plot(DoD,
        breaks = c(-5,-1,-0.5,0.5,1,10),
        col= difCol5,
        axes=FALSE,
        alpha=0.3,
        add =T)

# crop by designating two opposite corners
cropbox1 <- drawExtent()  

After executing the drawExtent() function, we now physically click on the plot at the two opposite corners of the box you want to crop to. You should see a box drawn on the plot at this point.

When we call this new object, we can view the new extent.

# view the extent of the cropbox1
cropbox1

## class      : Extent 
## xmin       : 473814 
## xmax       : 474982 
## ymin       : 4434537 
## ymax       : 4435390

It is a good idea to write this new extent down, so that you can use the extent again the next time you use the script.

Method 2: Define the cropbox

If you know the desired extent of the object you can also use it to crop the box, by creating an object that is a vector of the x min, x max, y min, and y max of the desired area.

# desired coordinates of the box
cropbox2<-c(473792.6,474999,4434526,4435453)

Once you have the crop box defined, either manually or using the coordinate, you then crop the desired layer to the crop box.

# crop desired layers to this cropbox
DTM_pre_crop <- crop(DTM_pre, cropbox2)
DTM_post_crop <- crop(DTM_post, cropbox2)
DTMpre_hill_crop <- crop(DTM_pre_hillshade,cropbox2)
DTMpost_hill_crop <- crop(DTM_post_hillshade,cropbox2)
DoD_crop <- crop(DoD, cropbox2)

# plot all again using the cropped layers

# PRE
plot(DTMpre_hill_crop,
        col=grey(1:90/100),  # create a color ramp of grey colors
        legend=FALSE,
        main="Four Mile Canyon Creek, Boulder County\nPre-Flood",
        axes=FALSE)
# note \n in the title forces a line break in the title
plot(DTM_pre_crop, 
        axes=FALSE,
        alpha=0.3,
        add=T)

# POST
# plot Post-flood w/ hillshade
plot(DTMpost_hill_crop,
        col=grey(1:90/100),  # create a color ramp of grey colors
        legend=FALSE,
        main="Four Mile Canyon Creek, Boulder County\nPost-Flood",
        axes=FALSE)

plot(DTM_post_crop, 
        axes=FALSE,
        alpha=0.3,
        add=T)

# CHANGE - DoD
plot(DTMpost_hill_crop,
        col=grey(1:90/100),  # create a color ramp of grey colors
        legend=FALSE,
        main="Elevation Change Post Flood\nFour Mile Canyon Creek, Boulder County",
        axes=FALSE)

plot(DoD_crop,
        breaks = c(-5,-1,-0.5,0.5,1,10),
        col= difCol5,
        axes=FALSE,
        alpha=0.3,
        add =T)

Now you have a graphic of your particular area of interest.

Additional Resources

• How could you create a DoD, if you only have LiDAR from a single time point but historical maps? Check out Allen James et al. 2012. Geomorphic change detection using historic maps and DEM differencing: The temporal dimension of geospatial analysis. Geomorphology 137:181-198.