Introduction

About this R Notebook

R Notebooks are a ‘flavor’ of R markdown that combine plain text and R commands in code chunks. You can download the Rmd file from the ‘code’ button at the top of the page, then open it in RStudio. You run code chunks in the document line-by-line, and the output appears immediately below the code chunk. When you save the Rmd file, it automatically creates a HTML file showing results from all the code that has been run so far.

If you’re in RStudio, you can minimize the console window (and probably close the right-hand panes as well). You won’t need it, because when you run R commands in a R Notebook the output appears below the code chunk (not the console). This takes some getting used to.

Keyboard shortcuts:
- run the current line of R: ctrl + enter
- run everything in the current code chunk: ctrl + shift + enter
- insert a new code chunk: ctrl + alt + i

Setup

Install caladaptr:

if (!require(devtools)) {install.packages("devtools"); library(devtools)}

Loading required package: devtools
Loading required package: usethis
The working directory was changed to D:/GitHub/cal-adapt/caladaptr inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

# 
# devtools::install_github("ucanr-igis/caladaptr")
library(caladaptr)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
caladaptr (version 0.4.5)
URL: https://ucanr-igis.github.io/caladaptr
Bug reports: https://github.com/ucanr-igis/caladaptr/issues

## Make sure you have version 0.4.1 caladaptr or later. 
packageVersion("caladaptr")

[1] ‘0.4.5’

Now we can load the other packages we’ll be using below:

pkgs_req <- c("ggplot2", "dplyr", "tidyr", "lubridate", "conflicted")
pkgs_missing <- pkgs_req[!(pkgs_req %in% installed.packages()[,"Package"])]
if (length(pkgs_missing)) install.packages(pkgs_missing, dependencies=TRUE)

library(units)
library(ggplot2)
library(dplyr)
library(conflicted)
library(tidyr)
library(tmap)
library(lubridate)
library(sf)

The last setup task is to define your preferences when you use an ambiguous function name (i.e., a function that exists in more than one package). This is particularly needed with a few common generic functions from dplyr:

conflict_prefer("filter", "dplyr", quiet = TRUE)
conflict_prefer("count", "dplyr", quiet = TRUE)
conflict_prefer("select", "dplyr", quiet = TRUE)

Example #1: Get Projected Temperature for a Point

Goal: Make a time series plot of projected maximum annual temperature for a single point, using the four recommended GCMs for California under RCP 4.5.

1. Create the API Request

The first step in getting climate variables back is to create a Cal-Adapt API request object. This involves stringing together a series of functions that specify the pieces of the request. The following constants can help you specify pieces of the request:

cap1 <- ca_loc_pt(coords = c(-121.4687, 38.5938)) %>%
  ca_gcm(gcms[1:4]) %>%                                 
  ca_scenario(c("rcp45", "rcp85")) %>%
  ca_period("year") %>%
  ca_years(start = 2030, end = 2080) %>%
  ca_cvar(c("tasmax"))

Entering the name of an API request object at the console will display its components:

cap1

Cal-Adapt API Request
Location(s): 
  x: -121.469
  y: 38.594
Variable(s): tasmax
Temporal aggregration period(s): year
GCM(s): HadGEM2-ES, CNRM-CM5, CanESM2, MIROC5
Scenario(s): rcp45, rcp85
Dates: 2030-01-01 to 2080-12-31
 

Tip: You don’t have to memorize a bunch of keywords. caladaptr provides several constants that contain the values you can pass to API construction functions, including gcms, scenarios, cvars, and periods.

## Available global change models
## Note, the first 4 GCMs are the recommended priority models for California
gcms

 [1] "HadGEM2-ES" "CNRM-CM5"   "CanESM2"    "MIROC5"     "ACCESS1-0"  "CCSM4"      "CESM1-BGC"  "CMCC-CMS"  
 [9] "GFDL-CM3"   "HadGEM2-CC" "ens32avg"   "ens32max"   "ens32min"  

## Available emissions scenarios
scenarios

[1] "rcp45"      "rcp85"      "historical"

## Available climate variables
cvars

[1] "tasmax" "tasmin" "pr"     "swe"   

