Part 1 of Chapter 7

This includes 7.1 and 7.2

Setup

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.1     ✔ purrr   0.3.4
## ✔ tibble  3.2.1     ✔ dplyr   1.1.1
## ✔ tidyr   1.2.1     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(terra)

## terra 1.7.3
## 
## Attaching package: 'terra'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(sf)

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(tmap)

Get Data

Get the land cover data for Washington state downloaded from https://www.mrlc.gov/viewer/. Just take the 2001 data.

wa_lcd = rast("data/data_2001.tiff")

Look

Use tmap to take a look at what we have.

Solution

load("thurston_sf.Rdata")
tm_shape(wa_lcd) +
  tm_raster() +
tm_shape(thurston_sf) +
  tm_borders(lwd = 3, col = "green") +
tm_layout(legend.show = F)

## stars object downsampled to 1039 by 963 cells. See tm_shape manual (argument raster.downsample)

## Warning: Duplicated levels found. They have been omitted

Crop to Thurston County

thurston_sf = st_transform(thurston_sf,crs(wa_lcd))
thurston_rast = crop(wa_lcd,thurston_sf)
tm_shape(thurston_rast) +
  tm_raster() +
tm_shape(thurston_sf) +
  tm_borders(lwd = 3,col = "green") +
tm_layout(legend.show = F)

## stars object downsampled to 1076 by 929 cells. See tm_shape manual (argument raster.downsample)

Use ggplot2

First use rasterdf.R to convert the raster to a dataframe.

Solution

source("rasterdf.R")
thurston_df <- rasterdf(thurston_rast)

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

Get Land Cover Codes

LCcodes <- unique(thurston_rast)[, 1]
LCcodes

##  [1] 11 21 22 23 24 31 41 42 43 52 71 81 82 90 95

Copy LCnames

Also get the official colors matching these names.

Solution

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

nlcdcols <- data.frame(coltab(thurston_rast))
nlcdcols <- nlcdcols[LCcodes + 1,]
LCcolors <- rgb(red = nlcdcols$red,
                green = nlcdcols$green,
                blue = nlcdcols$blue,
                names = as.character(LCcodes),
                maxColorValue = 255)
LCcolors

##        11        21        22        23        24        31        41        42 
## "#466B9F" "#DEC5C5" "#D99282" "#EB0000" "#AB0000" "#B3AC9F" "#68AB5F" "#1C5F2C" 
##        43        52        71        81        82        90        95 
## "#B5C58F" "#CCB879" "#DFDFC2" "#DCD939" "#AB6C28" "#B8D9EB" "#6C9FB8"

ggplot2 with Proper Colors

Plot the data with these colors.

Solution

ggplot(data = thurston_df) +
  geom_raster(aes(x = x, 
                  y = y, 
                  fill = as.character(value))) + 
  scale_fill_manual(name = "Land cover",
                    values = LCcolors,
                    labels = LCnames,
                    na.translate = FALSE) +
  coord_sf(expand = FALSE) +
  theme_void()

Section 7.2

In this section, we reduce the number of land cover types. This will reduce the visual clutter.

First decide how to collapse the existing codes into a smaller set and build a table we can use to convert the raster data. Follow the author of our text. We also assign names to the collapsed categories and assign an appropriate color to each.

Solution

LCcodes

##  [1] 11 21 22 23 24 31 41 42 43 52 71 81 82 90 95

##  [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(LCcodes, newclas)

newnames <- c("Water",
               "Developed",
               "Barren",
               "Forest",
               "GrassShrub",
               "Cropland",
               "Wetland")
newcols <- c("mediumblue", 
             "firebrick2", 
             "gray60", 
             "darkgreen", 
             "yellow2", 
             "orange4", 
             "paleturquoise2")

Now we apply the table to thurston_rast and produce thurston_rc. Note that the function classify is from the terra package.

thurston_rc <- classify(thurston_rast, lookup)

PLot

Now plot the raster with the collapsed values.

Solution

thurston_rc_df <- rasterdf(thurston_rc)

ggplot(data = thurston_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 = FALSE) +
  theme_void()

A Stack of Binaries

Convert our data to a stack of binaries and then to a dataframe. Don’t eliminate water. Look at the dataframe with str().

Solution

