solutions 

Getting Started with caladaptR

Overview

This Notebook will demonstrate how you can use caladaptR to:

• create an API Request using a preset area-of-interest (county)
• use pivot_wider to go from a long to wide format
  
• add a trend line to a plot
  
• explore the Cal-Adapt data catalog
  
• specify Livneh datasets
  
• wrangle results for different summaries and visualizations

Setup

The first thing we do is to load caladaptR and the other package we’re going to need. (If you haven’t installed these yet, see this setup script).

library(caladaptr)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
caladaptr (version 0.6.4)
URL: https://ucanr-igis.github.io/caladaptr
Bug reports: https://github.com/ucanr-igis/caladaptr/issues

library(units)

udunits database from C:/Users/Andy/Documents/R/win-library/4.1/units/share/udunits/udunits2.xml

library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(sf)

Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1

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

Part I. Retrieve County Level Data

In this section, we’ll fetch and retrieve average daily minimum temperature by year for a single county, and create a time-series plot showing the difference between RCP85 and RCP45:

Preset Areas-of-Interest

For this exercise, we need to use a Preset Area-of-Interest.

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"            

Pro Tip:

• all of the Preset AOIs can be imported into R as sf objects with ca_aoipreset_geom().

1. Find the FIPS code for Kings County

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 View the attribute table of the counties preset layer:

# ca_aoipreset_geom("counties") %>% View()

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.

cnty_cap <- 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")

cnty_cap

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
 

cnty_cap %>% ca_preflight()

General issues
 - none found
Issues for querying values
 - none found
Issues for downloading rasters
 - none found

Pro Tip:

• An API request can have more than one location. If using AOI Presets, pass multiple values to idval, or omit idval completely and all features in the layer will be queried.

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

plot(cnty_cap, locagrid = TRUE)

3. Fetch Data

Fetch data with ca_getvals_tbl():

cnty_tbl <- cnty_cap %>% ca_getvals_tbl(quiet = TRUE)

cnty_tbl %>% head()

4. Wrangle Results for Plotting

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:

cnty_diff_rcp85_45_tbl <- cnty_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(cnty_diff_rcp85_45_tbl)

5. Plot

Now we’re ready to make the plot. Since we’re mainly interested in the trend, we’ll add a smoothing line using geom_smooth():

ggplot(data = cnty_diff_rcp85_45_tbl, aes(x = as.Date(dt), y = as.numeric(rcp85 - rcp45))) +
  geom_line(aes(color=gcm)) +
  geom_smooth(method=lm, formula = y ~ x) +
  labs(title = "Difference between RCP8.5 and RCP4.5 in the Average Daily \nMinimum Temperature for Kings County", x = "year", y = "temp (F)")

This plot shows that as time goes on, the difference between RCP4.5 and RCP8.5 gets bigger and bigger.

Part II. Finding Data

Half the battle of working with climate data is finding the name of the dataset you’re interested in.

Browsing the Cal-Adapt Data Catalog

caladaptR comes with a copy of the Cal-Adapt raster series data catalog, which you can access with ca_catalog_rs():

# ca_catalog_rs() %>% View()

Challenge 1

How many raster datasets are there in the Cal-Adapt catalog? [Answer]

ca_catalog_rs() %>% nrow()

 - using raster series catalog from cache

[1] 959

Pro Tip:

• you can download a fresh copy of the Cal-Adapt raster series catalog with ca_catalog_fetch()

Searching for a Dataset

One way you can search for datasets is to use the filter boxes above each column in the RStudio View pane. For example search for layers whose name contains the word ‘snow’.

caladaptR also has a search function ca_catalog_search(). You can this function to find datasets using a key word or phrase. You can also use this function to see the details for a specific slug, example:

ca_catalog_search("swe_day_ens32avg_rcp45")

swe_day_ens32avg_rcp45
  name: LOCA VIC daily snow water equivalent for RCP 4.5, derived from 32 LOCA models ensemble average
  url: https://api.cal-adapt.org/api/series/swe_day_ens32avg_rcp45/
  tres: daily
  begin: 2006-01-01T00:00:00Z
  end: 2098-12-31T00:00:00Z
  units: mm
  num_rast: 1
  id: 638
  xmin: -124.5625
  xmax: -113.375
  ymin: 31.5625
  ymax: 43.75

Challenge 2

How many datasets are from gridMET? [Answer]

ca_catalog_search("gridmet", quiet = TRUE) %>% nrow()

[1] 6

## Ans. 6

Specifying Datasets

When you construct an API Request object, you can mix-and-match functions to specify the dataset you want:

LOCA downscaled modeled climate variables (including all Scripps) and their derivatives (i.e., VIC) can be specified using the constructor functions: ca_gcm() + ca_scenario() + ca_cvar() + ca_period().

Livneh data (observed historical variables based on spatially interpolated measurements) can be specified with ca_livneh() + ca_cvar() + ca_period().

everything else can be specified by slug, using ca_slug().

Challenge 3

What is the raster dataset with the slug tmmn_day_gridmet? For which years is it available? What are the units? [Hint] | [Answer]

Ans. gridMET daily minimum temperature historical.

1979-2020

Kelvin

ca_catalog_search("tmmn_day_gridmet")

tmmn_day_gridmet
  name: gridMET daily minimum temperature historical
  url: https://api.cal-adapt.org/api/series/tmmn_day_gridmet/
  tres: daily
  begin: 1979-01-01T00:00:00Z
  end: 2020-12-31T00:00:00Z
  units: K
  num_rast: 1
  id: 943
  xmin: -124.579167
  xmax: -113.370833
  ymin: 31.545833
  ymax: 43.754167

Explore the Climate Future of Lindcove Research and Extension Center