## Available temporal aggregation periods
periods

[1] "day"    "month"  "year"   "30yavg"

Note: Cal-Adapt has data for many but by no means all combinations of the above constants.

You can use these values as arguments to build up the API request:

To verify the location in an API request, you can plot it. (Note we still haven’t fetched any climate data yet, this just shows you the location the request will ask for.)

plot(cap1)

2. Fetch Data

Now it’s time to fetch data with ca_getvals_tbl(). The object returned by ca_getvals_tbl() is tibble (data frame):

cap1_tbl <- cap1 %>% 
  ca_getvals_tbl(quiet = TRUE)

cap1_tbl

3. Munge the Results

To produce the desired time series plot, we need to i) pull out just values for RCP 4.5, and ii) convert degrees to °F. For the unit conversion, we can use the handy set_units function from the units package.

cap1_rcp45_tbl <- cap1_tbl %>%
  filter(scenario == "rcp45") %>%
  mutate(temp_f = set_units(val, degF))

cap1_rcp45_tbl

Plot the Time Series

Plot these with ggplot:

ggplot(data = cap1_rcp45_tbl, aes(x = as.Date(dt), y = as.numeric(temp_f), group = gcm)) +
  geom_line(aes(color=gcm)) +
  labs(title = "Average Annual Maximum Temperature for RCP4.5", x = "year", y = "temp (F)")

YOUR TURN

Modify the above to create a similar plot for RCP 8.5.

## Plot of RCP 8.5 goes here

Example #2: Retrieve County Data

Goal: produce a time series plot showing the difference between RCP85 and RCP45 on maximum annual temperature for a single county.

The Cal-Adapt API has a number of ‘preset’ areas-of-interest (also called boundary layers) that you can use cookie-cutter style when retrieving climate data. The advantage of using an AOI Preset is that you don’t need to import a GIS layer to query according to these common features. You just need to know the name or ID number of the feature(s) you’re interested in.

The following AOI Presets are available:

aoipreset_types

 [1] "censustracts"      "counties"          "cdistricts"        "ccc4aregions"      "climregions"      
 [6] "hydrounits"        "irwm"              "electricutilities" "wecc-load-area"    "evtlocations"     
[11] "place"            

1. Find the FIPS code for your county of interest

To use an AOI Preset, you need to specify which feature(s) you’re interested in by providing the value(s) for one of the id fields. The specific columns available for identifying features vary according to preset. You can view the id columns and their values for an AOI Preset using the built-in aoipreset_idval constant. For example the counties layer allows you to specify a county by name, fips code, or id. Remember that everything in R is case-sensitive!

aoipreset_idval$counties

You can find county fips codes on Google - just be sure to use the 6-character version which includes the state. Alternately, you can plot the county preset layer and click on your county of interest:

counties_sf <- ca_aoipreset_geom("counties")

Reading layer `counties' from data source `C:\Users\Andy\AppData\Local\R\cache\R\caladaptr\counties.gpkg' using driver `GPKG'
Simple feature collection with 87 features and 54 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13871160 ymin: 3833648 xmax: -12625080 ymax: 5416187
projected CRS:  WGS 84 / Pseudo-Mercator

tmap_mode("view")
tm_shape(counties_sf) + tm_polygons()

For this example, we’ll look at Kings County (FIPS = 06031).

2. Create the API Request

Let’s create an API request for Kings County. Note below the inclusion of ca_options() to specify how we want to aggregate the pixels that fall within the country. This is required whenever you query polygon features.

cap2 <- ca_loc_aoipreset(type="counties", idfld = "fips", idval = "06031") %>%
  ca_gcm(gcms[1:4]) %>%
  ca_scenario(c("rcp45", "rcp85")) %>%
  ca_period("year") %>%
  ca_years(start = 2030, end = 2080) %>%
  ca_cvar(c("tasmin")) %>%
  ca_options(spatial_ag = "max")

cap2

Cal-Adapt API Request
Location(s): 
  AOI Preset: counties
  fips(s): 06031
