Skip to content

Commit

Permalink
Merge pull request #33 from afsc-gap-products/dev
Browse files Browse the repository at this point in the history 
Update with 2022 EBS temperature data and preliminary NBS data
  • Loading branch information
@sean-rohan-NOAA
sean-rohan-NOAA authored Sep 9, 2022
2 parents 1804a78 + 94b6448 commit 951e2e9
Show file tree
Hide file tree
Showing 68 changed files with 15,815 additions and 14,673 deletions.
16 changes: 11 additions & 5 deletions 0_update_cold_pool_index.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ ebs_csv_path <- here::here("data", paste0("index_hauls_temperature_data.csv"))
nbs_ebs_csv_path <- here::here("data", paste0("ebs_nbs_temperature_full_area.csv"))
### UPDATE!!! ###
nbs_years <- c(2010, 2017, 2018, 2019, 2021)
nbs_years <- c(2010, 2017, 2018, 2019, 2021, 2022)
```

# 2. Retrieve temperature data
Expand All @@ -51,7 +51,7 @@ if(update_sysdata) {
channel <- get_connected()
# Get temperature data and write csvs to data directory ----
coldpool::get_data(channel = channel)
coldpool:::get_data(channel = channel, include_preliminary_data = "nbs")
}
```

Expand All @@ -68,15 +68,15 @@ Writes:
interpolation_wrapper(temp_data_path = ebs_csv_path,
proj_crs = proj_crs,
cell_resolution = 5000, # 5x5 km grid resolution
select_years = 1982:2021,
select_years = 1982:2022,
interp_variable = "gear_temperature",
select_region = "sebs")
# Interpolate surface temperature and write rasters for SEBS
interpolation_wrapper(temp_data_path = ebs_csv_path,
proj_crs = proj_crs,
cell_resolution = 5000, # 5x5 km grid resolution
select_years = 1982:2021,
select_years = 1982:2022,
interp_variable = "surface_temperature",
select_region = "sebs")
```
Expand All @@ -86,7 +86,7 @@ interpolation_wrapper(temp_data_path = ebs_csv_path,
interpolation_wrapper(temp_data_path = nbs_ebs_csv_path,
proj_crs = proj_crs,
cell_resolution = 5000, # 5x5 km grid resolution
select_years = nbs_years,
select_years = nbs_years,
interp_variable = "gear_temperature",
select_region = "ebs") # Full EBS+NBS
Expand All @@ -110,6 +110,8 @@ Writes:
bottom_temp_files <- list.files(here::here("output", "raster", "sebs", "gear_temperature"),
full.names = TRUE)
bottom_temp_files <- bottom_temp_files[grep(pattern = "ste_", x = bottom_temp_files)]
bottom_temp_files <- bottom_temp_files[-grep(pattern = ".aux.xml", x = bottom_temp_files)]
bt_df <- data.frame(YEAR = numeric(length = length(bottom_temp_files)),
AREA_LTE2_KM2 = numeric(length = length(bottom_temp_files)),
AREA_LTE1_KM2 = numeric(length = length(bottom_temp_files)),
Expand Down Expand Up @@ -149,6 +151,8 @@ for(i in 1:length(bottom_temp_files)) {
surface_temp_files <- list.files(here::here("output", "raster", "sebs", "surface_temperature"),
full.names = TRUE)
surface_temp_files <- surface_temp_files[grep(pattern = "ste_", x = surface_temp_files)]
surface_temp_files <- surface_temp_files[-grep(pattern = ".aux.xml", x = surface_temp_files)]
sst_df <- data.frame(YEAR = numeric(length = length(bottom_temp_files)),
MEAN_SURFACE_TEMPERATURE = numeric(length = length(surface_temp_files)))
Expand All @@ -173,11 +177,13 @@ nbs_ebs_bt_files <- list.files(here::here("output", "raster", "ebs", "gear_tempe
# Select Stein's Matern and remove 2018 unplanned extension (didn't have full spatial coverage)
nbs_ebs_bt_files <- nbs_ebs_bt_files[grep(pattern = "ste_", x = nbs_ebs_bt_files)]
nbs_ebs_bt_files <- nbs_ebs_bt_files[-grep(pattern = 2018, x = nbs_ebs_bt_files)]
nbs_ebs_bt_files <- nbs_ebs_bt_files[-grep(pattern = ".aux.xml", x = nbs_ebs_bt_files)]
nbs_ebs_sst_files <- list.files(here::here("output", "raster", "ebs", "surface_temperature"),
full.names = TRUE)
nbs_ebs_sst_files <- nbs_ebs_sst_files[grep(pattern = "ste_", x = nbs_ebs_sst_files)]
nbs_ebs_sst_files <- nbs_ebs_sst_files[-grep(pattern = 2018, x = nbs_ebs_sst_files)]
nbs_ebs_sst_files <- nbs_ebs_sst_files[-grep(pattern = ".aux.xml", x = nbs_ebs_sst_files)]
nbs_mean_temperature <- data.frame(YEAR = numeric(length = length(nbs_ebs_bt_files)),
MEAN_GEAR_TEMPERATURE = numeric(length = length(nbs_ebs_bt_files)),
Expand Down
159 changes: 131 additions & 28 deletions 1_cold_pool_index.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ if(!dir.exists("plots")) {dir.create("plots")}
fig_res <- 600
cpa_palette <- c("#21dae7", "#0071ff", "#0000e3", "#000040")
max_year <- 2021
max_year <- 2022
min_year <- max_year - 19 # Default for 20 panel cold pool maps
coldpool:::cold_pool_index$YEAR[order(coldpool:::cold_pool_index$AREA_LTE2_KM2)]
Expand Down Expand Up @@ -100,28 +100,24 @@ area_df <- coldpool:::cold_pool_index %>%
proportion = value/sebs_layers$survey.area$AREA_KM2)
stacked_area_plot <- ggplot() +
geom_area(data = area_df %>%
dplyr::filter(YEAR < 2020),
geom_area(data = filter(area_df, YEAR < 2020),
mapping = aes(x = YEAR,
y = proportion,
fill = variable)) +
geom_area(data = filter(area_df, YEAR > 2020),
mapping = aes(x = YEAR,
y = proportion,
fill = variable)) +
geom_point(data = area_df %>%
dplyr::filter(YEAR == 2021),
mapping = aes(x = YEAR,
y = proportion,
fill = variable),
shape = 21,
size = rel(3)) +
scale_fill_manual(name = "Temperature",
values = cpa_palette) +
scale_y_continuous(name = "Proportion of EBS Shelf (Plus NW) Survey Area",
limits = c(0, 1),
expand = c(0, 0),
breaks = seq(0,1,0.1)) +
scale_x_continuous(name = "Year",
limits = c(1982, 2021.25),
limits = c(1982, max_year + 0.25),
expand = c(0, 0),
breaks = seq(1982,2021,4)) +
breaks = seq(1982, max_year, 4)) +
theme_bw() +
theme(panel.grid.minor = element_blank(),
axis.title = element_text(face = "bold"),
Expand All @@ -131,33 +127,29 @@ stacked_area_plot <- ggplot() +
legend.position = c(0.2, 0.8),
legend.title = element_blank())
png(file = here::here("plots", paste0(max_year, "stacked_temperature.png")), height = 4, width = 6, units = "in", res = 600)
png(file = here::here("plots", paste0(max_year, "_stacked_temperature.png")), height = 4, width = 6, units = "in", res = 600)
print(stacked_area_plot)
dev.off()
stacked_area_simple_label_plot <- ggplot() +
geom_area(data = area_df %>%
dplyr::filter(YEAR < 2020),
geom_area(data = filter(area_df, YEAR < 2020),
mapping = aes(x = YEAR,
y = proportion,
fill = variable)) +
geom_point(data = area_df %>%
dplyr::filter(YEAR == 2021),
y = proportion,
fill = variable)) +
geom_area(data = filter(area_df, YEAR > 2020),
mapping = aes(x = YEAR,
y = proportion,
fill = variable),
shape = 21,
size = rel(3)) +
y = proportion,
fill = variable)) +
scale_fill_manual(name = "Temperature",
values = cpa_palette) +
scale_y_continuous(name = "Proportion of EBS Shelf Survey Area",
limits = c(0, 1),
expand = c(0, 0),
breaks = seq(0,1,0.1)) +
scale_x_continuous(name = "Year",
limits = c(1982, 2021.25),
limits = c(1982, max_year + 0.25),
expand = c(0, 0),
breaks = seq(1982,2021,4)) +
breaks = seq(1982, max_year,4)) +
theme_bw() +
theme(panel.grid.minor = element_blank(),
axis.title = element_text(face = "bold"),
Expand Down Expand Up @@ -358,7 +350,7 @@ panel_extent <- data.frame(y = c(53, 67),
akgfmaps::transform_data_frame_crs(out.crs = coldpool:::ebs_proj_crs)
nbs_ebs_temp_breaks <- c(-Inf, seq(-1,8,1), Inf)
nbs_ebs_viridis_option <- "H" # viridis turbo palette
nbs_ebs_viridis_option <- "C" # viridis turbo palette
n_temp_breaks <- length(nbs_ebs_temp_breaks)-1
ebs_nbs_temperature_map <- ggplot2::ggplot() +
Expand Down Expand Up @@ -450,7 +442,7 @@ ebs_nbs_temperature_contour_map <- ggplot2::ggplot() +
plot.margin = unit(c(5.5, 5.5,-25,5.5), units = "pt"))
temp_map_cbar <- coldpool::legend_discrete_cbar(breaks = nbs_ebs_temp_breaks,
colors = viridis::viridis_pal(option = "H")(n_temp_breaks),
colors = viridis::viridis_pal(option = nbs_ebs_viridis_option)(n_temp_breaks),
legend_direction = "horizontal",
font_size = 3,
width = 0.1,
Expand Down Expand Up @@ -486,6 +478,117 @@ dev.off()
png(filename = here::here("plots", paste0(max_year, "_nbs_ebs_temperature_contour_map.png")), width = 6, height = 6, units = "in", res = fig_res)
print(ebs_nbs_contour_map)
dev.off()
# Four panel bottom temperature map
max_year <- 2022
skip_year <- 2020
n_panels <- 4
grid_layout <- c(2,2)
nbs_years <- c(2010, 2017, 2019, 2021, 2022)
plot_years <- (max_year-n_panels+ifelse(max_year-n_panels < skip_year, 0, 1)):max_year
plot_years <- plot_years[!(plot_years %in% skip_year)]
ebs_years <- which(!(plot_years %in% nbs_years))
fig_res <- 600
nbs_ebs_layers <- akgfmaps::get_base_layers(select.region = "ebs",
set.crs = coldpool:::ebs_proj_crs)
coords <- raster::coordinates(coldpool:::nbs_ebs_bottom_temperature)
bt_year_df <- data.frame()
for(ii in 1:length(plot_years)) {
if(plot_years %in% ebs_years) {
sel_raster <- raster::mask(coldpool:::nbs_ebs_bottom_temperature,
nbs_ebs_layers$survey.area |>
filter(SURVEY == "EBS_SHELF"))
} else {
sel_raster <- coldpool:::nbs_ebs_bottom_temperature
}
sel_layer_df <- data.frame(x = coords[,1],
y = coords[,2],
temperature = sel_raster@data@values[,grep(pattern = plot_years[ii], x = names(sel_raster))])
sel_layer_df <- sel_layer_df[!is.na(sel_layer_df$temperature),]
sel_layer_df$year <- plot_years[ii]
bt_year_df <- dplyr::bind_rows(bt_year_df, sel_layer_df)
}
# Union to combine strata 31, 32 into 30, etc.
nbs_ebs_agg_strata <- nbs_ebs_layers$survey.strata %>%
dplyr::mutate(agg_stratum = Stratum) %>%
dplyr::mutate(agg_stratum = replace(agg_stratum, agg_stratum %in% c(31,32), 30),
agg_stratum = replace(agg_stratum, agg_stratum %in% c(41,42,43), 40),
agg_stratum = replace(agg_stratum, agg_stratum %in% c(61,62), 60)) %>%
dplyr::group_by(agg_stratum) %>%
dplyr::summarise()
nbs_ebs_temp_breaks <- c(-Inf, seq(-1,8,1), Inf)
nbs_ebs_viridis_option <- "B" # viridis turbo palette
ebs_nbs_temperature_map <- ggplot2::ggplot() +
ggplot2::geom_sf(data = nbs_ebs_layers$akland,
fill = "grey70",
color = "black") +
ggplot2::geom_sf(data = nbs_ebs_layers$survey.area, fill = "grey65") +
ggplot2::geom_tile(data = bt_year_df,
aes(x = x,
y = y,
fill = cut(temperature, breaks = nbs_ebs_temp_breaks))) +
ggplot2::facet_wrap(~year, ncol = grid_layout[1], nrow = grid_layout[2]) +
ggplot2::geom_sf(data = nbs_ebs_agg_strata,
fill = NA,
color = "black") +
ggplot2::geom_sf(data = nbs_ebs_layers$graticule,
alpha = 0.3) +
ggplot2::coord_sf(xlim = nbs_ebs_layers$plot.boundary$x,
ylim = nbs_ebs_layers$plot.boundary$y) +
ggplot2::scale_x_continuous(name = "Longitude",
breaks = nbs_ebs_layers$lon.breaks) +
ggplot2::scale_y_continuous(name = "Latitude",
breaks = nbs_ebs_layers$lat.breaks) +
ggplot2::scale_fill_manual(values = viridis_pal(option = nbs_ebs_viridis_option)(length(nbs_ebs_temp_breaks)-1)) +
coldpool::theme_multi_map_blue_strip() +
theme(legend.position = "none",
plot.margin = unit(c(5,5,-5,5), units = "mm"),
axis.title = element_blank(),
axis.text = element_text(size = 9))
temp_map_cbar <- coldpool::legend_discrete_cbar(breaks = nbs_ebs_temp_breaks,
colors = viridis::viridis_pal(option = nbs_ebs_viridis_option),
legend_direction = "horizontal",
font_size = 3,
width = 0.1,
expand_size.x = 0.3,
expand_size.y = 0.3,
expand.x = 0.3,
expand.y = 0.9,
spacing_scaling = 1.2,
text.hjust = 0.5,
font.family = "sans",
neat.labels = FALSE) +
annotate("text",
x = 1.15,
y = 3.5,
label = expression(bold("Bottom Temperature"~(degree*C))),
size = rel(3.2)) +
theme(plot.margin = unit(c(0,0, 0, 5), units = "mm"))
ebs_nbs_map_grid <- cowplot::plot_grid(ebs_nbs_temperature_map,
temp_map_cbar,
nrow = 2,
rel_heights = c(0.85,0.15))
png(filename = here::here("plots", paste0(max_year, "_nbs_ebs_temperature_map_grid.png")), width = 6, height = 6, units = "in", res = fig_res)
print(ebs_nbs_map_grid)
dev.off()
```

