Loading R/plots.R 0 → 100644 +547 −0 Original line number Diff line number Diff line # CF trend plot plot_CF_trend_map <- function(CFs, CF_trends, EHA){ # plot background (greyscale shaded relief) read_stars("./data/Natural Earth/US_MSR_10M/US_MSR.tif", RasterIO = list(nBufXSize = 600, nBufYSize = 600)) -> NA_shaded_relief # major rivers intersecting the US st_read("./data/Natural Earth/10m_physical/ne_10m_rivers_lake_centerlines.shp") |> st_transform(st_crs(NA_shaded_relief)) -> major_rivers major_rivers[unlist(st_intersects(lower48_shp, major_rivers)),] -> major_rivers_US # US state boundaries st_read("./data/Natural Earth/10m_cultural/ne_10m_admin_1_states_provinces_lakes.shp") |> filter(admin == "United States of America", !name_en %in% c("Alaska", "Hawaii")) |> st_transform(st_crs(NA_shaded_relief)) -> lower48_shp # shapes for four major hydropower regions lower48_shp |> subset(name_en == "California") -> shp_CA lower48_shp |> subset(name_en %in% c("Washington", "Oregon", "Idaho")) |> st_union() -> shp_PNW lower48_shp |> subset(name_en %in% c( "Tennessee", "Alabama", "Georgia", "South Carolina")) |> st_union() -> shp_SE lower48_shp |> subset(name_en %in% c( "New York", "Vermont", "Maine", "New Hampshire")) |> st_union() -> shp_NE # get CF trends and split into categories CF_trends |> filter(analysis == "923 period") |> mutate(trend_CFpp_per_decade = sens_slope * 10 * 100) |> select(COMPLXID, trend_CFpp_per_decade, p_value) -> trend_data_all_plants trend_data_all_plants |> filter(trend_CFpp_per_decade < 0 & p_value < 0.05) |> pull(COMPLXID) -> dams_sig_down trend_data_all_plants |> filter(trend_CFpp_per_decade > 0 & p_value < 0.05) |> pull(COMPLXID) -> dams_sig_up trend_data_all_plants |> filter(p_value > 0.05) |> pull(COMPLXID) -> dams_nonsig trend_data_all_plants |> filter(is.na(p_value)) |> pull(COMPLXID) -> dams_nodata st_read(EHA) |> mutate(COMPLXID = str_replace(EHA_PtID, "\\_..*", "")) |> filter(COMPLXID %in% CFs$COMPLXID) |> st_transform(st_crs(NA_shaded_relief)) |> mutate(size_cat = case_when( CH_MW > 1000 ~ "> 1 ", CH_MW > 500 & CH_MW <= 1000 ~ "0.5 - 1 ", CH_MW > 100 & CH_MW <= 500 ~ "0.1 - 0.5 ", CH_MW >= 5 & CH_MW <= 100 ~ "0.05 - 0.1 " )) |> mutate(size_cat = factor(size_cat, levels = c( "0.05 - 0.1 ", "0.1 - 0.5 ", "0.5 - 1 ", "> 1 " ))) |> left_join(trend_data_all_plants, by = join_by(COMPLXID)) |> filter(!State %in% c("AK", "HI")) -> dam_points filter(dam_points, COMPLXID %in% dams_nonsig) -> dam_points_nosig filter(dam_points, COMPLXID %in% dams_sig_down) -> dam_points_sig_down filter(dam_points, COMPLXID %in% dams_sig_up) -> dam_points_sig_up filter(dam_points, COMPLXID %in% dams_nodata) -> dam_points_nodata # main map plot ggplot() + geom_stars(data = NA_shaded_relief |> st_crop(lower48_shp), show.legend = FALSE) + scale_fill_gradient(low = "grey70", high = "grey95", na.value = NA) + coord_equal() + geom_sf(data = lower48_shp, fill = NA, col = "grey20") + geom_sf(data = shp_CA, col = "black",fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = shp_PNW, col = "black", fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = shp_SE, col = "black", fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = shp_NE, col = "black", fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = major_rivers_US, col = "#0070B9", lwd = 0.05) + theme_void() + new_scale_fill() + geom_sf(data = dam_points_nodata, aes(size = size_cat), shape = 21, col = "black", fill = NA) + geom_sf(data = dam_points_nosig, aes(size = size_cat), shape = 21, col = "black", alpha = 0.8, fill = "white") + geom_sf(data = dam_points_sig_down, aes(size = size_cat, fill = trend_CFpp_per_decade), shape = 21, col = "black", alpha = 0.8) + geom_sf(data = dam_points_sig_up, aes(size = size_cat), shape = 21, col = "black", alpha = 0.8, fill = "dodgerblue") + scale_fill_viridis(option = "A") + scale_size_manual(values = c( "> 1 " = 7, "0.5 - 1 " = 5, "0.1 - 0.5 " = 2.5, "0.05 - 0.1 " = 1 )) + theme(legend.position = c(0.5,0.99), legend.direction = "horizontal", legend.box = "vertical", legend.background = element_rect(fill = "white", color = "white")) + labs( fill = "Trend in Capacity Factor (PPPD) \n (neg. trend w/ p<0.05 only)", size = "Nameplate in 2021 (GW) ", title = "a") + guides(size = guide_legend(override.aes=list(fill=NA), order = 1)) + annotate(geom = "point", size = 5, x = -90e5, y = 665e4, fill = "dodgerblue", col = "black", alpha = 0.8, pch = 21) + annotate("text", label = "pos. trend (p<0.05)", x = -89e5, y = 665e4, hjust = 0, size = 3) + annotate(geom = "point", size = 5, x = -90e5, y = 635e4, fill = "white", alpha = 0.8, pch = 21, col = "black") + annotate("text", label = "no significant trend", x = -89e5, y = 635e4, hjust = 0, size = 3) + annotate(geom = "point", size = 5, x = -90e5, y = 605e4, fill = "grey80", pch = 21, col = "black") + annotate("text", label = "insufficient data", x = -89e5, y = 605e4, hjust = 0, size = 3) + annotate("segment", x = -141e5, xend = -141e5, y = 560e4, yend = 240e4) + annotate("segment", x = -141e5, xend = -139e5, y = 560e4, yend = 560e4) + annotate("segment", x = -130e5, xend = -122e5, y = 375e4, yend = 240e4) + annotate("segment", x = -77e5, xend = -77e5, y = 535e4, yend = 240e4) + annotate("segment", x = -97e5, xend = -97e5, y = 345e4, yend = 240e4) + NULL -> map ### CA #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State == "CA")$COMPLXID) -> dam_gen_CA dam_gen_CA |> filter(!is.na(gen_MWh), !is.na(cap_MWh)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_CA senSlope(CF ~ year, data = CF_CA, intercept = "A1m") -> trend_CA CF_CA |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_CA$coefficients[2], intercept = trend_CA$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "c", subtitle = "[CF] California") -> p1_CA #### ### PNW #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State %in% c("WA", "OR", "ID"))$COMPLXID) -> dam_gen_PNW dam_gen_PNW |> filter(!is.na(gen_MWh), !is.na(cap_MWh)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_PNW senSlope(CF ~ year, data = CF_PNW, intercept = "A1m") -> trend_PNW CF_PNW |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_PNW$coefficients[2], intercept = trend_PNW$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "b", subtitle = "[CF] Pacific NW") -> p1_PNW #### ### SE #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State %in% c("TN", "AL", "GA", "SC"))$COMPLXID) -> dam_gen_SE dam_gen_SE |> filter(!is.na(gen_MWh), !is.na(cap_MWh)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_SE senSlope(CF ~ year, data = CF_SE, intercept = "A1m") -> trend_SE CF_SE |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_SE$coefficients[2], intercept = trend_SE$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "d", subtitle = "[CF] Southeast") -> p1_SE #### ### NE #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State %in% c("NY", "VT", "NH", "ME"))$COMPLXID) -> dam_gen_NE dam_gen_NE |> filter(!is.na(gen_MWh), !is.na(cap_MWh), !is.na(CF)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_NE senSlope(CF ~ year, data = CF_NE, intercept = "A1m") -> trend_NE CF_NE |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_NE$coefficients[2], intercept = trend_NE$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "e", subtitle = "[CF] Northeast") -> p1_NE #### # combine plots and save as pdf (p1_PNW + p1_CA + p1_SE + p1_NE + plot_layout(ncol = 4)) -> low_plots map + plot_spacer() + low_plots + plot_layout(ncol = 1, heights = c(0.75, -0.1 ,0.25)) -> final_plot ggsave( "./output/plots/CF_trend_map.pdf", final_plot, height = 8, width = 8 ) ggsave( "./output/plots/CF_trend_map.png", final_plot, height = 8, width = 8 ) return(final_plot) } plot_main_driver_of_trend <- function(CFs, CF_trends, EHA){ flow_CF_models_scored |> filter(flow_type == "c") |> filter(rsq > 0.2, abs(bias) < 0.1) |> arrange(-rsq) |> pull(COMPLXID) -> viable_flow_models CF_trends |> filter(analysis == "flow period", !is.na(p_value)) |> filter(p_value < 0.05) -> dams_with_significant_change flow_CF_models |> filter(COMPLXID %in% viable_flow_models) |> select(COMPLXID, year, CF, flow_type, CF_pred) |> pivot_wider(names_from = "flow_type", values_from = "CF_pred") |> split(~COMPLXID) |> map_dfr(function(x){ smwrStats::senSlope(CF ~ year, x) -> sen_CF smwrStats::senSlope(c ~ year, x) -> sen_flowC smwrStats::senSlope(n ~ year, x) -> sen_flowN bind_rows( tibble( slope = sen_CF$coefficients[2], data = "CF" ), tibble( slope = sen_flowC$coefficients[2], data = "Assimilated flow" ), tibble( slope = sen_flowN$coefficients[2], data = "Naturalized flow" ) ) |> mutate(COMPLXID = first(x$COMPLXID)) }) -> all_slopes all_slopes |> pivot_wider(names_from = "data", values_from = slope) |> mutate(explained_by_water = `Assimilated flow` / CF, explained_by_climate = `Naturalized flow` / CF) |> select(COMPLXID, CF_trend = CF, explained_by_water, explained_by_climate) |> mutate(water_not_climate = if_else( explained_by_water > 0.5 & explained_by_climate < 0.3, TRUE, FALSE)) |> mutate( flow_attrib = case_when( explained_by_water >= 0.8 ~ "Strongly driven by water availability", explained_by_water < 0.8 & explained_by_water >= 0.6 ~ "Mostly driven by water availability", explained_by_water < 0.6 & explained_by_water >= 0.4 ~ "~ Half of trend driven by water availability", explained_by_water < 0.4 & explained_by_water >= 0.2 ~ "Mostly driven by other factors", explained_by_water < 0.2 ~ "Strongly driven by other factors" )) |> select(COMPLXID, water_not_climate, flow_attrib) -> flow_cats dams_with_significant_change |> left_join(flow_cats, by = join_by(COMPLXID)) |> filter(sens_slope > 0) -> positive_cases dams_with_significant_change |> left_join(flow_cats, by = join_by(COMPLXID)) |> filter(sens_slope < 0) -> negative_cases st_read(EHA) |> mutate(COMPLXID = str_replace(EHA_PtID, "\\_..*", "")) |> filter(COMPLXID %in% dams_with_significant_change$COMPLXID) |> st_transform(st_crs(NA_shaded_relief)) |> mutate(size_cat = case_when( CH_MW > 1000 ~ "> 1 ", CH_MW > 500 & CH_MW <= 1000 ~ "0.5 - 1 ", CH_MW > 100 & CH_MW <= 500 ~ "0.1 - 0.5 ", CH_MW >= 5 & CH_MW <= 100 ~ "0.05 - 0.1 " )) |> mutate(size_cat = factor(size_cat, levels = c( "0.05 - 0.1 ", "0.1 - 0.5 ", "0.5 - 1 ", "> 1 " ))) |> filter(!State %in% c("AK", "HI")) -> dam_points dam_points |> filter(COMPLXID %in% positive_cases$COMPLXID) |> left_join(positive_cases) -> dam_points_positive dam_points |> filter(COMPLXID %in% negative_cases$COMPLXID) |> left_join(negative_cases) -> dam_points_negative ## map backgrounds ## # plot background (greyscale shaded relief) read_stars("./data/Natural Earth/US_MSR_10M/US_MSR.tif", RasterIO = list(nBufXSize = 600, nBufYSize = 600)) -> NA_shaded_relief # major rivers intersecting the US st_read("./data/Natural Earth/10m_physical/ne_10m_rivers_lake_centerlines.shp") |> st_transform(st_crs(NA_shaded_relief)) -> major_rivers major_rivers[unlist(st_intersects(lower48_shp, major_rivers)),] -> major_rivers_US # US state boundaries st_read("./data/Natural Earth/10m_cultural/ne_10m_admin_1_states_provinces_lakes.shp") |> filter(admin == "United States of America", !name_en %in% c("Alaska", "Hawaii")) |> st_transform(st_crs(NA_shaded_relief)) -> lower48_shp ggplot() + geom_stars(data = NA_shaded_relief |> st_crop(lower48_shp), show.legend = FALSE) + scale_fill_gradient(low = "grey70", high = "grey95", na.value = NA) + coord_equal() + geom_sf(data = lower48_shp, fill = NA, col = "grey20") + geom_sf(data = major_rivers_US, col = "#0070B9", lwd = 0.05) + theme_void() + new_scale_fill() + geom_sf(data = dam_points_negative, aes(fill = flow_attrib, size = size_cat), shape = 25) + geom_sf(data = dam_points_positive, aes(fill = flow_attrib, size = size_cat), shape = 24) + scale_size_manual(values = c( "> 1 " = 7, "0.5 - 1 " = 5, "0.1 - 0.5 " = 2.5, "0.05 - 0.1 " = 1 )) + scale_fill_manual(values = c( "Strongly driven by water availability" = "#440154FF", "Mostly driven by water availability" = "#3B528BFF", "~ Half of trend driven by water availability" = "#21908CFF", "Mostly driven by other factors" = "#5DC863FF", "Strongly driven by other factors" = "#FDE725FF" )) + NULL } _targets.R +24 −1 Original line number Diff line number Diff line Loading @@ -14,10 +14,16 @@ tar_option_set( "furrr", # ^^ apply in parallel "future", # parallel processing "sf", # geospatial tools "stars", # raster tools "foreign", # read old dbf data "zoo", # interpolation functions "smwrStats", # Sen's slope computation "stringr" # string manipulation functions "stringr", # string manipulation functions "terrainr", # map plotting "ggplot2", # plotting "ggnewscale",# plot support "patchwork", # plot combining "viridis" # color palettes ) ) Loading Loading @@ -104,5 +110,22 @@ list( CFs = dam_annual_gen_cap_CF_1970_2021, flows = DayFlow), format = "parquet" ), tar_target( flow_CF_models, simulate_annual_CF(data = model_data_ready), format = "parquet" ), tar_target( flow_CF_models_scored, evaluate_models(model_results = flow_CF_models), format = "parquet" ), tar_target( CF_trend_map, plot_CF_trend_map(CFs = dam_annual_gen_cap_CF_1970_2021, CF_trends = CF_trends, EHA = EHA), format = "rds" ) ) Loading
R/plots.R 0 → 100644 +547 −0 Original line number Diff line number Diff line # CF trend plot plot_CF_trend_map <- function(CFs, CF_trends, EHA){ # plot background (greyscale shaded relief) read_stars("./data/Natural Earth/US_MSR_10M/US_MSR.tif", RasterIO = list(nBufXSize = 600, nBufYSize = 600)) -> NA_shaded_relief # major rivers intersecting the US st_read("./data/Natural Earth/10m_physical/ne_10m_rivers_lake_centerlines.shp") |> st_transform(st_crs(NA_shaded_relief)) -> major_rivers major_rivers[unlist(st_intersects(lower48_shp, major_rivers)),] -> major_rivers_US # US state boundaries st_read("./data/Natural Earth/10m_cultural/ne_10m_admin_1_states_provinces_lakes.shp") |> filter(admin == "United States of America", !name_en %in% c("Alaska", "Hawaii")) |> st_transform(st_crs(NA_shaded_relief)) -> lower48_shp # shapes for four major hydropower regions lower48_shp |> subset(name_en == "California") -> shp_CA lower48_shp |> subset(name_en %in% c("Washington", "Oregon", "Idaho")) |> st_union() -> shp_PNW lower48_shp |> subset(name_en %in% c( "Tennessee", "Alabama", "Georgia", "South Carolina")) |> st_union() -> shp_SE lower48_shp |> subset(name_en %in% c( "New York", "Vermont", "Maine", "New Hampshire")) |> st_union() -> shp_NE # get CF trends and split into categories CF_trends |> filter(analysis == "923 period") |> mutate(trend_CFpp_per_decade = sens_slope * 10 * 100) |> select(COMPLXID, trend_CFpp_per_decade, p_value) -> trend_data_all_plants trend_data_all_plants |> filter(trend_CFpp_per_decade < 0 & p_value < 0.05) |> pull(COMPLXID) -> dams_sig_down trend_data_all_plants |> filter(trend_CFpp_per_decade > 0 & p_value < 0.05) |> pull(COMPLXID) -> dams_sig_up trend_data_all_plants |> filter(p_value > 0.05) |> pull(COMPLXID) -> dams_nonsig trend_data_all_plants |> filter(is.na(p_value)) |> pull(COMPLXID) -> dams_nodata st_read(EHA) |> mutate(COMPLXID = str_replace(EHA_PtID, "\\_..*", "")) |> filter(COMPLXID %in% CFs$COMPLXID) |> st_transform(st_crs(NA_shaded_relief)) |> mutate(size_cat = case_when( CH_MW > 1000 ~ "> 1 ", CH_MW > 500 & CH_MW <= 1000 ~ "0.5 - 1 ", CH_MW > 100 & CH_MW <= 500 ~ "0.1 - 0.5 ", CH_MW >= 5 & CH_MW <= 100 ~ "0.05 - 0.1 " )) |> mutate(size_cat = factor(size_cat, levels = c( "0.05 - 0.1 ", "0.1 - 0.5 ", "0.5 - 1 ", "> 1 " ))) |> left_join(trend_data_all_plants, by = join_by(COMPLXID)) |> filter(!State %in% c("AK", "HI")) -> dam_points filter(dam_points, COMPLXID %in% dams_nonsig) -> dam_points_nosig filter(dam_points, COMPLXID %in% dams_sig_down) -> dam_points_sig_down filter(dam_points, COMPLXID %in% dams_sig_up) -> dam_points_sig_up filter(dam_points, COMPLXID %in% dams_nodata) -> dam_points_nodata # main map plot ggplot() + geom_stars(data = NA_shaded_relief |> st_crop(lower48_shp), show.legend = FALSE) + scale_fill_gradient(low = "grey70", high = "grey95", na.value = NA) + coord_equal() + geom_sf(data = lower48_shp, fill = NA, col = "grey20") + geom_sf(data = shp_CA, col = "black",fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = shp_PNW, col = "black", fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = shp_SE, col = "black", fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = shp_NE, col = "black", fill = "limegreen", lwd = 0.8, alpha = 0.2) + geom_sf(data = major_rivers_US, col = "#0070B9", lwd = 0.05) + theme_void() + new_scale_fill() + geom_sf(data = dam_points_nodata, aes(size = size_cat), shape = 21, col = "black", fill = NA) + geom_sf(data = dam_points_nosig, aes(size = size_cat), shape = 21, col = "black", alpha = 0.8, fill = "white") + geom_sf(data = dam_points_sig_down, aes(size = size_cat, fill = trend_CFpp_per_decade), shape = 21, col = "black", alpha = 0.8) + geom_sf(data = dam_points_sig_up, aes(size = size_cat), shape = 21, col = "black", alpha = 0.8, fill = "dodgerblue") + scale_fill_viridis(option = "A") + scale_size_manual(values = c( "> 1 " = 7, "0.5 - 1 " = 5, "0.1 - 0.5 " = 2.5, "0.05 - 0.1 " = 1 )) + theme(legend.position = c(0.5,0.99), legend.direction = "horizontal", legend.box = "vertical", legend.background = element_rect(fill = "white", color = "white")) + labs( fill = "Trend in Capacity Factor (PPPD) \n (neg. trend w/ p<0.05 only)", size = "Nameplate in 2021 (GW) ", title = "a") + guides(size = guide_legend(override.aes=list(fill=NA), order = 1)) + annotate(geom = "point", size = 5, x = -90e5, y = 665e4, fill = "dodgerblue", col = "black", alpha = 0.8, pch = 21) + annotate("text", label = "pos. trend (p<0.05)", x = -89e5, y = 665e4, hjust = 0, size = 3) + annotate(geom = "point", size = 5, x = -90e5, y = 635e4, fill = "white", alpha = 0.8, pch = 21, col = "black") + annotate("text", label = "no significant trend", x = -89e5, y = 635e4, hjust = 0, size = 3) + annotate(geom = "point", size = 5, x = -90e5, y = 605e4, fill = "grey80", pch = 21, col = "black") + annotate("text", label = "insufficient data", x = -89e5, y = 605e4, hjust = 0, size = 3) + annotate("segment", x = -141e5, xend = -141e5, y = 560e4, yend = 240e4) + annotate("segment", x = -141e5, xend = -139e5, y = 560e4, yend = 560e4) + annotate("segment", x = -130e5, xend = -122e5, y = 375e4, yend = 240e4) + annotate("segment", x = -77e5, xend = -77e5, y = 535e4, yend = 240e4) + annotate("segment", x = -97e5, xend = -97e5, y = 345e4, yend = 240e4) + NULL -> map ### CA #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State == "CA")$COMPLXID) -> dam_gen_CA dam_gen_CA |> filter(!is.na(gen_MWh), !is.na(cap_MWh)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_CA senSlope(CF ~ year, data = CF_CA, intercept = "A1m") -> trend_CA CF_CA |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_CA$coefficients[2], intercept = trend_CA$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "c", subtitle = "[CF] California") -> p1_CA #### ### PNW #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State %in% c("WA", "OR", "ID"))$COMPLXID) -> dam_gen_PNW dam_gen_PNW |> filter(!is.na(gen_MWh), !is.na(cap_MWh)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_PNW senSlope(CF ~ year, data = CF_PNW, intercept = "A1m") -> trend_PNW CF_PNW |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_PNW$coefficients[2], intercept = trend_PNW$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "b", subtitle = "[CF] Pacific NW") -> p1_PNW #### ### SE #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State %in% c("TN", "AL", "GA", "SC"))$COMPLXID) -> dam_gen_SE dam_gen_SE |> filter(!is.na(gen_MWh), !is.na(cap_MWh)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_SE senSlope(CF ~ year, data = CF_SE, intercept = "A1m") -> trend_SE CF_SE |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_SE$coefficients[2], intercept = trend_SE$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "d", subtitle = "[CF] Southeast") -> p1_SE #### ### NE #### CFs |> filter(year %in% 1980:2021, COMPLXID %in% filter(dam_points, State %in% c("NY", "VT", "NH", "ME"))$COMPLXID) -> dam_gen_NE dam_gen_NE |> filter(!is.na(gen_MWh), !is.na(cap_MWh), !is.na(CF)) |> summarise(gen = sum(gen_MWh), cap = sum(cap_MWh), CF_75 = quantile(CF, 0.75), CF_25 = quantile(CF, 0.25), .by = year) |> mutate(CF = gen / cap) -> CF_NE senSlope(CF ~ year, data = CF_NE, intercept = "A1m") -> trend_NE CF_NE |> ggplot(aes(year, CF)) + geom_ribbon(aes(ymax = CF_75, ymin = CF_25), fill = "grey90") + geom_line(lwd = 1) + theme_classic() + geom_abline(slope = trend_NE$coefficients[2], intercept = trend_NE$coefficients[1], lwd = 1) + geom_hline(yintercept = 1, linetype = 2) + geom_hline(yintercept = 0) + scale_y_continuous(expand = c(0, 0), limits = c(0, 1), breaks = c(0,1)) + scale_x_continuous(breaks = c(1980, 2000, 2020)) + theme(#axis.line = element_blank(), axis.text.y = element_text(size = 10), #axis.ticks.y = element_blank(), axis.title = element_blank(), axis.text.x = element_text(size = 10)) + labs(title = "e", subtitle = "[CF] Northeast") -> p1_NE #### # combine plots and save as pdf (p1_PNW + p1_CA + p1_SE + p1_NE + plot_layout(ncol = 4)) -> low_plots map + plot_spacer() + low_plots + plot_layout(ncol = 1, heights = c(0.75, -0.1 ,0.25)) -> final_plot ggsave( "./output/plots/CF_trend_map.pdf", final_plot, height = 8, width = 8 ) ggsave( "./output/plots/CF_trend_map.png", final_plot, height = 8, width = 8 ) return(final_plot) } plot_main_driver_of_trend <- function(CFs, CF_trends, EHA){ flow_CF_models_scored |> filter(flow_type == "c") |> filter(rsq > 0.2, abs(bias) < 0.1) |> arrange(-rsq) |> pull(COMPLXID) -> viable_flow_models CF_trends |> filter(analysis == "flow period", !is.na(p_value)) |> filter(p_value < 0.05) -> dams_with_significant_change flow_CF_models |> filter(COMPLXID %in% viable_flow_models) |> select(COMPLXID, year, CF, flow_type, CF_pred) |> pivot_wider(names_from = "flow_type", values_from = "CF_pred") |> split(~COMPLXID) |> map_dfr(function(x){ smwrStats::senSlope(CF ~ year, x) -> sen_CF smwrStats::senSlope(c ~ year, x) -> sen_flowC smwrStats::senSlope(n ~ year, x) -> sen_flowN bind_rows( tibble( slope = sen_CF$coefficients[2], data = "CF" ), tibble( slope = sen_flowC$coefficients[2], data = "Assimilated flow" ), tibble( slope = sen_flowN$coefficients[2], data = "Naturalized flow" ) ) |> mutate(COMPLXID = first(x$COMPLXID)) }) -> all_slopes all_slopes |> pivot_wider(names_from = "data", values_from = slope) |> mutate(explained_by_water = `Assimilated flow` / CF, explained_by_climate = `Naturalized flow` / CF) |> select(COMPLXID, CF_trend = CF, explained_by_water, explained_by_climate) |> mutate(water_not_climate = if_else( explained_by_water > 0.5 & explained_by_climate < 0.3, TRUE, FALSE)) |> mutate( flow_attrib = case_when( explained_by_water >= 0.8 ~ "Strongly driven by water availability", explained_by_water < 0.8 & explained_by_water >= 0.6 ~ "Mostly driven by water availability", explained_by_water < 0.6 & explained_by_water >= 0.4 ~ "~ Half of trend driven by water availability", explained_by_water < 0.4 & explained_by_water >= 0.2 ~ "Mostly driven by other factors", explained_by_water < 0.2 ~ "Strongly driven by other factors" )) |> select(COMPLXID, water_not_climate, flow_attrib) -> flow_cats dams_with_significant_change |> left_join(flow_cats, by = join_by(COMPLXID)) |> filter(sens_slope > 0) -> positive_cases dams_with_significant_change |> left_join(flow_cats, by = join_by(COMPLXID)) |> filter(sens_slope < 0) -> negative_cases st_read(EHA) |> mutate(COMPLXID = str_replace(EHA_PtID, "\\_..*", "")) |> filter(COMPLXID %in% dams_with_significant_change$COMPLXID) |> st_transform(st_crs(NA_shaded_relief)) |> mutate(size_cat = case_when( CH_MW > 1000 ~ "> 1 ", CH_MW > 500 & CH_MW <= 1000 ~ "0.5 - 1 ", CH_MW > 100 & CH_MW <= 500 ~ "0.1 - 0.5 ", CH_MW >= 5 & CH_MW <= 100 ~ "0.05 - 0.1 " )) |> mutate(size_cat = factor(size_cat, levels = c( "0.05 - 0.1 ", "0.1 - 0.5 ", "0.5 - 1 ", "> 1 " ))) |> filter(!State %in% c("AK", "HI")) -> dam_points dam_points |> filter(COMPLXID %in% positive_cases$COMPLXID) |> left_join(positive_cases) -> dam_points_positive dam_points |> filter(COMPLXID %in% negative_cases$COMPLXID) |> left_join(negative_cases) -> dam_points_negative ## map backgrounds ## # plot background (greyscale shaded relief) read_stars("./data/Natural Earth/US_MSR_10M/US_MSR.tif", RasterIO = list(nBufXSize = 600, nBufYSize = 600)) -> NA_shaded_relief # major rivers intersecting the US st_read("./data/Natural Earth/10m_physical/ne_10m_rivers_lake_centerlines.shp") |> st_transform(st_crs(NA_shaded_relief)) -> major_rivers major_rivers[unlist(st_intersects(lower48_shp, major_rivers)),] -> major_rivers_US # US state boundaries st_read("./data/Natural Earth/10m_cultural/ne_10m_admin_1_states_provinces_lakes.shp") |> filter(admin == "United States of America", !name_en %in% c("Alaska", "Hawaii")) |> st_transform(st_crs(NA_shaded_relief)) -> lower48_shp ggplot() + geom_stars(data = NA_shaded_relief |> st_crop(lower48_shp), show.legend = FALSE) + scale_fill_gradient(low = "grey70", high = "grey95", na.value = NA) + coord_equal() + geom_sf(data = lower48_shp, fill = NA, col = "grey20") + geom_sf(data = major_rivers_US, col = "#0070B9", lwd = 0.05) + theme_void() + new_scale_fill() + geom_sf(data = dam_points_negative, aes(fill = flow_attrib, size = size_cat), shape = 25) + geom_sf(data = dam_points_positive, aes(fill = flow_attrib, size = size_cat), shape = 24) + scale_size_manual(values = c( "> 1 " = 7, "0.5 - 1 " = 5, "0.1 - 0.5 " = 2.5, "0.05 - 0.1 " = 1 )) + scale_fill_manual(values = c( "Strongly driven by water availability" = "#440154FF", "Mostly driven by water availability" = "#3B528BFF", "~ Half of trend driven by water availability" = "#21908CFF", "Mostly driven by other factors" = "#5DC863FF", "Strongly driven by other factors" = "#FDE725FF" )) + NULL }
_targets.R +24 −1 Original line number Diff line number Diff line Loading @@ -14,10 +14,16 @@ tar_option_set( "furrr", # ^^ apply in parallel "future", # parallel processing "sf", # geospatial tools "stars", # raster tools "foreign", # read old dbf data "zoo", # interpolation functions "smwrStats", # Sen's slope computation "stringr" # string manipulation functions "stringr", # string manipulation functions "terrainr", # map plotting "ggplot2", # plotting "ggnewscale",# plot support "patchwork", # plot combining "viridis" # color palettes ) ) Loading Loading @@ -104,5 +110,22 @@ list( CFs = dam_annual_gen_cap_CF_1970_2021, flows = DayFlow), format = "parquet" ), tar_target( flow_CF_models, simulate_annual_CF(data = model_data_ready), format = "parquet" ), tar_target( flow_CF_models_scored, evaluate_models(model_results = flow_CF_models), format = "parquet" ), tar_target( CF_trend_map, plot_CF_trend_map(CFs = dam_annual_gen_cap_CF_1970_2021, CF_trends = CF_trends, EHA = EHA), format = "rds" ) )
