Chapter 7 Raster Data - Discrete Variables

The previous chapter covered raster datasets with continuous variable, where the data can take any value. This chapter covers raster datasets with discrete variables, where the data are classified into a limited number of values. We still use numbers to work with discrete data. For example, different types of land cover can be identified using integers. Water might be coded as 11, built-up areas as 23, and deciduous forests as 41. In this situation, the numbers themselves are not meaningful, so summary statistics such as sum() and mean() would not provided any useful information. Therefore, we will need to use some different approaches when visualizing and summarizing discrete raster data.

This chapter will use classified land cover data from the National Land Cover Database (NLCD). The NLCD produces a standardized land cover data product for the United States. Data are available from 2001-2016 at 2-3 year intervals. These datasets have been generated from Landsat imagery using consistent methods that make the resulting products suitable for change analysis. More information on the NLCD can be found at the Multi-Resolution Land Characteristics Consortium (MRLC) http://www.mrlc.gov/. Data can be downloaded as a single, large dataset for the conterminous US from the MRLC website, or in smaller chunks via MRLC viewer at https://www.mrlc.gov/viewer/. For this exercise, we will work with a small portion of the NLCD encompassing Walton County, GA. This rural county is located on the outskirts of the Atlanta metropolitan area and as a result is under considerable development pressure.

As usual, we start by loading the required packages.

library(raster)         
library(rgdal)          
library(RColorBrewer) 
library(rasterVis)    
library(ggplot2)      
library(colorspace)
library(dplyr)
library(tidyr)
library(readr)
library(sf)

After the packages have been loaded, they are visible in the in the list of attached packages and objects, which can be viewed with the search() function.

search()
##  [1] ".GlobalEnv"           "package:colorspace"   "package:rasterVis"   
##  [4] "package:latticeExtra" "package:lattice"      "package:raster"      
##  [7] "package:gridExtra"    "package:RColorBrewer" "package:scales"      
## [10] "package:rgdal"        "package:sp"           "package:sf"          
## [13] "package:lubridate"    "package:tidyr"        "package:readxl"      
## [16] "package:readr"        "package:ggplot2"      "package:dplyr"       
## [19] "package:stats"        "package:graphics"     "package:grDevices"   
## [22] "package:utils"        "package:datasets"     "package:methods"     
## [25] "Autoloads"            "package:base"