In this example, we’ll explore daily projected climate data for a point location. We select the UC Lindcove Research and Extension Center (LREC), a field station in the San Joaquin Valley which has been a leading center for citrus research for decades. The trees are getting old and need to be replaced soon so its worth asking - is citrus production still going to viable in this part of California at the end of the century?

Create the API Request

Let’s start by getting 20 years worth of the daily maximum temperature for the 4 priority GCMs and 2 RCPs.

lrec_tasmax_prj_cap <- ca_loc_pt(coords = c(-119.060, 36.359), id = 1) %>% 
  ca_period("day") %>% 
  ca_gcm(gcms[1:4]) %>% 
  ca_scenario(c("rcp45", "rcp85")) %>% 
  ca_cvar("tasmax") %>% 
  ca_years(start = 2080, end = 2099)
  
lrec_tasmax_prj_cap %>% ca_preflight()

General issues
 - none found
Issues for querying values
 - none found
Issues for downloading rasters
 - none found

plot(lrec_tasmax_prj_cap)

Fetch Data

Now we’re ready to fetch data:

lrec_tasmax_prj_tbl <- lrec_tasmax_prj_cap %>% ca_getvals_tbl(quiet = TRUE)

## backup: lrec_tasmax_prj_tbl <- readRDS("data/lrec_tasmax_prj_tbl.rds")

dim(lrec_tasmax_prj_tbl)

[1] 58440     8

head(lrec_tasmax_prj_tbl)

Create a Box Plot of the Maximum Daily Temperature Values by Month

To make a box plot, we need to first add columns for Fahrenheit, month, and year:

library(lubridate)
lrec_tasmax_prj_tmpf_tbl <- lrec_tasmax_prj_tbl %>% 
  mutate(temp_f = set_units(val, degF), month = month(dt), year = year(dt))
head(lrec_tasmax_prj_tmpf_tbl)

For each month, let make a box plot of the temperature values for each emission scenario, treating all GCMs as equally likely:

ggplot(lrec_tasmax_prj_tmpf_tbl, aes(x = as.factor(month), y = as.numeric(temp_f))) + 
  geom_boxplot() +
  facet_grid(scenario ~ .) +
  labs(title = "Maximum Daily Temperature by Month", x = "month", y = "temp (F)",
       subtitle = "Lindcove REC, 4 GCMs combined, 2080-2099")

Count Extreme Heat Days

An extreme heat day is generally identified when the maximum temperature exceeds a threshold. The threshold can be chosen based on the historical range, or a biophysical process. For this example, we’ll select 105 °F.

Let’s count the total number of days the temperature exceeded 105 °F for each RCP. We start by adding a logical column whether the temperature exceeded our threshold.

lrec_hotday_tbl <- lrec_tasmax_prj_tmpf_tbl %>% 
  mutate(really_hot = (temp_f >= set_units(105, degF))) %>% 
  select(-spag, -val, -month)

head(lrec_hotday_tbl)

We can count the number of extreme heat days with a simple expression:

num_hot_days <- lrec_hotday_tbl %>%
  group_by(scenario, really_hot) %>% 
  count()

num_hot_days

We can improve the readability of this table by making each scenario a separate column. This is an example of pivoting, which you can handle using tidyr::pivot_wider.

num_hot_days %>% pivot_wider(names_from = scenario, values_from = n)

Challenge 4

Count the number of extreme heat days using a threshold of 110 Fahrenheit. Compute the number of extreme heat days per month by scenario. [Answer]

lrec_hothotday_tbl <- lrec_tasmax_prj_tmpf_tbl %>% 
  mutate(hot_hot = (temp_f >= set_units(110, degF))) %>% 
  select(-spag, -val)

head(lrec_hothotday_tbl)

lrec_hothotday_tbl %>% 
  filter(hot_hot) %>% 
  group_by(month, scenario) %>% 
  summarise(count = n()) %>% 
  pivot_wider(names_from = month, values_from = count)

`summarise()` has grouped output by 'month'. You can override using the `.groups` argument.

Pro Tip

• To see how to count the number consecutive heat days (i.e., heat spells), see this notebook.

Visualize the Distribution of Historical Observed Daily Precipitation by Decade

Our goal here is to look at the distribution of daily precipitation by making a histogram for each decade. We’ll use observed rainfall data from Livneh.

lrec_pr_livn_cap <- ca_loc_pt(coords = c(-119.060, 36.359), id = 1) %>% 
  ca_livneh(TRUE) %>% 
  ca_period("day") %>% 
  ca_cvar("pr") %>% 
  ca_years(start = 1950, end = 2009)

lrec_pr_livn_cap %>% ca_preflight()

General issues
 - none found
Issues for querying values
 - none found
Issues for downloading rasters
 - none found

lrec_pr_livn_tbl <- lrec_pr_livn_cap %>% 
  ca_getvals_tbl() %>% 
  rename(pr_mmday = val)  

## backup: lrec_pr_prj_tbl <- readRDS("./data/lrec_pr_prj_tbl.rds")

dim(lrec_pr_livn_tbl)

[1] 21915     7

head(lrec_pr_livn_tbl)

Finally, we’ll plot histograms of the precipitation values, logged because of the highly skewed distribution.

library(lubridate)

lrec_prdec_livn_tbl <- lrec_pr_livn_tbl %>% 
  mutate(pr_mmday_num = as.numeric(pr_mmday),
         decade = floor(year(dt) / 10) * 10) %>% 
  select(decade, pr_mmday_num)

ggplot(lrec_prdec_livn_tbl, aes(x=log(pr_mmday_num))) + 
  geom_histogram() +
  facet_wrap( ~ decade)

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 17073 rows containing non-finite values (stat_bin).

Conclusion

Once you get climate data in R as data frames, there’s a lot you can do with it!

---