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(tmap)
library(sf)

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

Get Elevation Data for Thurston County

To get elevation data, we’ll use the USDA site. It’s fairly easy to use. I had some problems with Chrome, but Safari worked smoothly.

Go to https://datagateway.nrcs.usda.gov and click on the Get Data Button.

This video will show you how to get elevation data. You’ll need to copy this link to your browser.

https://www.youtube.com/watch?v=5zh3VVO_ebA

Get a tif

The four tif files I downloaded from the USDA site are available in Moodle as T1.tif, etc.

Download all of these and create rasters T1,T2,T3,T4 using the terra package.

Solution

T1 = rast("T1.tif")
T2 = rast("T2.tif")
T3 = rast("T3.tif")
T4 = rast("T4.tif")

Create Thurston

Use terra::merge to combine the four rasters into a single raster called Thurston.

Solution

Thurston = terra::merge(T1,T2,T3,T4)

Examine Thurston

Use the methods presented in section 6.1 to examine the characteristics of Thurston.

Solution

Thurston

## class       : SpatRaster 
## dimensions  : 7432, 5168, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 422485.1, 577525.1, 5093929, 5316889  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / UTM zone 10N (EPSG:26910) 
## source(s)   : memory
## name        :         T1 
## min value   :  -10.39135 
## max value   : 2543.97705

summary(Thurston)

## Warning: [summary] used a sample

##        T1          
##  Min.   :  -1.806  
##  1st Qu.:  72.064  
##  Median : 158.066  
##  Mean   : 325.648  
##  3rd Qu.: 459.016  
##  Max.   :2451.856  
##  NA's   :2260

ncell(Thurston)

## [1] 38408576

Tmap

Use tmap to see Thurston. Play with the breaks in tm_raster() until you get something you like.

Solution

breaks = c(-20,20,50,75,100,200,500,2500)
tm_shape(Thurston) +
  tm_raster(breaks = breaks,midpoint = NA)

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

## Warning: Values have found that are higher than the highest break

Histogram

Create a histogram of the data.

Solution

hist(Thurston)

## Warning: [hist] a sample of3% of the cells was used (of which 2% was NA)

The data has a strong right skew, which can be problematic for graphics. A lot of the data has been forced into a fairly small area to the left in order to accomodate the rare data much farther to the right.

Using a square root transformation can be useful in this situation.

Transformation.

Create a square root transformation of the data in Thurston. Note that this will make values less than or equal to zero NA.

Solution

SQThurston = sqrt(Thurston)
summary(SQThurston)

## Warning: [summary] used a sample

##        T1        
##  Min.   : 0.000  
##  1st Qu.: 8.515  
##  Median :12.584  
##  Mean   :15.054  
##  3rd Qu.:21.446  
##  Max.   :49.516  
##  NA's   :2415

class(SQThurston)

## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"

hist(SQThurston)

## Warning: [hist] a sample of3% of the cells was used (of which 2% was NA)

Note that there is more visible detail in the histogram of the transformed data.

An Alternative Graphic

Instead of specifying breaks, use style = “cont”.

Do this for both the raw Thurston and the transformed SQThurston.

Compare these two side-by-side using tmap_arrange.

Solution

ml = tm_shape(SQThurston) +
  tm_raster(style = "cont",
            palette = "PuRd") +
  tm_layout(legend.show = F,
            title = "SQRT")

m = tm_shape(Thurston) +
  tm_raster(style = "cont",
            palette = "PuRd") +
  tm_layout(legend.show = F,
            title = "Raw") 

tmap_arrange(m,ml)

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

See Thurston

Load the file “thurston_sf.Rdata”. It contains the boundary data for Thurston County. Use it to overlay the raster data.

Solution

load("thurston_sf.Rdata")

tm_shape(SQThurston) +
  tm_raster(style = "cont",
            palette = "PuRd") +
  tm_shape(thurston_sf) +
     tm_borders() +
  tm_layout(legend.show = F,
            title = "SQRT")

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

Thurston Itself

Crop the raster data to the boundaries of Thurston County.

Solution

thurston_sf = st_transform(thurston_sf,crs(SQThurston))
cropped = crop(SQThurston,thurston_sf)

tm_shape(cropped) +
  tm_raster(style = "cont",
            palette = "PuRd") +
  tm_shape(thurston_sf) + 
  tm_borders(lwd = 3) +
  tm_layout(legend.show = F)

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

Section 6.3

Run all of the data-building code in this section.

Solution

# Read the yearly temperature data.

temp2012 <- rast("NLDAS_FORA0125_M.A201207.002.grb.temp.tif")
temp2013 <- rast("NLDAS_FORA0125_M.A201307.002.grb.temp.tif")
temp2014 <- rast("NLDAS_FORA0125_M.A201407.002.grb.temp.tif")
temp2015 <- rast("NLDAS_FORA0125_M.A201507.002.grb.temp.tif")
temp2016 <- rast("NLDAS_FORA0125_M.A201607.002.grb.temp.tif")
temp2017 <- rast("NLDAS_FORA0125_M.A201707.002.grb.temp.tif")

# Build the stack.

tempstack <- c(temp2012, temp2013, temp2014, 
               temp2015, temp2016, temp2017)
names(tempstack) <- c("July 2012", "July 2013", "July 2014", 
                      "July 2015", "July 2016", "July 2017")

# Convert the data from Kelvin to Celsius.

tempstack <- tempstack - 273.15

Make the Graphic

Use tmap to display a graphic similar to what the author did with ggplot2.

Solution

tm_shape(tempstack) +
  tm_raster(legend.is.portrait = F,
            style = "cont",
            palette = "PuRd") +
  tm_grid() + 
  tm_layout(title = "Mean July Temperature",
            legend.outside = T,
            legend.outside.position = "bottom")