Tutorial
Raster 06: Plot Raster Time Series Data in R Using RasterVis and Levelplot
Authors: Leah A. Wasser, Megan A. Jones, Zack Brym, Kristina Riemer, Jason Williams, Jeff Hollister, Mike Smorul
Last Updated: Apr 8, 2021
This tutorial covers how to improve plotting output using the rasterVis
package
in R. Specifically it covers using levelplot()
and adding meaningful custom
names to bands within a RasterStack
.
Learning Objectives
After completing this tutorial, you will be able to:
- Be able to assign custom names to bands in a RasterStack for prettier plotting.
- Understand advanced plotting of rasters using the
rasterVis
package andlevelplot
.
Things You’ll Need To Complete This Tutorial
You will need the most current version of R and, preferably, RStudio
loaded
on your computer to complete this tutorial.
Install R Packages
-
raster:
install.packages("raster")
-
rgdal:
install.packages("rgdal")
-
rasterVis:
install.packages("rasterVis")
-
RColorBrewer:
install.packages("RColorBrewer")
-
More on Packages in R – Adapted from Software Carpentry.
Data to Download
NEON Teaching Data Subset: Landsat-derived NDVI raster files
The imagery data used to create this raster teaching data subset were
collected over the
National Ecological Observatory Network's
Harvard Forest
and
San Joaquin Experimental Range
field sites.
The imagery was created by the U.S. Geological Survey (USGS) using a
multispectral scanner
on a Landsat Satellite.
The data files are Geographic Tagged Image-File Format (GeoTIFF).
Set Working Directory: This lesson assumes that you have set your working directory to the location of the downloaded and unzipped data subsets.
An overview of setting the working directory in R can be found here.
R Script & Challenge Code: NEON data lessons often contain challenges that reinforce learned skills. If available, the code for challenge solutions is found in the downloadable R script of the entire lesson, available in the footer of each lesson page.
Get Started
In this tutorial, we are working with the same set of rasters used in the
Raster Time Series Data in R
tutorial. These data are derived from the Landsat satellite and stored in
GeoTIFF
format. Each raster covers the
NEON Harvard Forest field site.
If you have not already created the RasterStack, originally created in Raster Time Series Data in R , please create it now.
# import libraries
library(raster)
library(rgdal)
library(rasterVis)
## Loading required package: terra
## terra version 1.1.4
##
## Attaching package: 'terra'
## The following objects are masked from 'package:tidyr':
##
## expand, fill, pack, separate
## The following object is masked from 'package:zoo':
##
## time<-
## The following object is masked from 'package:grid':
##
## depth
## The following object is masked from 'package:scales':
##
## rescale
## The following object is masked from 'package:ggmap':
##
## inset
## The following object is masked from 'package:rgdal':
##
## project
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, near
## The following object is masked from 'package:knitr':
##
## spin
## Loading required package: lattice
## Loading required package: latticeExtra
##
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
##
## layer
library(RColorBrewer)
# set working directory to ensure R can find the file we wish to import
wd <- "~/Git/data/" # this will depend on your local environment environment
# be sure that the downloaded file is in this directory
setwd(wd)
# Create list of NDVI file paths
all_NDVI_HARV <- list.files(paste0(wd,"NEON-DS-Landsat-NDVI/HARV/2011/NDVI"), full.names = TRUE, pattern = ".tif$")
# Create a time series raster stack
NDVI_HARV_stack <- stack(all_NDVI_HARV)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO",
## prefer_proj = prefer_proj): Discarded datum Unknown based on WGS84
## ellipsoid in CRS definition
# apply scale factor
NDVI_HARV_stack <- NDVI_HARV_stack/10000
Plot Raster Time Series Data
We can use the plot
function to plot our raster time series data.
# view a plot of all of the rasters
# nc specifies number of columns
plot(NDVI_HARV_stack, nc = 4)
Our plot is nice however, it's missing some key elements including, easily
readable titles. It also contains a legend that is repeated for each image. We
can use levelplot
from the rasterVis
package to make our plot prettier!
The syntax for the levelplot()
function is similar to that for the plot()
function. We use main="TITLE"
to add a title to the entire plot series.
# create a `levelplot` plot
levelplot(NDVI_HARV_stack,
main="Landsat NDVI\nNEON Harvard Forest")
Adjust the Color Ramp
Next, let's adjust the color ramp used to render the rasters. First, we
can change the red color ramp to a green one that is more visually suited to our
NDVI (greenness) data using the colorRampPalette()
function in combination with
colorBrewer
.
# use colorbrewer which loads with the rasterVis package to generate
# a color ramp of yellow to green
cols <- colorRampPalette(brewer.pal(9,"YlGn"))
# create a level plot - plot
levelplot(NDVI_HARV_stack,
main="Landsat NDVI -- Improved Colors \nNEON Harvard Forest Field Site",
col.regions=cols)
The yellow to green color ramp visually represents NDVI well given it's a measure of greenness. Someone looking at the plot can quickly understand that pixels that are more green, have a higher NDVI value.
- For all of the
brewer_pal
ramp names see the
Refine Plot & Tile Labels
Next, let's label each raster in our plot with the Julian day that the raster
represents. The current names come from the band (layer names) stored in the
RasterStack
and first the part each name is the Julian day.
To create a more meaningful label we can remove the "x" and replace it with
"day" using the gsub()
function in R. The syntax is as follows:
gsub("StringToReplace","TextToReplaceIt", Robject)
.
First let's remove "_HARV_NDVI_crop" from each label.
# view names for each raster layer
names(NDVI_HARV_stack)
## [1] "X005_HARV_ndvi_crop" "X037_HARV_ndvi_crop" "X085_HARV_ndvi_crop"
## [4] "X133_HARV_ndvi_crop" "X181_HARV_ndvi_crop" "X197_HARV_ndvi_crop"
## [7] "X213_HARV_ndvi_crop" "X229_HARV_ndvi_crop" "X245_HARV_ndvi_crop"
## [10] "X261_HARV_ndvi_crop" "X277_HARV_ndvi_crop" "X293_HARV_ndvi_crop"
## [13] "X309_HARV_ndvi_crop"
# use gsub to modify label names.that we'll use for the plot
rasterNames <- gsub("X","Day ", names(NDVI_HARV_stack))
# view Names
rasterNames
## [1] "Day 005_HARV_ndvi_crop" "Day 037_HARV_ndvi_crop"
## [3] "Day 085_HARV_ndvi_crop" "Day 133_HARV_ndvi_crop"
## [5] "Day 181_HARV_ndvi_crop" "Day 197_HARV_ndvi_crop"
## [7] "Day 213_HARV_ndvi_crop" "Day 229_HARV_ndvi_crop"
## [9] "Day 245_HARV_ndvi_crop" "Day 261_HARV_ndvi_crop"
## [11] "Day 277_HARV_ndvi_crop" "Day 293_HARV_ndvi_crop"
## [13] "Day 309_HARV_ndvi_crop"
# Remove HARV_NDVI_crop from the second part of the string
rasterNames <- gsub("_HARV_ndvi_crop","",rasterNames)
# view names for each raster layer
rasterNames
## [1] "Day 005" "Day 037" "Day 085" "Day 133" "Day 181" "Day 197"
## [7] "Day 213" "Day 229" "Day 245" "Day 261" "Day 277" "Day 293"
## [13] "Day 309"
Once the names for each band have been reassigned, we can render our plot with
the new labels using names.attr=rasterNames
.
# use level plot to create a nice plot with one legend and a 4x4 layout.
levelplot(NDVI_HARV_stack,
layout=c(4, 4), # create a 4x4 layout for the data
col.regions=cols, # add a color ramp
main="Landsat NDVI - Julian Days \nHarvard Forest 2011",
names.attr=rasterNames)
We can adjust the columns of our plot too using layout=c(cols,rows)
. Below
we adjust the layout to be a matrix of 5 columns and 3 rows.
# use level plot to create a nice plot with one legend and a 4x4 layout.
levelplot(NDVI_HARV_stack,
layout=c(5, 3), # create a 5x3 layout for the data
col.regions=cols, # add a color ramp
main="Landsat NDVI - Julian Days \nHarvard Forest 2011",
names.attr=rasterNames)
Finally, scales
allows us to modify the x and y-axis scale. Let's simply
remove the axis ticks from the plot with scales =list(draw=FALSE)
.
# use level plot to create a nice plot with one legend and a 4x4 layout.
levelplot(NDVI_HARV_stack,
layout=c(5, 3), # create a 5x3 layout for the data
col.regions=cols, # add a color ramp
main="Landsat NDVI - Julian Days \nHarvard Forest 2011",
names.attr=rasterNames,
scales=list(draw=FALSE )) # remove axes labels & ticks
- Create a plot and label each tile "Julian Day" with the julian day value following.
- Change the colorramp to a divergent brown to green color ramp to represent the data. Hint: Use the brewerpal page to help choose a color ramp.
Questions: Does having a divergent color ramp represent the data better than a sequential color ramp (like "YlGn")? Can you think of other data sets where a divergent color ramp may be best?