When a function is called, R goes through all available packages in memory to find the right one. If there are functions with the same name in more than one package, then R will run the function from the first package found in the search list. The other functions are “masked,” meaning they are not called by default. For example, the raster library has a select() function and so does dplyr. If dplyr comes before raster in the search list, then the dplyr select() will always be run, and the raster select() function is always masked. If you want to choose a function from a particular library, you will have to call it explicitly using the double-colon::operator; e.g.dplyr::select()orraster::select()`. Note that the order of packages in the search list is the opposite of the order that they are loaded - the most recently loaded packages mask previously loaded packages.

These function conflicts are a common source of errors in R programming. If you think you are using a the raster select() function to process a raster object and the dplyr select() function is invoked instead, it will return an error. One way to minimize these errors is to load your most important packages last instead of first. Also, if you are using a function with a generic name like select() that is found in multiple packages, it is good practice to call is explicitly with the :: operator. To see if a particular function is present in multiple packages, you can use the help() function with the package name as an argument. If that function is present in two or more loaded packages, RStudio will list them in the Help window. Try this out with help(select).

7.1 Importing and mapping classified satellite imagery

We start by reading the Walton County data and creating a raster object.

nlcd16 <- raster("NLCD_2016_Land_Cover_WaltonGA.tiff")
class(nlcd16)
## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"
nlcd16
## class      : RasterLayer 
## dimensions : 1501, 1708, 2563708  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 1098705, 1149945, 1238565, 1283595  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source     : E:/work/projects/gdswr-book/NLCD_2016_Land_Cover_WaltonGA.tiff 
## names      : NLCD_2016_Land_Cover_WaltonGA 
## values     : 11, 95  (min, max)

Mapping these raster data with ggplot() will require the rasterdf() custom function to convert them to a data frame. To avoid having to paste the code for this function into every script, it can be stored a separate file and run using the source() function. The file containing the function code needs to be in the working directory, or else the file path needs to be explicitly specified.

source("rasterdf.R")
rasterdf
## function (x, aggregate = 1) 
## {
##     resampleFactor <- aggregate
##     inputRaster <- x
##     inCols <- ncol(inputRaster)
##     inRows <- nrow(inputRaster)
##     resampledRaster <- raster(ncol = (inCols/resampleFactor), 
##         nrow = (inRows/resampleFactor))
##     extent(resampledRaster) <- extent(inputRaster)
##     y <- resample(inputRaster, resampledRaster, method = "ngb")
##     coords <- xyFromCell(y, seq_len(ncell(y)))
##     dat <- stack(as.data.frame(getValues(y)))
##     names(dat) <- c("value", "variable")
##     dat <- cbind(coords, dat)
##     dat
## }
nlcd16_df <- rasterdf(nlcd16)

ggplot(data = nlcd16_df) +
  geom_raster(aes(x = x, y = y, fill = value))

This map and the accompanying legend are not particularly useful. The values stored in the land cover raster are numeric codes. For example, a value of 11 represents water wheres a value of 21 represents low-density residential development. The default chart created by ggplot() maps these values on a continuous scale rather than as categories, and the resulting symbology and legend have no real meaning. A detailed list of all the NLCD codes can be found at the following link: https://www.mrlc.gov/data/legends/national-land-cover-database-2011-nlcd2011-legend.

To generate a better-looking map, will start by creating two vectors - one containing the numeric codes for all 16 NLCD classes that occur in the conterminous United States, and another containing a human-readable names in the same order as the codes.

LCcodes <- c(11,12,21,22,23,24,31,41,
             42,43,52,71,81,82,90,95)

LCnames <-c(
  "Water",
  "IceSnow",
  "DevelopedOpen",
  "DevelopedLow",
  "DevelopedMed",
  "DevelopedHigh",
  "Barren",
  "DeciduousForest",
  "EvergreenForest",
  "MixedForest",
  "ShrubScrub",
  "GrassHerbaceous",
  "PastureHay",
  "CultCrops",
  "WoodyWetlands",
  "EmergentHerbWet")

Next, we’ll use a trick to extract the default NLCD symbology from the color table that is embedded in the TIF dataset. Then the NLCD codes are assigned as the names attribute for this vector of colors.

LCcolors <- attr(nlcd16, "legend")@colortable[LCcodes + 1]
names(LCcolors) <- as.character(LCcodes)
LCcolors
##        11        12        21        22        23        24        31        41 
## "#466B9F" "#D1DEF8" "#DEC5C5" "#D99282" "#EB0000" "#AB0000" "#B3AC9F" "#68AB5F" 
##        42        43        52        71        81        82        90        95 
## "#1C5F2C" "#B5C58F" "#CCB879" "#DFDFC2" "#DCD939" "#AB6C28" "#B8D9EB" "#6C9FB8"

Now the ggplot() function can be used to plot the NLCD land cover raster. In the aes() function, the raster value is converted to a character so that ggplot() will recognize is as a categorical variable. Thescale_fill_manual() function is then used to specify the colors that match the NLCD codes. We are removing the second item of LCnames because there are no “Perennial Ice/Snow” areas in Georgia, so we do not want this class to show up in the legend. (You can check this out by running unique(nlcd11) and noting that “12” is not present.) The theme() function is used to remove the axis titles and change the appearance of the background.

ggplot(data = nlcd16_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = LCcolors,
                    labels = LCnames[-2],
                    na.translate = FALSE) +
  coord_sf(expand = F) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

In the last tutorial, we reduced the size of a raster dataset by defining an extent and then using the crop() function. Here, we zoom in and plot a portion of Walton County by setting the x and y limits of the plot using the coord_fixed() function.

ggplot(data = nlcd16_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = LCcolors,
                    labels = LCnames[-2],
                    na.translate = FALSE) +
  coord_sf(expand = F, xlim = c(1114000, 1125000), 
           ylim = c(1260000, 1270000)) + 
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

7.2 Reclassifying raster data

Often, land cover patterns can be seen more clearly if the land cover classes are aggregated into a smaller number of broader land cover types. After listing the classes in the Walton County dataset, a vector is created with the codes for the new land cover classes. These two vectors are combined to create a lookup table that associates the old and new codes. The subs() function is then used to assign a new land cover class to every pixel in Walton County. New vectors are laso created with the names and display colors for the new classes.

oldclas <- unique(nlcd16)
oldclas
##  [1] 11 21 22 23 24 31 41 42 43 52 71 81 82 90 95
newclas <- c(1, 2, 2, 2, 2, 3, 4, 4, 4, 5, 5, 5, 6, 7, 7)
lookup <- data.frame(oldclas, newclas)
nlcd16_rc <- subs(nlcd16, lookup)
newnames <- c("Water",
               "Developed",
               "Barren",
               "Forest",
               "GrassShrub",
               "Cropland",
               "Wetland")
newcols <- c("mediumblue", 
             "red2", 
             "gray60", 
             "darkgreen", 
             "yellow2", 
             "orange4", 
             "paleturquoise2")

The same ggplot() code as before can be used to generate a categorical map with the new reclassified raster dataset.

nlcd16_rc_df <- rasterdf(nlcd16_rc)
ggplot(data = nlcd16_rc_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = newcols,
                    labels = newnames,
                    na.translate = FALSE) +
  coord_sf(expand = F) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

It can be difficult to discern the patterns in a map with too many bright, saturated colors. There are couple of things that can be done about this problem. One is to select different color for each class.A list of all the “named” colors in R is available at the following link: http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf Alternately, a function such as desaturate() from the colorspace package can be used to modify your the colors. The amount argument controls how much saturation is reduced for all the colors in the palette.

newcols2 <- desaturate(newcols, amount = 0.3)

ggplot(data = nlcd16_rc_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = newcols2,
                    labels = newnames,
                    na.translate = FALSE) +
  coord_sf(expand = F) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

Another useful visualization technique is the “small multiples” approach, in which each land cover class is displayed in a separate map. To generate this map, a separate binary raster is generated for each of the land cover types (not including water) and then combined these rasters into a stack, which is an object containing multiple rasters that have the same extent, spatial resolution, and thematic content.

Developed <- nlcd16_rc == 2
Barren <- nlcd16_rc == 3
Forest <- nlcd16_rc == 4
GrassShrub <- nlcd16_rc == 5
Cropland <- nlcd16_rc == 6
Wetland <- nlcd16_rc == 7

nlcd16_stk <- stack(Developed, Barren, Forest, GrassShrub, Cropland, Wetland)
names(nlcd16_stk) <- newnames[-1]
nlcd16_stk
## class      : RasterStack 
## dimensions : 1501, 1708, 2563708, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 1098705, 1149945, 1238565, 1283595  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## names      : Developed, Barren, Forest, GrassShrub, Cropland, Wetland 
## min values :         0,      0,      0,          0,        0,       0 
## max values :         1,      1,      1,          1,        1,       1

Conveniently, the rasterdf() function works on raster stacks as well as individual raster layers. In the resulting data frame, the “variable” column contains the name of each raster layer in the stack. The facet_wrap() function can be used to display each of these layers as a separate map.

nlcd16_stk_df <- rasterdf(nlcd16_stk)
ggplot(data = nlcd16_stk_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) +
  scale_fill_manual(name = "Present", 
                    values = c("gray80", "gray20"),
                    labels = c("No", "Yes"),
                    na.translate = FALSE) +  
  facet_wrap(~ variable, ncol = 3) +
  coord_sf(expand = F) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

7.3 Focal analysis of raster data

In situations where land cover types are relative rare and are distributed in small patches, they can be hard to see in the categorical map as well as the small multiple map. One way to visualize the concentrations of these classes to use a focal analysis to summarize the density of pixels within a larger circular window. First a weights object is defined using the focalWeight() function. Here, all cells within 500 m of the focal cell are included in the summary. The weight of each cell within this window is equal to its proportion of the total window size. Then the focal() function is called with the the weights objects and a summary function as arguments. By summing across the weights, the focal() function computes the proportion of each land cover class within a 500 m radius.

fwts500c <- focalWeight(Developed, d=500, type = "circle")

Developed_f <- focal(Developed, w=fwts500c, fun=sum)
Barren_f <- focal(Barren, w=fwts500c, fun=sum)
Forest_f <- focal(Forest, w=fwts500c, fun=sum)
GrassShrub_f <- focal(GrassShrub, w=fwts500c, fun=sum)
Cropland_f <- focal(Cropland, w=fwts500c, fun=sum)
Wetland_f <- focal(Wetland, w=fwts500c, fun=sum)

focal_stk <- stack(Developed_f, Barren_f, Forest_f, GrassShrub_f, 
                   Cropland_f, Wetland_f)
names(focal_stk) <- names(nlcd16_stk)

focal_stk_df <- rasterdf(focal_stk)
ggplot(data = focal_stk_df) +
  geom_raster(aes(x = x, y = y, fill = value)) +
  scale_fill_gradientn(name = "Density", colors=brewer.pal(9, "YlOrRd")) +
  facet_wrap(~ variable, ncol = 3) +
  coord_sf(expand = F) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

Changing the size of the focal window changes the degree of smoothing in the resulting maps, similar to the bandwidth of a kernel density analysis.

fwts100c <- focalWeight(Developed, d=100, type = "circle")
fwts500c <- focalWeight(Developed, d=500, type = "circle")
fwts1000c <- focalWeight(Developed, d=1000, type = "circle")
fwts2000c <- focalWeight(Developed, d=2000, type = "circle")

dev_100 <- focal(Developed, w=fwts100c, fun=sum)
dev_500 <- focal(Developed, w=fwts500c, fun=sum)
dev_1000 <- focal(Developed, w=fwts1000c, fun=sum)
dev_2000 <- focal(Developed, w=fwts2000c, fun=sum)

focal_stk2 <- stack(dev_100, dev_500, dev_1000, dev_2000)
names(focal_stk2) <- c("100m", "500m", "1000m", "2000m")

focal_stk2_df <- rasterdf(focal_stk2)
ggplot(data = focal_stk2_df) +
  geom_raster(aes(x = x, y = y, fill = value)) +
  scale_fill_gradientn(name = "Density", colors=brewer.pal(9, "YlOrRd")) +
  facet_wrap(~ variable, ncol = 2) +
  coord_sf(expand = F) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

By default, the focal() function will return a value of NA if there are any NA values inside the window, or if the window extends outside the boundary of the raster dataset. These NA values are displayed in gray in the preceding figure, and the gray boundary around the outside of the dataset increases in width as the window size increases.

7.4 Land cover change analysis

Often, we are interested in understanding how land cover has changed over time as a result of urban growth, agricultural abandonment, clearcutting, and other land use activities. To explore land cover trends in Walton County, GA, we will import six additional years of NLCD data and combine with with the 2016 NLCD dataset into a raster stack. Then the entire stack is reclassified into the seven broader land cover classes using the subs() function.

nlcd13 <- raster("NLCD_2013_Land_Cover_WaltonGA.tiff")
nlcd11 <- raster("NLCD_2011_Land_Cover_WaltonGA.tiff")
nlcd08 <- raster("NLCD_2008_Land_Cover_WaltonGA.tiff")
nlcd06 <- raster("NLCD_2006_Land_Cover_WaltonGA.tiff")
nlcd04 <- raster("NLCD_2004_Land_Cover_WaltonGA.tiff")
nlcd01 <- raster("NLCD_2001_Land_Cover_WaltonGA.tiff")

nlcd_stk <- stack(nlcd01, nlcd04, nlcd06, nlcd08,
                  nlcd11, nlcd13, nlcd16)
nlcd_stk_rc <- subs(nlcd_stk, lookup)
names(nlcd_stk_rc) <- c("Land Cover 2001", 
                        "Land Cover 2004",
                        "Land Cover 2006",
                        "Land Cover 2008",
                        "Land Cover 2011",
                        "Land Cover 2013",
                        "Land Cover 2016")

The ggplot() function can be used with facet_wrap() to map the time series of land cover.

nlcd_stk_df  <- rasterdf(nlcd_stk_rc)

ggplot(data = nlcd_stk_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = newcols2,
                    labels = newnames,
                    na.translate = FALSE) +
  coord_sf(expand = F) +
  facet_wrap(~ variable, ncol = 4) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position="bottom", 
        panel.background = element_rect(fill = "white", color = "black"))

If you look closely at these maps, it is possible to see areas where land cover is changing. However, it is difficult to get a clear sense of the overall trends - which land cover classes are changing and by how much? It would be helpful to calculate the total area of each land cover class and plot the changes over time. However, the rectangular boundary of our map includes Walton County as well as other surrounding areas that fit within the rectange. Restricting these summaries to Walton County requires the following steps. 1. Read in a shapefile of Georgia counties. 2. Reproject the shapefile into the same coordinate system as the raster data. 3. Extract Walton county from the shapefile. 4. Use the mask() function to convert all pixels outside of Walton County to NA values. 5. Use the crop() function to trim down the size of the raster to the borders of Walton County.

gacounty <- st_read("GA_SHP.shp")
## Reading layer `GA_SHP' from data source `E:\work\projects\gdswr-book\GA_SHP.shp' using driver `ESRI Shapefile'
## Simple feature collection with 159 features and 18 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -85.60516 ymin: 30.35576 xmax: -80.75143 ymax: 35.00066
## geographic CRS: NAD83
gacounty <- st_transform(gacounty, crs(nlcd_stk_rc))
walton <- filter(gacounty, NAME10 == "Walton")

nlcd_stk_msk <- mask(nlcd_stk_rc, walton)
nlcd_stk_crp <- crop(nlcd_stk_msk, walton)

nlcd_crp_df <- rasterdf(nlcd_stk_crp)

ggplot(data = nlcd_crp_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = newcols2,
                    labels = newnames,
                    na.translate = FALSE) +
  coord_sf(expand = F) +
  facet_wrap(~ variable, ncol = 4) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position="bottom", 
        panel.background = element_rect(fill = "white", color = "black"))

The freq() function can be used to extract the the number of cells of each land cover class from each layer in the raster stack. The output can be converted to a data frame with land cover classes as rows and years as columns. Before plotting the data frame, a series of dplyr functions is used to select only the columns of interest, rename them, convert to a long format with one row for each combination of year and land cover class, and convert from cell county to square kilometers.

freq_list <- freq(nlcd_stk_crp)
freq_df <- as.data.frame(freq_list)

nlcd_chg <- freq_df %>%
  dplyr::select(c(1, seq(2,14,2))) %>%
  rename("Class" = 1,
         "Y2001" = 2,
         "Y2004" = 3,
         "Y2006" = 4,
         "Y2008" = 5,
         "Y2011" = 6,
         "Y2013" = 7,
         "Y2016" = 8) %>%
  pivot_longer(cols = Y2001:Y2016, 
               names_to = "yearchar", 
               values_to = "cells") %>%
  mutate(km2 = cells * 900 / 1000000,
         year = parse_number(yearchar)) %>%
  drop_na()

In the preceding code, note that the dplyr select() function was called explicitly as dplyr::select() to avoid confusion with the select() function in the raster package.

The following code generates a line graph of change in area for each land cover class.

ggplot(data = nlcd_chg) +
  geom_line(aes(x = year, y = km2, color = as.factor(Class))) +
  geom_point(aes(x = year, y = km2, color = as.factor(Class))) +
  scale_color_manual(name = "Land Cover Class",
                     values = newcols, 
                     labels = newnames) +
  labs(x = "Year", y = expression(Area~(km^2))) +
  theme_classic()

Another way to graph the changes is with a separate subplot for each land cover class. Before generating this graph, Class is convert to a factor with labels corresponding to the class names. This approach allows ensures that Class is treated like a categorical variable instead of a numeric variable, and displays the label for each factor level in the strip above each facet.

nlcd_chg <- nlcd_chg %>%
  mutate(Class = factor(Class,
                        levels = 1:7,
                        labels = newnames))

ggplot(data = nlcd_chg) +
  geom_line(aes(x = year, y = km2)) +
  facet_wrap(~ Class, ncol = 4) +
  labs(x = "Year", y = expression(Area~(km^2))) +
  theme_bw()

In both of the previous graphs, it is very difficult to see the trends of the less common land cover classes. To see their change more clearly, the scale of the y-axis can be changed to vary freely with the range of values for each class.

ggplot(data = nlcd_chg) +
  geom_line(aes(x = year, y = km2)) +
  facet_wrap(~ Class, scales = "free_y", ncol = 4) +
  labs(x = "Year", y = expression(Area~(km^2))) +
  theme_bw()

7.5 Land cover transition matrices

Another way to quantify land cover change between two time periods is through a transition matrix. This is an n x n matrix, where n is the number of land cover classes. In this example, the 2001 classes are displayed as rows and the 2011 classes are displayed as columns. So for a particular class in 2001 (row) we can look across the numbers in each column to see how many pixels of that class have transitioned to the other land cover classes. To generate the transition matrix we used the crosstab()function to calculate the values. Then we rename the columns, remove rows containing NA values, and convert the data frame into matrix format.

changemat <- crosstab(nlcd_stk_crp[["Land.Cover.2001"]], 
                      nlcd_stk_crp[["Land.Cover.2016"]])
class(changemat)
## [1] "xtabs" "table"
names(changemat)
## NULL
rownames(changemat) <- newnames
colnames(changemat) <- newnames
changemat
##                Land.Cover.2016
## Land.Cover.2001  Water Developed Barren Forest GrassShrub Cropland Wetland
##      Water        9221        31     18    155         94        0     114
##      Developed       0    111969      0      0          0        0       0
##      Barren         46      1357   3542     37         29        0       3
##      Forest       1826     14817     76 428736      23331        6     110
##      GrassShrub    181     10171    657  35163     256227      670      59
##      Cropland        0        90      0      0          0     2522       0
##      Wetland       580        23      7     21         35        0   47729

When interpreting a change matrix, note that the numbers along the diagonal from upper left to lower right represent cells that did not change, including water cells that remained water, developed cells that remained developed, etc. Calculating the total of each row will provide the numbers of cells in each class in 2001, and the total of each column will provide the numbers of cells in each class in 2016.

It is also useful to express the transition matrix as the percentage of each class in 2001 that changes to each class type in 2011. A quick and easy way to do this is to use a couple of matrix manipulation functions from base R. The apply() function allows us to calculate the total number of pixels in each land cover class in 2001 by summing across the tows. The sweep() function allows us to “sweep” this vector across the array such that each value in the array is divided by its corresponding row sum.

sum01 <- apply(changemat, MARGIN=1, FUN=sum)
percchange <- 100 * sweep(changemat, MARGIN=1, STATS=sum01, FUN="/")
options(scipen = 999)
print(percchange, digits = 1)
##                Land.Cover.2016
## Land.Cover.2001   Water Developed  Barren  Forest GrassShrub Cropland Wetland
##      Water       95.723     0.322   0.187   1.609      0.976    0.000   1.183
##      Developed    0.000   100.000   0.000   0.000      0.000    0.000   0.000
##      Barren       0.917    27.064  70.642   0.738      0.578    0.000   0.060
##      Forest       0.389     3.160   0.016  91.434      4.976    0.001   0.023
##      GrassShrub   0.060     3.355   0.217  11.600     84.528    0.221   0.019
##      Cropland     0.000     3.446   0.000   0.000      0.000   96.554   0.000
##      Wetland      1.198     0.048   0.014   0.043      0.072    0.000  98.624

Converting the change matrix into a data frame automatically transforms it into long format. Then the change matrix can be plotted as a grouped barchart, in which each group represents a 2001 land cover class, and each bar within a group represents a 2016 land cover class.

changedf <- data.frame(changemat)

changedf <- changedf %>%
  mutate(ha = Freq * 900/1000000) %>%
  rename(lc2001 = 1,
         lc2016 = 2)

ggplot(data = changedf) +
  geom_bar(aes(x = lc2001, y = ha, group = lc2016, fill = lc2016),
           color = "black",
           position = "dodge",
           stat = "identity") +
  scale_fill_manual(name = "2016 Land Cover", values = newcols) +
  labs(x = "2001 Land Cover", y = expression(Area~(km^2))) 

Plotting change as percent area (relative the the 2001 area of each land cover class) makes it easier to see the changes occuring in the less abundance classes.

percchangedf <- data.frame(percchange)

percchangedf <- percchangedf %>%
  rename(lc2001 = 1,
         lc2016 = 2)

ggplot(data = percchangedf) +
  geom_bar(aes(x = lc2001, y = Freq, group = lc2016, fill = lc2016),
           color = "black",
           position = "dodge",
           stat = "identity") +
  scale_fill_manual(name = "2016 Land Cover", values = newcols) +
  labs(x = "2001 Land Cover", y = "% Area") 

7.6 Mapping specific types of land cover change

To better see where specific types of change are occurring, new classified data characterizing change events can be generated. The example will focus on changes in forest cover. Forest loss is classified as pixels with forest in 2001 and no forest in 2011, and forest gain is classified as pixels with no forest in 2001 and forest in 2011. Raster algebra is used to create a new raster layer where 1 = no forest change, 2 = forest loss, and 3 = forest gain.

nlcd01_rc <- nlcd_stk_rc[["Land.Cover.2001"]]
nlcd16_rc <- nlcd_stk_rc[["Land.Cover.2016"]]
forloss <- nlcd01_rc == 4 & nlcd16_rc != 4
forgain <- nlcd01_rc != 4 & nlcd16_rc == 4
forchange <- 1 + forloss + forgain * 2
forchange_df <- rasterdf(forchange)

Now these forest changes can be mapped.

names_chg <- c("No Change", "Forest Loss", "Forest Gain")
cols_chg <- c("lightgrey", "red", "blue")

ggplot(forchange_df) +
  geom_raster(aes(x = x, y = y, fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = cols_chg,
                    labels = names_chg,
                    na.translate = FALSE) +
  coord_sf(expand = F) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white", color = "black"))

7.7 Practice

  1. Conduct focal analyses of forest land cover using 100m, 500m, 1000m, and 2000m circular window sizes. Make a stack of the resulting raster layers and display them as a faceted map.

  2. Generate two change Walton County change matrices, one for 2001-2008 and one from 2008-2016. Use these outputs to determine which time period had a greater amount of land cover change.

  3. Generate a classified change map highlighting areas of gain and loss of the Grass/Shrub land cover class.