thurston_stk <- segregate(thurston_rc)
names(thurston_stk) = newnames
thurston_stk_df = rasterdf(thurston_stk)
str(thurston_stk_df)

## 'data.frame':    33424167 obs. of  4 variables:
##  $ x       : num  -2054490 -2054460 -2054430 -2054400 -2054370 ...
##  $ y       : num  2977650 2977650 2977650 2977650 2977650 ...
##  $ value   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ variable: Factor w/ 7 levels "Water","Developed",..: 1 1 1 1 1 1 1 1 1 1 ...

Distribution of LC Codes

Use group_by() and summarize to see this distribution.

Solution

thurston_stk_df %>% 
  group_by(variable) %>% 
  summarise(mean = mean(value)) 

## # A tibble: 7 × 2
##   variable       mean
##   <fct>         <dbl>
## 1 Water      0.0712  
## 2 Developed  0.189   
## 3 Barren     0.00569 
## 4 Forest     0.433   
## 5 GrassShrub 0.241   
## 6 Cropland   0.000965
## 7 Wetland    0.0591

Note that this describes the entire area in our map, which goes beyond Thurston County.

Getting a County

Load the counties_us.Rdata file and look at the first few records.

load("counties_us.Rdata")
head(counties_us)

## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -97.97959 ymin: 36.02183 xmax: -76.52951 ymax: 45.93594
## Geodetic CRS:  WGS 84
##   STATEFP COUNTYFP COUNTYNS       AFFGEOID GEOID      NAME LSAD      ALAND
## 1      24      510 01702381 0500000US24510 24510 Baltimore   25  209650970
## 2      31      169 00835906 0500000US31169 31169    Thayer   06 1486153245
## 3      37      077 01008556 0500000US37077 37077 Granville   06 1377850155
## 4      46      091 01265759 0500000US46091 46091  Marshall   06 2170724143
## 5      39      075 01074050 0500000US39075 39075    Holmes   06 1094405864
## 6      39      033 01074029 0500000US39033 39033  Crawford   06 1040619651
##      AWATER                       geometry
## 1  28758714 MULTIPOLYGON (((-76.71131 3...
## 2   3025339 MULTIPOLYGON (((-97.82082 4...
## 3  14503653 MULTIPOLYGON (((-78.80252 3...
## 4 124586482 MULTIPOLYGON (((-97.97924 4...
## 5   3695230 MULTIPOLYGON (((-82.22066 4...
## 6   2359514 MULTIPOLYGON (((-83.11274 4...

Go to https://transition.fcc.gov/oet/info/maps/census/fips/fips.txt to determine the FIPS code for the state you’re interested in.

I found that the code for Washington is 53.

Loak at the records in counties_us with COUNTYFP equal “53”. You need to drop the geometry temporarily.

wa = counties_us %>% 
  st_drop_geometry() %>% 
  filter(STATEFP == "53") %>%
  select(STATEFP,COUNTYFP,GEOID,NAME) %>% 
  arrange(NAME)
wa

##    STATEFP COUNTYFP GEOID         NAME
## 1       53      001 53001        Adams
## 2       53      003 53003       Asotin
## 3       53      005 53005       Benton
## 4       53      007 53007       Chelan
## 5       53      009 53009      Clallam
## 6       53      011 53011        Clark
## 7       53      013 53013     Columbia
## 8       53      015 53015      Cowlitz
## 9       53      017 53017      Douglas
## 10      53      019 53019        Ferry
## 11      53      021 53021     Franklin
## 12      53      023 53023     Garfield
## 13      53      025 53025        Grant
## 14      53      027 53027 Grays Harbor
## 15      53      029 53029       Island
## 16      53      031 53031    Jefferson
## 17      53      033 53033         King
## 18      53      035 53035       Kitsap
## 19      53      037 53037     Kittitas
## 20      53      039 53039    Klickitat
## 21      53      041 53041        Lewis
## 22      53      043 53043      Lincoln
## 23      53      045 53045        Mason
## 24      53      047 53047     Okanogan
## 25      53      049 53049      Pacific
## 26      53      051 53051 Pend Oreille
## 27      53      053 53053       Pierce
## 28      53      055 53055     San Juan
## 29      53      057 53057       Skagit
## 30      53      059 53059     Skamania
## 31      53      061 53061    Snohomish
## 32      53      063 53063      Spokane
## 33      53      065 53065      Stevens
## 34      53      067 53067     Thurston
## 35      53      069 53069    Wahkiakum
## 36      53      071 53071  Walla Walla
## 37      53      073 53073      Whatcom
## 38      53      075 53075      Whitman
## 39      53      077 53077       Yakima

