Part 1 of Chapter 7

This includes 7.1 and 7.2

Setup

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.4 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ 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, by = c("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