## ----global_options, include=FALSE-------------------------------------------- knitr::opts_chunk$set(fig.width = 8, fig.height = 8, fig.align = 'centre', echo = TRUE, warning = FALSE, message = FALSE, eval = FALSE, tidy = FALSE) ## ----load-pkg, eval=FALSE----------------------------------------------------- # library(heatwaveR) # library(lubridate) # library(dplyr) # library(tidyr) # library(purrr) # library(ggplot2) # library(viridis) # library(trend) # library(kableExtra) # library(rerddap) # library(sf) ## ----MEOW, eval=FALSE--------------------------------------------------------- # # Load and process the MEOW shapefile # # NB: Change your file path to where the files were unzipped # MEOW_NZ <- st_read('~/Desktop/meow_ecos.shp') %>% # dplyr:::filter(PROVINCE %in% c('Northern New Zealand', 'Southern New Zealand', # 'Subantarctic New Zealand')) %>% # st_make_valid() %>% # st_transform(crs = 4326) %>% # st_shift_longitude() %>% # mutate(ECOREGION = factor(ECOREGION, levels = c("Kermadec Island", "Three Kings-North Cape", # "Northeastern New Zealand", "Central New Zealand", # "Chatham Island", "South New Zealand", # "Bounty and Antipodes Islands", # "Snares Island","Auckland Island","Campbell Island"))) # # # Find max lon/lat ranges # lon_range <- range(sf::st_coordinates(MEOW_NZ$geometry)[,1]) # lat_range <- range(sf::st_coordinates(MEOW_NZ$geometry)[,2]) # # # Expand a bit to make sure all necessary pixels are downloaded # lon_range <- c(lon_range[1]-0.25, lon_range[2]+0.25) # lat_range <- c(lat_range[1]-0.25, lat_range[2]+0.25) ## ----MEOW_pix, eval=FALSE----------------------------------------------------- # # Create OISST grid # lon_lat_OISST <- base::expand.grid(seq(0.125, 359.875, by = 0.25), # seq(-89.875, 89.875, by = 0.25)) %>% # dplyr::rename(lon = Var1, lat = Var2) %>% # dplyr::arrange(lon, lat) %>% # st_as_sf(coords = c('lon', 'lat'), crs = 4326) # # # Join OISST grid to MEOW spatial and filter out pixels not within NZ EEZ # lon_lat_NZ <- st_join(lon_lat_OISST, MEOW_NZ) %>% # dplyr::select(geometry, ECOREGION, PROVINCE) %>% # dplyr::filter(!is.na(ECOREGION)) # # # Convert the coordinates back to a tibble for further use # # NB: Should be 6,398 pixels # lon_lat_NZ_df <- lon_lat_NZ %>% # mutate(lon = sf::st_coordinates(.)[,1], # lat = sf::st_coordinates(.)[,2]) %>% # mutate(lon = as.numeric(lon), lat = as.numeric(lat)) %>% # data.frame() %>% # dplyr::select(-geometry) ## ----load-data, eval=FALSE---------------------------------------------------- # # Get SST for around NZ + Islands and for all years # OISST_sub_dl <- function(time_df){ # OISST_dat <- griddap(x = "ncdcOisst21Agg", # url = "https://coastwatch.pfeg.noaa.gov/erddap/", # time = c(time_df$start, time_df$end), # zlev = c(0, 0), # latitude = lat_range, # longitude = lon_range, # fields = c("sst"))$data %>% # mutate(time = as.Date(stringr::str_remove(time, "T00:00:00Z"))) %>% # dplyr::rename(t = time, temp = sst) %>% # dplyr::select(lon, lat, t, temp) %>% # na.omit() # } # # # The span of years to download # dl_years <- data.frame(date_index = 1:5, # start = as.Date(c("1982-01-01", "1990-01-01", # "1998-01-01", "2006-01-01", "2014-01-01")), # end = as.Date(c("1989-12-31", "1997-12-31", # "2005-12-31", "2013-12-31", "2021-12-31"))) # # # Download the data # # NB: Takes ~21 minutes on a desktop computer # # NB: Contains 175,208,352 rows, 4 columns, and uses ~15 GB of RAM # OISST_data <- dl_years %>% # group_by(date_index) %>% # group_modify(~OISST_sub_dl(.x)) %>% # ungroup() %>% # dplyr::select(lon, lat, t, temp) # # # Filter out only the NZ pixels # # NB: Contains 87,195,152 rows # OISST_data_NZ <- left_join(lon_lat_NZ_df, OISST_data, by = c("lon", "lat")) %>% drop_na() # # # Get the climatologies only # clim_only <- function(df){ # # First calculate the climatologies # clim <- ts2clm(data = df, climatologyPeriod = c("1982-01-01", "2011-12-31")) # # Then the events # event <- detect_event(data = clim) # # Return only the climatology metric dataframe of results # return(event$climatology) # } # # # Extract the climatology values # # NB: Takes ~21 minutes on a desktop computer # OISST_clim <- OISST_data_NZ %>% # group_by(ECOREGION, PROVINCE, lon, lat) %>% # group_modify(~clim_only(.x)) %>% # drop_na() ## ----summarize-metrics, eval=FALSE-------------------------------------------- # # Get summary # MHW_summary <- OISST_clim %>% # dplyr::select(-c(doy, thresh, threshCriterion, durationCriterion, event)) %>% # mutate(season_num = month(as.Date(floor_date(t, unit = "season"))), # year = year(t)) %>% # mutate(season = recode_factor(season_num, `12` = "Summer", `3` = "Autumn", # `6` = "Winter", `9` = "Spring")) %>% # mutate(intensity = temp-seas) %>% # group_by(ECOREGION, PROVINCE, lon, lat, year, season) %>% # summarise(nevents = length(unique(event_no)), # nMHWdays = length(t), # int_cumulative = sum(intensity), # int_mean = mean(intensity), # int_max = max(intensity), .groups = "drop") %>% # pivot_longer(-c(ECOREGION, PROVINCE, lon, lat, year, season), # names_to = 'metrics', values_to = 'values') # # # Save for future use # # NB: 10.8 MB # saveRDS(MHW_summary,"~/Desktop/MHW_summary_NZ_1982_2021_OISST.Rds") ## ----nz-trend----------------------------------------------------------------- # # Load data from previous code chunk # MHW_summary_NZ <- readRDS("~/Desktop/MHW_summary_NZ_1982_2021_OISST.Rds") # # # Get count of unique pixels # MHW_summary_NZ_npix <- MHW_summary_NZ %>% # dplyr::select(lon, lat) %>% # dplyr::distinct() # # # Get annual summaries # MHW_annual_summary_NZ <- MHW_summary_NZ %>% # pivot_wider(names_from = metrics, values_from = values) %>% # group_by(year) %>% # summarize(Number_MHW_days = sum(nMHWdays)/nrow(MHW_summary_NZ_npix), # Nevents = sum(nevents)/nrow(MHW_summary_NZ_npix), # Mean_Intensity = mean(int_mean), # Maximum_Intensity = mean(int_max), # Cumulative_Intensity = sum(int_cumulative)/nrow(MHW_summary_NZ_npix), .groups = "drop") %>% # complete(year = 1982:2021) %>% # pivot_longer(-year, names_to = 'Metrics', values_to = 'values') %>% # mutate(Metrics = factor(Metrics, levels = c("Number_MHW_days", "Nevents", "Mean_Intensity", # "Maximum_Intensity", "Cumulative_Intensity"))) # # # Prepare prettier labels # metrics_labs <- c(`Number_MHW_days` = "Number of MHW days", # `Nevents` = "Number of Events", # `Mean_Intensity` = "Mean Intensity (°C)", # `Maximum_Intensity` = "Maximum Intensity (°C)", # `Cumulative_Intensity` = "Cumulative Intensity (°C Days)") # # # Create figure # p <- ggplot(MHW_annual_summary_NZ, aes(x = year, y = values, col = Metrics)) + # geom_line(size = 0.5) + # geom_smooth(se = T) + # scale_colour_viridis(begin = 0, end = 0.75, option = "viridis", discrete = T) + # facet_wrap(Metrics ~ ., scales = 'free', labeller = as_labeller(metrics_labs), ncol = 2) + # labs(x = "Year", Y = NULL) + # theme_bw() + # theme(legend.position = "none") # p # # # Coerce to interactive plotly format # # plotly::ggplotly(p) ## ----nz-trend2---------------------------------------------------------------- # # Get annual trends # MHW_annual_trends_NZ <- MHW_annual_summary_NZ %>% # group_by(Metrics) %>% # nest() %>% # mutate(ts_out = purrr::map(data, ~ts(.x$values, start = 1982, end = 2021, frequency = 1))) %>% # mutate(sens = purrr::map(ts_out, ~sens.slope(.x, conf.level = 0.95)), # pettitt = purrr::map(ts_out, ~pettitt.test(.x)), # lm = purrr::map(data, ~lm(values ~ year, .x))) %>% # mutate(Sens_Slope = as.numeric(unlist(sens)[1]), # P_Value = as.numeric(unlist(sens)[3]), # Change_Point_Year = time(ts_out[[1]])[as.numeric(unlist(pettitt)[3])], # Change_Point_pvalue = as.numeric(unlist(pettitt)[4]), # lm_slope = unlist(lm)$coefficients.year) %>% # # Add step of cutting time series in 2 using Change_Point_Year # mutate(pre_ts = purrr::map(ts_out, ~window(.x, start = 1982, end = Change_Point_Year)), # post_ts = purrr::map(ts_out, ~window(.x, start = Change_Point_Year, end = 2021))) %>% # # Add step of calculating sen's slope and p-value to pre and post change point year # mutate(sens_pre = purrr::map(pre_ts, ~sens.slope(.x, conf.level = 0.95)), # Sens_Slope_pre = as.numeric(unlist(sens_pre)[1]), P_Value_pre = as.numeric(unlist(sens_pre)[3]), # sens_post = purrr::map(post_ts, ~sens.slope(.x, conf.level = 0.95)), # Sens_Slope_post = as.numeric(unlist(sens_post)[1]), # P_Value_post = as.numeric(unlist(sens_post)[3])) %>% # dplyr::select(Metrics, Sens_Slope, P_Value, Change_Point_Year, Change_Point_pvalue, lm_slope, # Sens_Slope_pre, P_Value_pre, Sens_Slope_post, P_Value_post) # # # Create table # kable(MHW_annual_trends_NZ, caption = "Table 1 - Sens's Slope and p-value.") %>% # kable_classic() ## ----ecoregions-trends-------------------------------------------------------- # # Count pixels per realm # MHW_summary_NZ_npix_realms <- MHW_summary_NZ %>% # group_by(lon, lat, ECOREGION) %>% # tally() %>% # group_by(ECOREGION) %>% # tally() %>% # rename(npix = n) # MHW_summary_NZ_npix_realms # # # Summarise by realm # MHW_summary_NZ_realms <- left_join(MHW_summary_NZ, MHW_summary_NZ_npix_realms, by = 'ECOREGION') %>% # pivot_wider(names_from = metrics,values_from = values) %>% # group_by(year, ECOREGION, season, npix) %>% # summarize(Number_MHW_days = sum(nMHWdays), # Nevents = sum(nevents), # Mean_Intensity = mean(int_mean), # Maximum_Intensity = mean(int_max), # Cumulative_Intensity = sum(int_cumulative), .groups = "drop") %>% # mutate(Number_MHW_days = Number_MHW_days/npix, # Nevents = Nevents/npix, # Cumulative_Intensity = Cumulative_Intensity/npix) %>% # group_by(ECOREGION, season) %>% # complete(year = 1982:2021) %>% # dplyr::select(-npix) %>% # pivot_longer(-c(year, ECOREGION, season), names_to = 'Metrics', values_to = 'values') %>% # mutate(Metrics = factor(Metrics, levels = c("Number_MHW_days", "Nevents", "Mean_Intensity", # "Maximum_Intensity", "Cumulative_Intensity"))) %>% # replace(is.na(.), 0) # # # Plot the results # p <- MHW_summary_NZ_realms %>% # dplyr::filter(Metrics == 'Number_MHW_days') %>% # ggplot(aes(x = year, y = values, colour = season)) + # geom_line(size = 0.5) + # geom_smooth(se = T) + # labs(x = "Year", y = "Number of MHW days") + # scale_colour_viridis(begin = 0, end = 0.75, option = "inferno", discrete = T) + # guides(colour = guide_legend(override.aes = list(shape = 15, size = 2))) + # facet_wrap(ECOREGION~., ncol = 3) + # theme_bw() + # theme(legend.position = c(0.8,0.05), # legend.title = element_blank()) # p # # # Coerce to plotly format # # plotly::ggplotly(p) ## ----ecoregions-trends2------------------------------------------------------- # # Get trends # MHW_trends_NZ_realms <- MHW_summary_NZ_realms %>% # group_by(Metrics, ECOREGION, season) %>% # nest() %>% # mutate(ts_out = purrr::map(data, ~ts(.x$values, start = 1982, end = 2021, frequency = 1))) %>% # mutate(sens = purrr::map(ts_out, ~sens.slope(.x, conf.level = 0.95)), # pettitt = purrr::map(ts_out, ~pettitt.test(.x)), # lm = purrr::map(data, ~lm(values ~ year,.x))) %>% # mutate(Sens_Slope = as.numeric(unlist(sens)[1]), # P_Value = as.numeric(unlist(sens)[3]), # Change_Point_Year = time(ts_out[[1]])[as.numeric(unlist(pettitt)[3])], # Change_Point_pvalue = as.numeric(unlist(pettitt)[4]), # lm_slope = unlist(lm)$coefficients.year) %>% # # Add step of cutting time series in 2 using Change_Point_Year # mutate(pre_ts = purrr::map(ts_out, ~window(.x, start = 1982, end = Change_Point_Year)), # post_ts = purrr::map(ts_out, ~window(.x, start = Change_Point_Year, end = 2021))) %>% # # Add step of calculating sen's slope and p-value to pre and post change point year # mutate(sens_pre = purrr::map(pre_ts, ~sens.slope(.x, conf.level = 0.95)), # Sens_Slope_pre = as.numeric(unlist(sens_pre)[1]), # P_Value_pre = as.numeric(unlist(sens_pre)[3]), # sens_post = purrr::map(post_ts, ~sens.slope(.x, conf.level = 0.95)), # Sens_Slope_post = as.numeric(unlist(sens_post)[1]), # P_Value_post = as.numeric(unlist(sens_post)[3])) %>% # ungroup() %>% # To remove Season column # dplyr::select(ECOREGION, Metrics, season, Sens_Slope, P_Value, Change_Point_Year, Change_Point_pvalue, # lm_slope, Sens_Slope_pre, P_Value_pre, Sens_Slope_post, P_Value_post) %>% # dplyr::filter(Metrics == 'Number_MHW_days') # # # Create table # kable(MHW_trends_NZ_realms, caption = "Table 2 - Sens's Slope and p-value.") %>% # kable_classic() ## ----pixel, eval=FALSE-------------------------------------------------------- # # Simple event detection # MHW_WA <- detect_event(ts2clm(sst_WA, climatologyPeriod = c("1982-01-01", "2011-12-31"))) # # # Get annual metrics # MHW_WA_metrics_year <- MHW_WA$climatology %>% # dplyr::filter(!is.na(event_no)) %>% # mutate(year = year(t), # intensity = temp-seas) %>% # group_by(event, year) %>% # summarise(nevents = length(unique(event_no)), # nMHWdays = length(t), # int_cumulative = sum(intensity), # int_mean = mean(intensity), # int_max = max(intensity), .groups = "drop") %>% # # Need to complete time series in case of years with no MHWs in order to get full ts in trend analysis # complete(year = 1982:2018) %>% # dplyr::select(-event) %>% # pivot_longer(-c(year), names_to = 'Metrics', values_to = 'values') %>% # replace(is.na(.), 0) # # # Plot the data # p <- ggplot(MHW_WA_metrics_year, aes(x = year, y = values)) + # geom_point() + # geom_line() + # geom_smooth(se = T, method = 'lm') + # labs(x = "Year") + # facet_wrap(~Metrics, scales = 'free') + # theme_bw() + # theme(legend.position = "none") # p # # # Coerce to plotly # # plotly::ggplotly(p) ## ----eval=FALSE--------------------------------------------------------------- # MHW_WA_trends <- MHW_WA_metrics_year %>% # group_by(Metrics) %>% # nest() %>% # mutate(ts_out = purrr::map(data, ~ts(.x$values, start = 1982, end = 2018, frequency = 1))) %>% # mutate(sens = purrr::map(ts_out, ~sens.slope(.x, conf.level = 0.95)), # pettitt = purrr::map(ts_out, ~pettitt.test(.x)), # lm = purrr::map(data,~lm(values ~ year, .x))) %>% # mutate(Sens_Slope = as.numeric(unlist(sens)[1]), # P_Value = as.numeric(unlist(sens)[3]), # Change_Point_Year = time(ts_out[[1]])[as.numeric(unlist(pettitt)[3])], # Change_Point_pvalue = as.numeric(unlist(pettitt)[4]), # lm_slope = unlist(lm)$coefficients.year) %>% # # Add step of cutting time series in 2 using Change_Point_Year # mutate(pre_ts = purrr::map(ts_out,~window(.x, start = 1982, end = Change_Point_Year)), # post_ts = purrr::map(ts_out,~window(.x, start = Change_Point_Year, end = 2018))) %>% # # Add step of calculating sen's slope and p-value to pre and post change point year # mutate(sens_pre = purrr::map(pre_ts, ~sens.slope(.x, conf.level = 0.95)), # Sens_Slope_pre = as.numeric(unlist(sens_pre)[1]), # P_Value_pre = as.numeric(unlist(sens_pre)[3]), # sens_post = purrr::map(post_ts, ~sens.slope(.x, conf.level = 0.95)), # Sens_Slope_post = as.numeric(unlist(sens_post)[1]), # P_Value_post = as.numeric(unlist(sens_post)[3])) %>% # dplyr::select(Metrics, Sens_Slope, P_Value, Change_Point_Year, Change_Point_pvalue, # lm_slope, Sens_Slope_pre, P_Value_pre, Sens_Slope_post, P_Value_post) # # # Create table # kable(MHW_WA_trends, caption = "Table 3 - Sens's Slope and p-value.") %>% # kable_classic()