Terms and Definitions
CDD: “cooling degree days”; measured as the number of degrees above 65 degrees Farenheit, eg. If it’s 75 degrees Farenheit, CDD = 10. If it’s 60 degrees Farenheit, CDD = 0
Heat Index: a function of temperature and humidity
“Extreme Danger Days”: according to weather.gov, when the heat index is above 126
library(googleCloudStorageR)
library(ncdf4)
library(lubridate)##
## Attaching package: 'lubridate'## The following object is masked from 'package:base':
##
## datelibrary(pracma)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(sf)## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0library(leaflet)
library(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following objects are masked from 'package:lubridate':
##
## hour, isoweek, mday, minute, month, quarter, second, wday, week,
## yday, yearcreate_geometry = function(ID = 1, lat = 45, lon = -100, cellsize = 1){
# create data frame with a corner points, the fifth corner is a repeat of the first (needs to be a closed shape)
df <- data.frame(
longitude = c(lon - cellsize/2, lon - cellsize/2, lon + cellsize/2, lon + cellsize/2, lon - cellsize/2),
latitude = c(lat - cellsize/2, lat + cellsize/2, lat + cellsize/2, lat - cellsize/2, lat - cellsize/2 )
)
# convert to a geometry
polygon <- df %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
summarize(geometry = st_combine(geometry)) %>%
st_cast("POLYGON")
return(data.frame(ID = ID, lat = lat, lon = lon, geometry = polygon$geometry))
}Read in the pre-processed heat index data
# ~ 4GB file as of 2021-01-12; should read in with fread in a few seconds
heat_index_data = fread('../processed_data/heat_index_data.csv') Determine the number of warning days in each heat index category for each grid point According to the Economic Risks of Climate Change we can expect the following:
agged_history <- heat_index_data %>%
group_by(ID, lat, lon, year) %>%
summarise(hist_caution = sum(heat_index >= 80),
hist_extreme_caution = sum(heat_index >= 91),
hist_danger = sum(heat_index >= 103),
hist_extreme_danger = sum(heat_index >= 126),
Caution_rcp85_2050 = sum(heat_index_rcp85_2050 >= 80),
Extreme_Caution_rcp85_2050 = sum(heat_index_rcp85_2050 >= 91),
Danger_rcp85_2050 = sum(heat_index_rcp85_2050 >= 103),
Extreme_Danger_rcp85_2050 = sum(heat_index_rcp85_2050 >= 126),
cdd_total = sum(cdd),
cdd_avg = mean(cdd),
cdd_2040_rcp85_total = sum(cdd_2040), # note - these are 2040 rcp8.5
cdd_2040_rcp_85_avg = mean(cdd_2040),
cdd_hist_rcp85_delta = sum(cdd_2040) - sum(cdd),
cdd_avg_delta = mean(cdd_2040) - mean(cdd),
demand_delta = 0.00742*(sum(cdd_2040) - sum(cdd)),
Danger_change_rcp45_2040 = sum(heat_index_rcp45_2040 >= 103) - sum(heat_index >= 103),
Caution_change_rcp45_2040 = sum(heat_index_rcp45_2040 >= 80) - sum(heat_index >= 80),
Danger_change_rcp85_2050 = sum(heat_index_rcp85_2050 >= 103) - sum(heat_index >= 103),
Caution_change_rcp85_2050 = sum(heat_index_rcp85_2050 >= 80) - sum(heat_index >= 80),
Ext_Caution_change_rcp85_2050 = sum(heat_index_rcp85_2050 >= 91) - sum(heat_index >= 91)
) %>%
dplyr::ungroup()## `summarise()` regrouping output by 'ID', 'lat', 'lon' (override with `.groups` argument)record_summary <- agged_history %>%
dplyr::select(ID, lat, lon, hist_caution, hist_danger, hist_extreme_caution, Danger_change_rcp85_2050, Caution_change_rcp85_2050, Ext_Caution_change_rcp85_2050, cdd_total, cdd_hist_rcp85_delta) %>%
group_by(ID, lat, lon) %>%
summarise(Avg_Caution = mean(hist_caution),
Avg_Danger = mean(hist_danger),
Avg_Extreme_Caution = mean(hist_extreme_caution),
Max_Caution = max(hist_caution),
Avg_CDD = mean(cdd_total),
Max_CDD = max(cdd_total),
Avg_Danger_Change_rcp85_2050 = mean(Danger_change_rcp85_2050),
Avg_Ext_Caution_Change_rcp85_2050 = mean(Ext_Caution_change_rcp85_2050),
Avg_Caution_Change_rcp85_2050 = mean(Caution_change_rcp85_2050),
Avg_CDD_Change_rcp85_2040 = mean(cdd_hist_rcp85_delta)
) %>%
dplyr::ungroup()## `summarise()` regrouping output by 'ID', 'lat' (override with `.groups` argument)land_df = fread('../processed_data/landmask_df.csv')
cellsize = 1
coords_sf <- purrr::map(.x = 1:nrow(land_df),
.f = ~ create_geometry(ID = land_df$ID[.x], lat = land_df$lat[.x],
lon = land_df$lon[.x],cellsize = cellsize)) %>%
dplyr::bind_rows()
# add geometry
summary_sf <- coords_sf %>%
left_join(record_summary) %>%
sf::st_as_sf()## Joining, by = c("ID", "lat", "lon")# convert NAs to 0
summary_sf[is.na(summary_sf)] <- 0
summary_sf <- summary_sf %>% filter(lat < 49, lat > 25) # trim the grid to leave off parts of Canada and Mexicopal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_Caution_Change_rcp85_2050)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_Caution_Change_rcp85_2050,
fillColor = ~ pal(Avg_Caution_Change_rcp85_2050), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_Caution_Change_rcp85_2050,
pal = pal,
title = "Additional Caution Days<br> rcp85 2050")pal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_Caution)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_Caution,
fillColor = ~ pal(Avg_Caution), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_Caution,
pal = pal,
title = "Historical Average<br> Caution Days")summary_sf <- summary_sf %>% mutate(pct_change_Caution_rcp85_2050 = case_when(Avg_Caution_Change_rcp85_2050 == Avg_Caution ~ 0,
Avg_Caution == 0 ~ 1,
Avg_Caution > 0 ~ Avg_Caution_Change_rcp85_2050/ Avg_Caution))
pal <- colorNumeric(palette = "YlOrRd",
domain = c(0, 1))
summary_sf %>%
filter(pct_change_Caution_rcp85_2050 < 1) %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~pct_change_Caution_rcp85_2050,
fillColor = ~ pal(pct_change_Caution_rcp85_2050), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = c(0,1),
pal = pal,
title = "Percent Change Caution Days<br> RCP85 2050")I’ve trimmed out all the grid points where percent change exceeded 100%; due to some grid points having an average of say .02 caution days per year in the historical record and expecting to increase to an average of 3 caution days per year under RCP85 2050, the plot starts to get hard to interperet.
# TODO: see if I can recreate this plot with some way to trim out the "outliers" in percent changepal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_Ext_Caution_Change_rcp85_2050)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_Ext_Caution_Change_rcp85_2050,
fillColor = ~ pal(Avg_Ext_Caution_Change_rcp85_2050), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_Ext_Caution_Change_rcp85_2050,
pal = pal,
title = "Additional Extreme Caution Days<br> rcp85 2050")pal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_Extreme_Caution)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_Extreme_Caution,
fillColor = ~ pal(Avg_Extreme_Caution), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_Extreme_Caution,
pal = pal,
title = "Historical Average<br> Extreme Caution Days")summary_sf <- summary_sf %>% mutate(pct_change_Ext_Caution_rcp85_2050 = case_when(Avg_Ext_Caution_Change_rcp85_2050 == Avg_Extreme_Caution ~ 0,
Avg_Extreme_Caution == 0 ~ 1,
Avg_Extreme_Caution > 0 ~ Avg_Ext_Caution_Change_rcp85_2050/ Avg_Extreme_Caution))
pal <- colorNumeric(palette = "YlOrRd",
domain = c(3, 94))
summary_sf %>%
filter(pct_change_Ext_Caution_rcp85_2050 >= 3) %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~pct_change_Ext_Caution_rcp85_2050,
fillColor = ~ pal(pct_change_Ext_Caution_rcp85_2050), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = c(3,94),
pal = pal,
title = "Percent Change <br> Extreme Caution Days <br> RCP85 2050")pal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_CDD_Change_rcp85_2040)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_CDD_Change_rcp85_2040,
fillColor = ~ pal(Avg_CDD_Change_rcp85_2040), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_CDD_Change_rcp85_2040,
pal = pal,
title = "Additional CDD Days<br> rcp85 2040")pal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_CDD)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_CDD,
fillColor = ~ pal(Avg_CDD), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_CDD,
pal = pal,
title = "Average Historical<br> CDD Days")summary_sf <- summary_sf %>% mutate(pct_change_cdd_rcp85_2040 = case_when(Avg_CDD_Change_rcp85_2040 == Avg_CDD ~ 0,
Avg_CDD == 0 ~ 1,
Avg_CDD > 0 ~ Avg_CDD_Change_rcp85_2040/ Avg_CDD))
pal <- colorNumeric(palette = "YlOrRd",
domain = c(0, 1.3))
summary_sf %>%
filter(pct_change_cdd_rcp85_2040 < 1.3) %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~pct_change_cdd_rcp85_2040,
fillColor = ~ pal(pct_change_cdd_rcp85_2040), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = c(0,1.3),
pal = pal,
title = "Percent Change Historical CDD<br> RCP85 2040")There doesn’t seem to be a ton of change in Danger days, compared to caution.
pal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_Danger_Change_rcp85_2050)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_Danger_Change_rcp85_2050,
fillColor = ~ pal(Avg_Danger_Change_rcp85_2050), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_Danger_Change_rcp85_2050,
pal = pal,
title = "Additional Danger Days<br> rcp85 2050")Historically, it doesn’t look like heat_indexes of “Danger” have registered all that frequently
pal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$Avg_Danger)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~Avg_Danger,
fillColor = ~ pal(Avg_Danger), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$Avg_Danger,
pal = pal,
title = "Historical Average<br> Danger Days")Compile a daily average of CDD for each gridpoint
daily_summer_history <- heat_index_data %>%
filter(month %in% c(6,7,8,9)) %>%
group_by(ID, lat, lon) %>%
summarise(cdd_historical_avg = mean(cdd),
cdd_2040_rcp_85_avg = mean(cdd_2040),
cdd_avg_delta = mean(cdd_2040_rcp_85_avg) - mean(cdd_historical_avg)
) %>%
dplyr::ungroup()## `summarise()` regrouping output by 'ID', 'lat' (override with `.groups` argument)Add geometry
summary_sf <- coords_sf %>%
left_join(daily_summer_history) %>%
sf::st_as_sf()## Joining, by = c("ID", "lat", "lon")# convert NAs to 0
summary_sf[is.na(summary_sf)] <- 0
summary_sf <- summary_sf %>% filter(lat < 49, lat > 25) # trim the grid to leave off parts of Canada and Mexicopal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$cdd_historical_avg)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~cdd_historical_avg,
fillColor = ~ pal(cdd_historical_avg), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$cdd_historical_avg,
pal = pal,
title = "Average Daily <br> Summer CDD<br> rcp85 2050")pal <- colorNumeric(palette = "YlOrRd",
domain = summary_sf$cdd_avg_delta)
summary_sf %>%
leaflet() %>%
addTiles() %>%
addPolygons(layerId = ~cdd_avg_delta,
fillColor = ~ pal(cdd_avg_delta), color = "black",
weight = 0.2, fillOpacity = 0.8
) %>%
addLegend(position = "topleft",
values = summary_sf$cdd_avg_delta,
pal = pal,
title = "Additional Daily <br> Summer CDD<br> rcp85 2050")