Working with National Land Cover Database (NLCD) in R: Raster package and treemaps with ggplot2

This tutorial accompanies the Small Town Big Data blog post Towns, they are a’changin’. Read it here.

Welcome to R. It’s the GIS you never knew you needed. Let’s get set up and do some math with maps.

This tutorial offers an introduction to working with National Land Cover Database rasters in R. We’ll cover:

Loading and plotting a series of rasters
Animating a series of rasters in a GIF file
Basic dataframe manipulation using a for loop to summarize raster data
Treemaps in Ggplot2

Software: If you’re a beginner, I recommend downloading RStudio. Open up a new .R file and start cutting and pasting!

Data: You can peruse and download NLCD data here

The projection is a little funky, and cropping and clipping to a county boundary can take a lot of processing time. To give you a little head start, I’ve projected, clipped and trimmed seven timesteps of Routt County NLCD for you, and included it in this repository. But feel free to download your own area of interest and let me know what you find!

Let’s jump right in.

Loading raster data

# Set your working directory to the /data folder included in the repository using setwd(). 
# Install the following packages for this tutorial: 
# install.packages(c("FedData", "rasterVis", "raster","sp", "rgdal", "reshape2", "treemapify", "ggplot2", "kableExtra", "animation" "scales"))
library(FedData)
library(rasterVis)
library(raster)
library(sp)
library(rgdal)
library(reshape2)
library(treemapify)
library(ggplot2)
library(kableExtra)
library(animation)
library(scales)

# Load a single raster using the raster() function This is the  
r<-raster("nlcd_1.tif")

# What does a raster object look like in R? What are the key attributes of this data format? 
r

## class      : RasterLayer 
## dimensions : 4311, 2513, 10833543  (nrow, ncol, ncell)
## resolution : 30.1, 29.6  (x, y)
## extent     : 289422, 365063.3, 4417034, 4544639  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs 
## source     : C:/Users/ncinglis/Desktop/stbd/nlcd/nlcd_1.tif 
## names      : nlcd_1 
## values     : 11, 95  (min, max)

# We have seven different years of NLCD data. The stack() function takes as many rasters as you give it and stacks them into one object.
# If they have the same resolution, extent and spatial reference, this is an easy way to conduct operations on multiple layers at the same time.

# We can load in the NLCD series straight into a stack.I do this so often that I wrote a little function that takes a regular expression and creates a raster stack from the files that match it.

list2stack <- function(x) {
  z<-list.files(pattern = x)
  y<-raster::stack(z)
  return(y)
}

# If you want to get into regular expressions, there are entire books on it. Just trust me this one works: 
p<-"^nlcd_.*\\.tif$"
s<-list2stack(p)

#Let's plot this stack and take a look at it. 
plot(s)

Making maps: Plotting an NLCD raster with the correct colors

National Land Cover Database has a standardized legend with specific colors that correspond to land cover classes. We want to visualize Routt County land cover in a manner consistent with other NLCD maps. This isn’t exactly trivial in R, but with a few steps, we can plot a land cover map with code that we can use forever and ever!

# Load in the NLCD legend, colors and descriptions from package FedData. 

legend<-pal_nlcd()
legend

##         class code                  description   color
## 1       water   11                   Open Water #476BA0
## 2       water   12           Perennial Ice/Snow #D1DDF9
## 3   developed   21        Developed, Open Space #DDC9C9
## 4   developed   22     Developed, Low Intensity #D89382
## 5   developed   23  Developed, Medium Intensity #ED0000
## 6   developed   24    Developed, High Intensity #AA0000
## 7      barren   31 Barren Land (Rock/Sand/Clay) #B2ADA3
## 8      forest   41             Deciduous Forest #68AA63
## 9      forest   42             Evergreen Forest #1C6330
## 10     forest   43                 Mixed Forest #B5C98E
## 11  shrubland   51                  Dwarf Scrub #A58C30
## 12  shrubland   52                  Scrub/Shrub #CCBA7C
## 13 herbaceous   71         Grassland/Herbaceous #E2E2C1
## 14 herbaceous   72            Sedge/Herbaceuous #C9C977
## 15 herbaceous   73                      Lichens #99C147
## 16 herbaceous   74                         Moss #77AD93
## 17    planted   81                  Pasture/Hay #DBD83D
## 18    planted   82             Cultivated Crops #AA7028
## 19   wetlands   90               Woody Wetlands #BAD8EA
## 20   wetlands   95 Emergent Herbaceous Wetlands #70A3BA

# Love that. I wish I had found FedData 4 years ago. 
# Now, we make a vector of all the values we have in our study area and select those from the legend object. 

vals<-unique(s[[6]])
df<-legend[legend$code %in% vals,]


# Alright, now let's make this plot pretty, working with 2016 data. We're going to use leveplot() from package rasterVis. First, we need to recognize it as a categorical raster using ratify(): 
rat<-ratify(s[[6]])

# I used some code from the creator of rasterVis to make a custom legend:
myKey <- list(rectangles=list(col = df$color),
              text=list(lab=df$description),
              space='left',
              columns=1,
              size=2,
              cex=.6)

# And here we plot it minus the axes, ticks and labels. 
levelplot(rat, att='ID', 
                     col.regions=df$color,
                     par.settings = list(axis.line = list(col = "transparent"), 
                                        strip.background = list(col = 'transparent'), 
                                        strip.border = list(col = 'transparent')), 
                     scales = list(col = "transparent"),
                     colorkey=F,
                     key=myKey)

# Gorgeous. Let's animate it.

Animation

We can animate a series of rasters to show land cover change over time. All we need to do is take that levelplot() function and loop through each layer of our raster stack. We use the package animation to save as a GIF.

# Make a list of years in string format for titles
years.list<-list("2001", "2004", "2006", "2008", "2011", "2013", "2016")


saveGIF(
  {
    for(i in c(1:nlayers(s))){
      rat<-ratify(s[[i]])
      a<-levelplot(rat, att='ID', 
                col.regions=df$color,
                par.settings = list(axis.line = list(col = "transparent"), 
                                    strip.background = list(col = 'transparent'), 
                                    strip.border = list(col = 'transparent')), 
                scales = list(col = "transparent"),
                main=paste0("Routt County land cover ", years.list[[i]]),
                colorkey=F,
)
      print(a)
    }
  }, interval=0.8, movie.name="szoom.gif", ani.width = 600)


# Check your working directory for the gif file. It should look like this:

Working with National Land Cover Database (NLCD) in R: Raster package and treemaps with ggplot2

Small Town Big Data

January 24, 2021

Welcome to R. It’s the GIS you never knew you needed. Let’s get set up and do some math with maps.

Loading raster data

Making maps: Plotting an NLCD raster with the correct colors

Animation