Let’s get Cowlitz County and restore its geometry. Finally, get the bounding box. This will allow you to mark the correct area in the MRLC Viewer.

cowlitz_sf = counties_us %>% 
  st_drop_geometry() %>% 
  filter(GEOID == "53015") %>% 
  left_join(counties_us )

## Joining with `by = join_by(STATEFP, COUNTYFP, COUNTYNS, AFFGEOID, GEOID, NAME,
## LSAD, ALAND, AWATER)`

cowlitz_sf = st_as_sf(cowlitz_sf)
st_bbox(cowlitz_sf)

##       xmin       ymin       xmax       ymax 
## -123.21795   45.85054 -122.23991   46.38804

Comparison Over Time

First, we need to get the data for several years, reclassify it as before and assign names to the rasters in the stack.

thurston_stk <- rast(c("data/data_2001.tiff",
                   "data/data_2004.tiff",
                   "data/data_2006.tiff",
                   "data/data_2008.tiff",
                   "data/data_2011.tiff",
                   "data/data_2013.tiff",
                   "data/data_2016.tiff",
                   "data/data_2019.tiff"))

thurston_rc <- classify(thurston_stk, lookup)
names(thurston_rc) <- c("2001", "2004", "2006", "2008", "2011", "2013", "2016", "2019")

A Graphic

I ran the code in the text exactly as I found it after fixing one type (newcols2 instead of newcols).

I couldn’t see any real differences, so I trimmed the list down to just 2001 and 2019. I still don’t notice a difference.

thurston_rc_df  <- rasterdf(thurston_rc[[c("2001","2019")]] )

ggplot(data = thurston_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 = TRUE) +
  facet_wrap(facets = vars(variable), ncol = 2) +
  theme_void() +
  theme(strip.text.x = element_text(size = 12, face="bold"),
        legend.position="bottom")

Another Approach

I decided to produce the time-series graph showing all of the coverage types in one panel.

I copied the code from the text and replaced “NLCD” with thurston everywhere.

freq_df <- freq(thurston_rc, usenames=TRUE)
glimpse(freq_df)

## Rows: 56
## Columns: 3
## $ layer <chr> "2001", "2001", "2001", "2001", "2001", "2001", "2001", "2004", …
## $ value <dbl> 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1…
## $ count <dbl> 793735, 1871085, 53422, 4134774, 2063068, 12571, 500985, 788387,…

thurston_chg <- freq_df %>%
  mutate(km2 = count * 900 / 1000000,
         class = factor(value,
                        levels = 1:7,
                        labels = newnames),
         year = as.numeric(layer))
head(freq_df)

##   layer value   count
## 1  2001     1  793735
## 2  2001     2 1871085
## 3  2001     3   53422
## 4  2001     4 4134774
## 5  2001     5 2063068
## 6  2001     6   12571

ggplot(data = thurston_chg) +
  geom_line(aes(x = year, y = km2, color = class)) +
  geom_point(aes(x = year, y = km2, color = class),size=3) +
  scale_color_manual(name = "Land Cover Class",
                     values = newcols) +
  labs(x = "Year", y = expression("Area(km"^2*")")) +
  theme_classic()

This is a very interesting graph. There is clearly a relationship betwee grass/shrub and forest. The two curves are essentialy mirror images. At first, forest declines and grass/shrub grows. After 2010, the opposite is true. How would you explain this?

The best explanation is that the early time period involved some clear-cut logging. Apparently, the forest regrew, probably with the help of a reseeding program.

Comparing 2001 and 2019

I want to do this numerically. I need to extract the data for these two year.

thurston_ends = thurston_chg %>% 
  filter(layer == "2001" | layer == "2019") %>% 
  select(layer, km2, class) %>% 
  arrange(class,layer)

head(thurston_ends)

##   layer       km2     class
## 1  2001  714.3615     Water
## 2  2019  715.7817     Water
## 3  2001 1683.9765 Developed
## 4  2019 1794.4515 Developed
## 5  2001   48.0798    Barren
## 6  2019   36.1125    Barren