```{r fig4_mean_temperature, message=FALSE, warning=FALSE, include=FALSE}
Expand Down Expand Up @@ -639,7 +742,7 @@ cp_summary <- dplyr::bind_rows(
dplyr::select(YEAR, diff, symbol, z, col) |>
dplyr::mutate(var = "Cold pool area")) |>
dplyr::mutate(group = YEAR < 2020,
var = factor(cp_summary$var,
var = factor(var,
levels = c("Bottom temperature",
"Sea surface temperature",
"Cold pool area")))
Expand Down
4 changes: 2 additions & 2 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: coldpool
Type: Package
Title: Generate GAP's EBS temperature products
Version: 1.4
Version: 2.0
Authors@R: c(person("Sean", "Rohan", email = "sean.rohan@noaa.gov", role = c("aut", "cre")),
person("Lewis", "Barnett", email = "lewis.barnett@noaa.gov", role = c("aut", "ctb")),
person("Emily", "Markowitz", role = c("ctb")))
Expand All @@ -28,7 +28,7 @@ Description: This package calculates the area of the cold pool in the eastern Be
License: GPL (>=2)
Encoding: UTF-8
LazyData: false
RoxygenNote: 7.1.2
RoxygenNote: 7.2.1
Imports: colorspace,
viridis
Suggests: knitr, rmarkdown, RODBC, getPass, readr, testthat, cowplot, metR, cowplot
1 change: 1 addition & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

export(cold_pool_index)
export(compare_cpa_station_filter)
export(convert_ddm_to_dd)
export(cpa_from_raster)
export(cpa_pre2021)
export(ebs_bottom_temperature)
Expand Down
Loading

0 comments on commit 951e2e9

Please sign in to comment.