Variable(s): tasmin
Temporal aggregration period(s): year
GCM(s): HadGEM2-ES, CNRM-CM5, CanESM2, MIROC5
Scenario(s): rcp45, rcp85
Dates: 2030-01-01 to 2080-12-31
Options:
  spatial ag: max
 

As before, we can plot the API request to double-check we got the right location:

plot(cap2)

3. Fetch Data

Fetch data with ca_getvals_tbl():

cap2_tbl <- cap2 %>% 
  ca_getvals_tbl(quiet = TRUE) 

cap2_tbl

4. Munge Results

To compute the difference between RCP 8.5 and RCP 4.5, we need to split them into separate columns. This is an example of going from a ‘long’ format to a ‘wide’ format. Fortunately, the tidyr package has a function called pivot_wider that can do this in single command. While we’re at it we’ll convert the temp to degrees Fahrenheit:

rcp85_minus_rcp45_tbl <- cap2_tbl %>% 
  mutate(temp_f = set_units(val, degF)) %>% 
  select(fips, gcm, scenario, dt, temp_f) %>% 
  pivot_wider(names_from = scenario, values_from = temp_f)

head(rcp85_minus_rcp45_tbl)

5. Plot

Now we’re ready to make the plot:

ggplot(data = rcp85_minus_rcp45_tbl, aes(x = as.Date(dt), y = as.numeric(rcp85 - rcp45), group = gcm)) +
  geom_line(aes(color=gcm)) +
  labs(title = "Difference between RCP8.5 and RCP4.5 in the Maximum Annual Minimum \nTemperature for Kings County", x = "year", y = "temp (F)")

Example #3: Analyze Extreme Heat in a Census Tract

Goal: Compare extreme heat occurrence under RCP 4.5 and RCP 8.5 for a single census tract over a 20-year time span. We’ll measure ‘extreme heat’ in two ways, i) the proportion of total days when the max temp exceeds a certain threshold, and ii) the number of ‘heat spells’ defined as a period of n or more consecutive days when the temperature exceeds a threshold.

1. Find the GEOID of a Census Tract in Kings County

Census tracts have 10-digit GEOID numbers rather than 6-digit FIPS codes (more info). To find the GEOID for a census tract in Kings County, we could:

1. Ask Google.
   
2. Plot all the census tracts (as we did above with counties), and click on the one we’re interested in.
   
3. Look it up using the tidycensus package.
   
4. Take advantage of the fact that the GEOID of census tracts in Kings County all start with its FIPS code (6031).

We’ll use method #4, and grab the census tracts GEOIDs from aoipreset_idval:

tracts_kings <- grep(pattern = "^6031", aoipreset_idval$censustracts$tract, value = TRUE)
tracts_kings

 [1] "6031000100" "6031000200" "6031000300" "6031000402" "6031000403" "6031000404" "6031000405" "6031000500"
 [9] "6031000601" "6031000602" "6031000701" "6031000702" "6031000800" "6031000900" "6031001001" "6031001002"
[17] "6031001003" "6031001100" "6031001200" "6031001300" "6031001401" "6031001402" "6031001500" "6031001601"
[25] "6031001602" "6031001701" "6031981800"

To plot just these tracts we can feed these GEOID values into filter():

census_tracts_sf <- ca_aoipreset_geom("censustracts")

Reading layer `censustracts' from data source `C:\Users\Andy\AppData\Local\R\cache\R\caladaptr\censustracts.gpkg' using driver `GPKG'
Simple feature collection with 8034 features and 63 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -13850030 ymin: 3833633 xmax: -12705030 ymax: 5162403
projected CRS:  WGS 84 / Pseudo-Mercator

tmap_mode("view")
tm_shape(census_tracts_sf %>% filter(tract %in% tracts_kings)) + tm_polygons()

For the rest of this exercise we’ll use 6031001701.

2. Construct the API Request

The following will give us 20 years of daily maximum temperature for our census tract. To keep it simple we’ll just use one GCM.

cap3 <- ca_loc_aoipreset(type="censustracts", idfld = "tract", idval = 6031001701) %>%
  ca_gcm("MIROC5") %>%
  ca_scenario(c("rcp45", "rcp85")) %>%
  ca_period("day") %>%
  ca_years(start = "2070-01-01", end = "2089-12-31") %>%
  ca_cvar("tasmax") %>% 
  ca_options(spatial_ag = "max")

