Merge pull request #30 from afsc-gap-products/dev

Added historical temperature data

.gitignore

@@ -50,3 +50,4 @@ vignettes/*.pdf
 /Presentation/
 /output/
 *.tiff
+/plots/cpchange

0_update_cold_pool_index.Rmd

@@ -231,7 +231,7 @@ if(update_sysdata) {
                                                         file_name_contains = "ste_",
                                                         file_type = ".tif")
   
-  cpa_pre2021 <- coldpool:::cpa_pre2021
+  cpa_pre2021 <- read.csv(file = here::here("inst", "extdata", "old_method_cpa_temperature_2021.csv"))
   ebs_proj_crs <- coldpool:::ebs_proj_crs
   cold_pool_index <- output_df
     

1_cold_pool_index.Rmd

@@ -609,6 +609,74 @@ print(plot_sebs_average_temperature)
 dev.off()
 ```
 
+```{r zscore_temperature}
+cp_summary <- dplyr::bind_rows(
+  dplyr::mutate(cold_pool_index, 
+                diff = MEAN_GEAR_TEMPERATURE - mean(MEAN_GEAR_TEMPERATURE)) |> 
+    dplyr::mutate(sign = sign(diff),
+                  z = diff/sd(MEAN_GEAR_TEMPERATURE)) |>
+    dplyr::inner_join(data.frame(sign = c(-1,1),
+                                 symbol = c("-","+"),
+                                 col = c(1,2))) |>
+    dplyr::select(YEAR, diff, symbol, z, col) |>
+    dplyr::mutate(var = "Bottom temperature"),
+  dplyr::mutate(cold_pool_index, 
+                diff = MEAN_SURFACE_TEMPERATURE - mean(MEAN_SURFACE_TEMPERATURE)) |> 
+    dplyr::mutate(sign = sign(diff),
+                  z = diff/sd(MEAN_SURFACE_TEMPERATURE)) |>
+    dplyr::inner_join(data.frame(sign = c(-1,1),
+                                 symbol = c("-","+"),
+                                 col = c(1,2))) |>
+    dplyr::select(YEAR, diff, symbol, z, col) |>
+    dplyr::mutate(var = "Sea surface temperature"),
+  dplyr::mutate(cold_pool_index, 
+                diff = AREA_LTE2_KM2 - mean(AREA_LTE2_KM2)) |> 
+    dplyr::mutate(sign = sign(diff),
+                  z = diff/sd(AREA_LTE2_KM2)) |>
+    dplyr::inner_join(data.frame(sign = c(-1,1),
+                                 symbol = c("-","+"),
+                                 col = c(2,1))) |>
+    dplyr::select(YEAR, diff, symbol, z, col) |>
+    dplyr::mutate(var = "Cold pool area")) |>
+  dplyr::mutate(group = YEAR < 2020,
+                var = factor(cp_summary$var, 
+                             levels = c("Bottom temperature",
+                                        "Sea surface temperature",
+                                        "Cold pool area")))
+
+zscore_plot <- ggplot(data = cp_summary, 
+       aes(x = YEAR, 
+           y = z, 
+           group = group)) +
+  geom_hline(yintercept = c(0), linetype = 1) +
+  geom_hline(yintercept = c(-1,1), linetype = 2) +
+  geom_hline(yintercept = c(-2,2), linetype = 3) +
+  geom_point() +
+  geom_line() +
+  geom_text(data = cp_summary,
+            aes(x = YEAR, y = 2.5, label = symbol, color = factor(col))) +
+  facet_wrap(~var, nrow = 3) +
+  scale_x_continuous(name = "Year") +
+  scale_y_continuous(name = "Anomaly") +
+  scale_color_manual(values = c("blue", "red")) +
+  theme_bw() +
+  theme(axis.text = element_text(color = "black"),
+        axis.ticks = element_line(color = "black"),
+        panel.border = element_rect(color = "black", fill = NA),
+        panel.background = element_rect(color = "black", fill = NA),
+        strip.text = element_text(size = 9,
+                                  color = "white",
+                                  face = "bold",
+                                  margin = margin(0.5, 0, 0.5, 0, "mm")),
+        strip.background = element_rect(fill = "#0055a4",
+                                        color = NA),
+        legend.position = "none")
+
+png(file = here::here("plots", paste0(max_year, "_anomaly.png")), width = 6, height = 6, units = "in", res = fig_res)
+print(zscore_plot)
+dev.off()
+```
+
 Contributed by Sean Rohan^1^ and Lewis Barnett^1^   
 ^1^ Resource Assessment and Conservation Engineering Division, Alaska Fisheries Science Center, National Marine Fisheries Service, NOAA  
 **Contact**: sean.rohan@noaa.gov

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: coldpool
 Type: Package
 Title: Generate GAP's EBS temperature products
-Version: 1.3
+Version: 1.4
 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")))

NAMESPACE

@@ -3,6 +3,7 @@
 export(cold_pool_index)
 export(compare_cpa_station_filter)
 export(cpa_from_raster)
+export(cpa_pre2021)
 export(ebs_bottom_temperature)
 export(ebs_proj_crs)
 export(ebs_surface_temperature)

R/data.R

@@ -35,21 +35,29 @@
 
 #' Historical cold pool area calculations (cpa_pre2021)
 #' 
-#' Historical temperature cold pool and temperature band areas that were calculated prior to 2021 by Bob Lauth and Lyle Britt. Historical cold pool calculations were conducted in ArcGIS and used inverse distance weighting with a maximum of four nearest neighbor points for interpolation (search radius nmax = 4) and weighting function power equal to 2 (i.e., idp = 2). Default ArcGIS settings were used to generate rasters cells for interpolation (cell dimensions equal to the smaller horizontal or vertical dimension divided by 250). Protocols for including stations for interpolation varied among years but it was intended for all standard survey grid stations from the EBS survey to be included. In some years, it is possible that near-shore stations or red king crab resample stations were included in the calculations. Rasters were converted to polygon shapefiles and areas were calculated from polygons. The CRS used in ArcGIS .mxd files was North American Datum 1983 / Albers Equal Area Alaska (EPSG:3338)
+#' Historical temperature cold pool area  (AREA_SUM_KM2_LTE2) and other temperature isotherm (AREA_KM2_MINUS1, AREA_KM2_0, AREA_KM2_1, AREA_KM2_2, AREA_KM2_LTE0, AREA_KM2_LTE1) areas that were calculated prior to 2021 by Bob Lauth, Lyle Britt, Dan Nichol, and Rebecca Haehn. These calculations were obtained from raster surfaces generated using ArcMap by interpolating temperatures using inverse distance weighting with a maximum of four nearest neighbor points for interpolation (search radius nmax = 4) and weighting function power equal to 2 (i.e., idp = 2). Default ArcGIS settings were used to generate rasters cells for interpolation (cell dimensions equal to the smaller horizontal or vertical dimension divided by 250). Protocols for including stations for interpolation varied among years but it was intended for all standard survey grid stations from the EBS survey to be included. In some years, it is possible that near-shore stations or red king crab resample stations were included in the calculations. Rasters were converted to polygon shapefiles and areas were calculated from polygons. The CRS used in ArcGIS .mxd files was North American Datum 1983 / Albers Equal Area Alaska (EPSG:3338). Historical weighted average temperature data products (AVGBSBT_PLUSNW, AVGBSST_PLUSNW, AVGBSBT_STANDARD, AVGBSST_STANDARD, AVGBSBT_NBS, AVGBSST_NBS, AVGBSBT_LT100M, AVGBSST_LT100M) were calculated by weighting temperature observations from survey stations in proportion to the area of the stratum where they were collected relative to the total survey area. Historical mean temperature for stations at bottom depths < 100 m (AVGBSBT_LT100M, AVGBSST_LT100M) were calculated by averaging temperature observations from samples collected at bottom depths < 100 m.
 #' 
-#' @format A data frame with 38 rows and 9 columns:
+#' @format A data frame with 39 rows and 16 columns:
 #' \describe{
 #'      \item{YEAR}{Year}
-#'      \item{AVG_STRATA_WEIGHTED_BOTTEMP_STD_AREA}{Mean bottom temperature by stratum, weighted by stratum area, in square kilometers}
 #'      \item{AREA_KM2_MINUS1}{Area with bottom temperatures less than -1 celsius, in square kilometers}
 #'      \item{AREA_KM2_0}{Area with bottom temperatures between -1 and 0 celsius, in square kilometers}
 #'      \item{AREA_KM2_1}{Area with bottom temperatures between 0 and 1 celsius, in square kilometers}
 #'      \item{AREA_KM2_2}{Area with bottom temperatures between 1 and 2 celsius, in square kilometers}
 #'      \item{AREA_SUM_KM2_LTE2}{Cold Pool Index. Total area with bottom temperatures less than or equal to 2 celsius, in square kilometers}
 #'      \item{AREA_KM2_LTE1}{Total area with bottom temperatures less than or equal to 1 celsius, in square kilometers}
-#'       \item{AREA_KM2_LTE0}{Total area with bottom temperatures less than or equal to 0 celsius, in square kilometers}
+#'      \item{AREA_KM2_LTE0}{Total area with bottom temperatures less than or equal to 0 celsius, in square kilometers}
+#'      \item{AVGBSBT_PLUSNW}{Stratum-area weighted mean bottom temperature for EBS shelf standard and plusNW survey strata}
+#'      \item{AVGBSST_PLUSNW}{Stratum-area weighted mean sea surface temperature for EBS shelf standard and plusNW survey strata}
+#'      \item{AVGBSBT_STANDARD}{Stratum-area weighted mean bottom temperature for EBS shelf standard survey strata}
+#'      \item{AVGBSST_STANDARD}{Stratum-area weighted mean sea surface temperature for EBS shelf standard survey strata}
+#'      \item{AVGBSBT_NBS}{Stratum-area weighted mean bottom temperature for the NBS}
+#'      \item{AVGBSST_NBS}{Stratum-area weighted mean sea surface temperature for the NBS}
+#'      \item{AVGBSBT_LT100M}{Mean bottom temperature for EBS shelf stations at bottom depths < 100 m}
+#'      \item{AVGBSST_LT100M}{Mean sea surface temperature for EBS shelf stations at bottom depths < 100 m}
 #' }
 #' @source \url{https://www.fisheries.noaa.gov/contact/groundfish-assessment-program}
+#' @export
 "cpa_pre2021"
 
 #' CRS for eastern Bering Sea cold pool index

R/figs.R

@@ -28,7 +28,8 @@ plot_stn_doy <- function() {
                           dplyr::summarise()) %>%
     dplyr::inner_join(haul_dat %>%
                         dplyr::group_by(STATIONID) %>%
-                        dplyr::summarise(MEAN_DOY = mean(DOY)))
+                        dplyr::summarise(MEAN_DOY = mean(DOY),
+                                         MEAN_DEPTH = mean(BOTTOM_DEPTH, na.rm = TRUE)))
 
   # Define plot exent (through trial end error)
   panel_extent <- data.frame(x = c(-1326559.21, -87636.05),
@@ -113,8 +114,84 @@ plot_stn_doy <- function() {
                    legend.position = "right",
                    legend.background = element_blank())
 
+  doy_for_tm <- ggplot() +
+    geom_sf(data = ebs_layers$akland, 
+            fill = "grey70", 
+            color = "black") +
+    geom_sf(data = start_df, 
+            aes(fill = MEAN_DOY)) +
+    geom_sf_text(data = start_df, 
+            aes(label = round(MEAN_DOY))) +
+    ggplot2::geom_sf(data = agg_stratum,
+                     fill = NA,
+                     color = "black",
+                     size = rel(1.1)) +
+    shadowtext::geom_shadowtext(data = data.frame(agg_stratum = sf::st_centroid(agg_stratum)$agg_stratum, # Centroid of aggregate stratum polygons
+                                                  x = sf::st_coordinates(sf::st_centroid(agg_stratum))[,1], # Coordinates of centroid of aggregate stratum polygons
+                                                  y = sf::st_coordinates(sf::st_centroid(agg_stratum))[,2]),# Coordinates of centroid of aggregate stratum polygons 
+                                aes(x = x,
+                                    y = y, 
+                                    label = agg_stratum),
+                                size = rel(4.5),
+                                color = "black",
+                                bg.color = "white") +
+    shadowtext::geom_shadowtext(data = data.frame(x = -158.5, 
+                                                  y = 62.4, 
+                                                  lab = "Alaska") %>%
+                                  akgfmaps::transform_data_frame_crs(out.crs = coldpool:::ebs_proj_crs),
+                                mapping = aes(x = x,
+                                              y = y,
+                                              label = lab),
+                                size = rel(6),
+                                color = "black",
+                                bg.color = "white") +
+    shadowtext::geom_shadowtext(data = data.frame(x = -166.2, 
+                                                  y = 60.08, 
+                                                  lab = "Nunivak\nIsland") %>%
+                                  akgfmaps::transform_data_frame_crs(out.crs = coldpool:::ebs_proj_crs),
+                                mapping = aes(x = x,
+                                              y = y,
+                                              label = lab),
+                                size = rel(2.5),
+                                color = "black",
+                                bg.color = "white") +
+    shadowtext::geom_shadowtext(data = data.frame(x = c(-169, -159.5), 
+                                                  y = c(66.3, 57.7), 
+                                                  lab = c("Bering\nStrait", "Bristol\nBay")) %>%
+                                  akgfmaps::transform_data_frame_crs(out.crs = coldpool:::ebs_proj_crs),
+                                mapping = aes(x = x,
+                                              y = y,
+                                              label = lab),
+                                size = rel(3),
+                                color = "black",
+                                bg.color = "white") +
+    scale_fill_viridis_c(name = "Sample Mean\nDay of Year", 
+                         option = "B") +
+    scale_color_brewer() +
+    coord_sf(xlim = panel_extent$x,
+             ylim = panel_extent$y) +
+    ggplot2::scale_x_continuous(name = "Longitude", 
+                                breaks = ebs_layers$lon.breaks) + 
+    ggplot2::scale_y_continuous(name = "Latitude", 
+                                breaks = ebs_layers$lat.breaks) +
+    theme_bw() +
+    ggplot2::theme(axis.title = element_blank(),
+                   panel.border = element_rect(color = "black", fill = NA),
+                   legend.title = element_text(size = 9, color = "black"),
+                   legend.text = element_text(size = 8, color = "black"),
+                   axis.text = element_text(size = 8, color = "black"),
+                   panel.grid = element_blank(),
+                   panel.background = element_rect(color = "black", fill = "#bee8ff"),
+                   legend.margin = margin(-12,0,0,0),
+                   legend.position = "right",
+                   legend.background = element_blank())
+
   png(filename = here::here("plots", "ebs_nbs_survey_area.png"), width = 5, height = 5, units = "in", res = 600)
   print(plot_ebs_nbs_survey_stations)
   dev.off()
 
+  png(filename = here::here("plots", "doy_for_tm.png"), width = 10, height = 10, units = "in", res = 120)
+  print(doy_for_tm)
+  dev.off()
+
 }