Pivot

We want to do arithmetic with the 2001 and 2019 values of each class. To do this we need to pivot the data.

thurston_comp = thurston_ends %>% 
  pivot_wider(names_from = layer, 
              names_prefix = "km2_",
              values_from = km2)
thurston_comp

## # A tibble: 7 × 3
##   class      km2_2001 km2_2019
##   <fct>         <dbl>    <dbl>
## 1 Water         714.     716. 
## 2 Developed    1684.    1794. 
## 3 Barren         48.1     36.1
## 4 Forest       3721.    4077. 
## 5 GrassShrub   1857.    1393. 
## 6 Cropland       11.3     14.3
## 7 Wetland       451.     456.

Comparisons

Get the absolute and percentage changes.

thurston_comp = thurston_comp %>% 
  mutate(change = km2_2019 - km2_2001,
         pctchange = change/km2_2001)


thurston_comp %>% 
  select(class,change,pctchange)

## # A tibble: 7 × 3
##   class       change pctchange
##   <fct>        <dbl>     <dbl>
## 1 Water         1.42   0.00199
## 2 Developed   110.     0.0656 
## 3 Barren      -12.0   -0.249  
## 4 Forest      355.     0.0955 
## 5 GrassShrub -463.    -0.250  
## 6 Cropland      2.94   0.260  
## 7 Wetland       5.20   0.0115

Transitions

We can make good guesses at where area is coming from and going to. It looks like the lost GrassShrub goes mostly to Forest and then some to developed. However it would be good to verify this by creating a transitions matrix.

We can copy the code from the text and edit only for the name change.

changeras <- c(thurston_rc[["2001"]], 
               thurston_rc[["2019"]])
changeout <- crosstab(changeras)
class(changeout)

## [1] "xtabs" "table"

changedf <- as_tibble(changeout)
head(changedf)

## # A tibble: 6 × 3
##   X2001 X2019      n
##   <chr> <chr>  <int>
## 1 1     1     782828
## 2 2     1          0
## 3 3     1       5695
## 4 4     1       1010
## 5 5     1       2333
## 6 6     1         11

shortnames <- c("Wat", 
                "Dev", 
                "Bare", 
                "For", 
                "Grass", 
                "Crop", 
                "Wet")

changedf <- changedf %>%
  mutate(X2001 = factor(X2001,
                        levels = 1:7,
                        labels = shortnames),
         X2019 = factor(X2019,
                        levels = 1:7,
                        labels = shortnames),
         ha = n * 900/10000) %>%
  group_by(X2001) %>%
  mutate(tot2001 = sum(ha),
         perc = 100 * ha / tot2001)

changedf

## # A tibble: 49 × 6
## # Groups:   X2001 [7]
##    X2001 X2019       n        ha tot2001     perc
##    <fct> <fct>   <int>     <dbl>   <dbl>    <dbl>
##  1 Wat   Wat    782828  70455.    71436.  98.6   
##  2 Dev   Wat         0      0    168398.   0     
##  3 Bare  Wat      5695    513.     4808.  10.7   
##  4 For   Wat      1010     90.9  372130.   0.0244
##  5 Grass Wat      2333    210.   185676.   0.113 
##  6 Crop  Wat        11      0.99   1131.   0.0875
##  7 Wet   Wat      3436    309.    45089.   0.686 
##  8 Wat   Dev       730     65.7   71436.   0.0920
##  9 Dev   Dev   1871048 168394.   168398. 100.    
## 10 Bare  Dev      2965    267.     4808.   5.55  
## # … with 39 more rows

The Matrix

The author used land area in hectares (ha). I prefer the count of 30X30 cells (n).

changemat <- matrix(changedf$n, 
                    nrow = 7, 
                    ncol = 7)
rownames(changemat) <- shortnames
colnames(changemat) <- shortnames

changemat

##          Wat     Dev  Bare     For  Grass  Crop    Wet
## Wat   782828     730  1152     978    987     0   7060
## Dev        0 1871048    24       6      3     4      0
## Bare    5695    2965 33987    1456   4448     7   4864
## For     1010   51295  3518 3498850 579839    38    224
## Grass   2333   63267  1422 1028232 962475  4688    651
## Crop      11    1284     2       2    169 11097      6
## Wet     3436    3246    20      83    242     0 493958