cap3

Cal-Adapt API Request
Location(s): 
  AOI Preset: censustracts
  tract(s): 6031001701
Variable(s): tasmax
Temporal aggregration period(s): day
GCM(s): MIROC5
Scenario(s): rcp45, rcp85
Dates: 2070-01-01 to 2089-12-31
Options:
  spatial ag: max
 

plot(cap3)

3. Fetch Data

Now we’re ready to fetch data:

cap3_tbl <- cap3 %>% 
  ca_getvals_tbl(quiet = TRUE)

dim(cap3_tbl)

[1] 14610     8

head(cap3_tbl)

4. Munge the Results

The only ‘munging’ we need to do is to add a column for the temperature in Farenheit.

cap3_degf_tbl <- cap3_tbl %>% mutate(temp_f = set_units(val, degF))

head(cap3_degf_tbl)

5. Compare the Distribution of Maximum Daily Temp

Plot a histogram of maximum temperature over this 20-year period for each emission scenario:

hist(cap3_degf_tbl %>%  filter(scenario == "rcp45") %>% pull(temp_f),
     main = "RCP 4.5: Maximum Daily Temp, 2070-2089",
     xlab = "temp (F)") 

hist(cap3_degf_tbl %>%  filter(scenario == "rcp85") %>% pull(temp_f),
     main = "RCP 8.5: Maximum Daily Temp, 2070-2089",
     xlab = "temp (F)") 

6. Count the Number of Days over 110°F Projected Under each RCP

cap3_degf_tbl %>% 
  mutate(really_hot_tf = temp_f >= set_units(110, degF)) %>%
  group_by(scenario, really_hot_tf) %>% 
  count() %>% 
  pivot_wider(names_from = scenario, values_from = n)

7. Count the Number of Heat Spells Projected Under each RCP

A simple definition of ‘heat spell’ is consecutive hot days. Let’s count Count the number of heat spells consisting of 4 or more consecutive days of >110°F.

The first step is to pull out the values for just one RCP, sort them by date, add a logical column if the maximum temp that day exceeded our treshhold. We’ll then feed those values into rle(), which breaks up the time series into ‘runs’ of TRUE and FALSE.

## Use rle() to chop up a time series of TRUE of FALSE into 'runs'
rcp45_heat_runs <- cap3_degf_tbl %>% 
  filter(scenario == "rcp45") %>% 
  arrange(dt) %>% 
  mutate(really_hot_tf = temp_f >= set_units(110, degF)) %>% 
  pull(really_hot_tf) %>% 
  rle()

rcp45_heat_runs

Run Length Encoding
  lengths: int [1:241] 206 4 49 2 274 2 45 1 1 2 ...
  values : logi [1:241] FALSE TRUE FALSE TRUE FALSE TRUE ...

Next, we find the runs where two conditions are met: run length is >= 4 days, and the run value = TRUE (meaning its a run of really hot days).

sum(rcp45_heat_runs$values & rcp45_heat_runs$lengths >= 4)

[1] 36

We can do the same for RCP 8.5:

rcp85_heat_runs <- cap3_degf_tbl %>% 
  filter(scenario == "rcp85") %>% 
  arrange(dt) %>% 
  mutate(really_hot_tf = temp_f >= set_units(110, degF)) %>% 
  pull(really_hot_tf) %>% 
  rle()

sum(rcp85_heat_runs$values & rcp85_heat_runs$lengths >= 4)

[1] 45

Conclusion

Under RCP 4.5, the MIROC5 model predicts this census tract will experience 36 extreme heat spells (defined as 4 or more days where the temperature hit 110°F) during this 20 year period. In comparison, MIROC5 predicts 45 extreme heat spells under RCP 8.5.

YOUR TURN

How many heat spells will there be under each RCP if we define a heat spell to be at least 7 days over 100 °F?

## Your answer goes here:

Plot the distribution of the number of consecutive days over 100 °F (i.e., heat spell duration).

## Your answer goes here: