library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readxl)
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tmap)
library(RColorBrewer)
source("functions.R")
# reading the four degree files for different global warming levels
# compared to pre-industrial 1851-1900
# temperature and precipitation in the data reflect changes in
# degrees F and percentage change for precipitation
# compared to RECENT 1990-2020 values
# temperature changes are in degrees F
# precipitation changes are percentage increases or decreases
d1 = read_csv("NCA5_Atlas_Global_Warming_Level_1.5_deg_C.csv") %>%
tolow() %>%
pivot_longer(cols=6:22,names_to="varname",values_to="datavalue")
## Rows: 3111 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): NAME, STATE_NAME, STATE_ABBR, FIPS
## dbl (18): OBJECTID, pr_above_nonzero_99th_GWL1, prmax1day_GWL1, prmax5yr_GWL...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d2 = read_csv("NCA5_Atlas_Global_Warming_Level_2_deg_C.csv") %>%
tolow() %>%
pivot_longer(cols=6:22,names_to="varname",values_to="datavalue")
## Rows: 3111 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): NAME, STATE_NAME, STATE_ABBR, FIPS
## dbl (18): OBJECTID, pr_above_nonzero_99th_GWL2, prmax1day_GWL2, prmax5yr_GWL...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d3 = read_csv("NCA5_Atlas_Global_Warming_Level_3_deg_C.csv") %>%
tolow() %>%
pivot_longer(cols=6:22,names_to="varname",values_to="datavalue")
## Rows: 3111 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): NAME, STATE_NAME, STATE_ABBR, FIPS
## dbl (18): OBJECTID, pr_above_nonzero_99th_GWL3, prmax1day_GWL3, prmax5yr_GWL...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d4 = read_csv("NCA5_Atlas_Global_Warming_Level_4_deg_C.csv") %>%
tolow() %>%
pivot_longer(cols=6:22,names_to="varname",values_to="datavalue")
## Rows: 3111 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): NAME, STATE_NAME, STATE_ABBR, FIPS
## dbl (18): OBJECTID, pr_above_nonzero_99th_GWL4, prmax1day_GWL4, prmax5yr_GWL...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
vardesc = read_excel("varnames1.xlsx")
# These Global Warming Levels (GWLs) correspond to global average temperature
# increases of 1.5, 2, 3, and 4 degrees Celsius
# above pre-industrial levels measured from 1851 to 1900.
# Historical analyses were based on NOAA's nClimGrid and GHCN-Daily datasets.
# On the Fahrenheit scale, these warming levels are
# 2.7 ?F (1.5 ?C),
# 3.6 ?F (2 ?C),
# 5.4 ?F (3 ?C), and
# 7.2 ?F (4 ?C).
# As of the 2020s, global average temperature has already increased
# around 2 ?F (1.1 ?C) above pre-industrial levels.
# The fgwl values below are increases from the HISTORIC 1850-1900
# global average temperature.
# The adjfgwl values below are increases from the RECENT
# global average temperature, which is 2? F higher than
# the HISTORIC fgwl values.
# creating single long file
# pulling substrings from the variable names
d5 = d1 %>%
bind_rows(d2) %>%
bind_rows(d3) %>%
bind_rows(d4) %>%
filter(!varname %in% c("shape_length", "shape_area")) %>%
mutate(geocode = paste0("g",fips),
coname = paste0(name,", ",state_abbr)) %>%
mutate(varname0 = str_sub(varname, end=-6)) %>%
mutate(gwl = str_sub(varname,start=-4)) %>%
mutate(fgwl = case_when(
gwl == "gwl1" ~ "fgwl27",
gwl == "gwl2" ~ "fgwl36",
gwl == "gwl3" ~ "fgwl54",
gwl == "gwl4" ~ "fgwl72")) %>%
mutate(adjfgwl = case_when(
gwl == "gwl1" ~ "adjfgwl07",
gwl == "gwl2" ~ "adjfgwl16",
gwl == "gwl3" ~ "adjfgwl34",
gwl == "gwl4" ~ "adjfgwl52")) %>%
mutate(gwlcombo = paste0(gwl,"-",fgwl)) %>%
select(geocode,coname,gwl,fgwl,adjfgwl,
varname,varname0,datavalue) %>%
left_join(vardesc)
## Joining with `by = join_by(varname0)`
# Here's a summary of the gwl equivalents
#
# Similar
# Historic Historic Recent Scenario
# Celsius Fahrenheit Fahrenheit for Year 2100
# Baseline Baseline Baseline Temperature
#
# gwl1 fgwl27 adjfgwl07 SSP1-RCP1.9
# gwl2 fgwl36 adjfgwl16 SSP1-RCP2.6
# gwl3 fgwl54 adjfgwl34 SSP2-RCP4.5
# gwl4 fgwl72 adjfgwl52 SSP3-RCP7.0
# gwl5 fgwl90 adjfgwl70 SSP5-RCP8.5
s1 = st_read("cb_2018_us_state_20m.shp") %>%
tolow()
## Reading layer `cb_2018_us_state_20m' from data source
## `/Users/apple/Desktop/UNI/FALL 2023/1. CCA/Class/class2523/cb_2018_us_state_20m.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 52 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
## Geodetic CRS: NAD83
c1 = st_read("cb_2018_us_county_20m.shp") %>%
tolow() %>%
mutate(geocode = paste0("g", geoid)) %>%
filter(!mid(geocode,2,2) %in% c("02","15","72")) %>%
select(geocode)
## Reading layer `cb_2018_us_county_20m' from data source
## `/Users/apple/Desktop/UNI/FALL 2023/1. CCA/Class/class2523/cb_2018_us_county_20m.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3220 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.1743 ymin: 17.91377 xmax: 179.7739 ymax: 71.35256
## Geodetic CRS: NAD83
# selecting tavg temperature average variable
# temperature is degree increase compared to 1991-2020
# calculating the temperature difference
# between gwl 3 degrees C and gwl 4 degrees C
# which is about the difference between
# RCP4.5 and RCP7.0 in 2100
t5 = d5 %>%
# filter(mid(geocode,1,3) == "g13") %>%
filter(varname0 == "tavg") %>%
select(geocode,varname0,adjfgwl,datavalue) %>%
pivot_wider(names_from="adjfgwl", values_from="datavalue") %>%
filter(!mid(geocode,2,2) %in% c("02","15","72")) %>%
mutate(cgwl34diff = adjfgwl34-adjfgwl16)
t6 = c1 %>%
right_join(t5)
## Joining with `by = join_by(geocode)`
# fixing a shapefile data incompatibility
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
tmap_mode("view")
## tmap mode set to interactive viewing
# mapping gwl4 = fgwl54 = adjfgwl34
# which is close to RCP 7.0, BAU
tm_basemap("OpenStreetMap") +
tm_shape(t6) +
tm_fill("adjfgwl34",alpha=.5, palette="YlOrRd") +
tm_borders(lwd=2) +
tm_layout(title="Change in Annual Temperature from Present with BAU")
# mapping the lower temperature difference
# between gwl 3 degrees C and gwl 4 degrees C
# which is about the difference between
# RCP4.5 and RCP7.0 in 2100
tm_basemap("OpenStreetMap") +
tm_shape(t6) +
tm_fill("cgwl34diff",alpha=.5, palette="PiYG") +
tm_shape(s1) +
tm_borders(lwd=2) +
tm_layout(title="Decrease in Annual Temperature from BAU with Climate Action")
# selecting pr_annual annual precipitation variable
# precip is increase compared to 1991-2020
p5 = d5 %>%
# filter(mid(geocode,1,3) == "g13") %>%
filter(varname0 == "pr_annual") %>%
select(geocode,varname0,adjfgwl,datavalue) %>%
pivot_wider(names_from="adjfgwl", values_from="datavalue") %>%
filter(!mid(geocode,2,2) %in% c("02","15","72"))
p6 = c1 %>%
right_join(p5)
## Joining with `by = join_by(geocode)`
sf_use_s2(FALSE)
# mapping precipitation increase
tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap") +
tm_shape(p6) +
tm_fill("adjfgwl34",alpha=.5, palette="Spectral",n=10) +
tm_shape(s1) +
tm_borders(lwd=2) +
tm_layout(title="Change in Annual Precipitation from Present with BAU")
## Variable(s) "adjfgwl34" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
# selecting annual precipitation
# calculating the precipitation difference
# between gwl 3 degrees C and gwl 4 degrees C
# which is about the difference between
# RCP4.5 and RCP7.0 in 2100
p5diff = d5 %>%
# filter(mid(geocode,1,3) == "g13") %>%
filter(varname0 == "pr_annual") %>%
select(geocode,varname0,adjfgwl,datavalue) %>%
pivot_wider(names_from="adjfgwl", values_from="datavalue") %>%
filter(!mid(geocode,2,2) %in% c("02","15","72")) %>%
mutate(cgwl34diff = adjfgwl52-adjfgwl34)
p6diff = c1 %>%
right_join(p5diff)
## Joining with `by = join_by(geocode)`
# mapping change in precipitation from BAU to climate action
sf_use_s2(FALSE)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap") +
tm_shape(p6diff) +
tm_fill("cgwl34diff",alpha=.5, palette="Spectral",n=10) +
tm_shape(s1) +
tm_borders(lwd=2) +
tm_layout(title="Change in Annual Precipitation from BAU with Climate Action")
## Variable(s) "cgwl34diff" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
# selecting maximum precipitation day in 5 years
# calculating the lower 5-year day max precipitation difference
# between gwl 3 degrees C and gwl 4 degrees C
# which is about the difference between
# RCP4.5 and RCP7.0 in 2100
p15 = d5 %>%
# filter(mid(geocode,1,3) == "g13") %>%
filter(varname0 == "prmax5yr") %>%
select(geocode,varname0,adjfgwl,datavalue) %>%
pivot_wider(names_from="adjfgwl", values_from="datavalue") %>%
filter(!mid(geocode,2,2) %in% c("02","15","72")) %>%
mutate(cgwl34diff = adjfgwl52-adjfgwl34)
p16 = c1 %>%
right_join(p15)
## Joining with `by = join_by(geocode)`
# mapping the 5-year day maximum precipitation with gwl 4 degrees
tm_basemap("OpenStreetMap") +
tm_shape(p6) +
tm_fill("adjfgwl52",alpha=.5, palette="Spectral",n=10) +
tm_shape(s1) +
tm_borders(lwd=2) +
tm_layout(title="Change in Highest Precipitation Day in Five Years from Present with BAU")
## Variable(s) "adjfgwl52" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
# mapping the 5-year day maximum precipitation difference
# between gwl 3 degrees C and gwl 4 degrees C
# which is about the difference between
# RCP4.5 and RCP7.0 in 2100
tm_basemap("OpenStreetMap") +
tm_shape(p16) +
tm_fill("cgwl34diff",alpha=.5, palette="Spectral",n=10) +
tm_shape(s1) +
tm_borders(lwd=2) +
tm_layout(title="Change in Highest Precipitation Day in Five Years from BAU with Climate Action")
## Variable(s) "cgwl34diff" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.