AIM : Show the utility and limitation of running a trend analysis in MHW metrics in conjonction of the heatwaveR package.

Document Note : The maps, plots and app within this document are interactive so make sure you give them a play like zooming in and out in the maps but also on the plots. Clicking on the legend allows to only select and display the time series needed.

Table of contents

  1. Regional case scenario

  2. Pixel-based case scenario

  3. Bibliography

Regional case scenario

Data - Downloading and Preparing NOAA OISST Data By Robert W Schlegel and AJ Smit From https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html

Event Detection - Uses heatwaveR:: (Schlegel and Smit (2018)), translation of GitHub python functions from Eric C.J. Oliver (published in paper Holbrook et al. (2019) A global assessment of marine heatwaves and their drivers) https://robwschlegel.github.io/heatwaveR/articles/gridded_event_detection.html

The file YearSeason_NZ_BIOREGION_MeanMetrics_MHWDays_Pixels_1982_2021_OISST.Rds results from the functions from Schlegel and Smit (2018), used on OISST v2.1 data for NZ EEZ, using climatologyPeriod = c(“1982-01-01,” “2011-01-01”). MHW metrics are grouped by season and year for each pixels following the methology of Thoral et al. (2022).

Here we load the resulting dataframe for sake of clarity.

Year of most intense MHW and MEOW classification 

mhw_metrics_mean_realm <- readRDS("C:/Users/thoralf/OneDrive - NIWA/Documents/Collab/NZ_Coastal_MHW_2022/Data/YearSeason_NZ_BIOREGION_MeanMetrics_MHWDays_Pixels_1982_2021_OISST.Rds")

most_extreme_MHW_year_intensity <- mhw_metrics_mean_realm %>% 
  dplyr::filter(metrics == 'int_cumulative') %>% 
  group_by(lon, lat, year) %>%
  summarise(max_intensity = max(values)) %>% 
  dplyr::filter(max_intensity==max(max_intensity))
## `summarise()` has grouped output by 'lon', 'lat'. You can override using the `.groups` argument.
pal <- c(brewer.pal(9,'Purples'),brewer.pal(9,'Blues'),brewer.pal(9,'Oranges'),brewer.pal(9,'Reds'))
p_intensity <- ggplot() + 
  #geom_tile(data=most_extreme_MHW_year_intensity,aes(x=lon,y=lat,fill=year)) + 
  geom_raster(data=most_extreme_MHW_year_intensity,aes(x=lon,y=lat,fill=year)) + 
  scale_fill_gradientn(colours=pal) + 
  ggtitle('Year of Most Intense MHW events (1981-2021)') + 
  borders('world2',xlim=c(170,200),ylim=c(-60,-20),fill='grey') + 
  #xlim(c(170,180)) + ylim(c(-40,-30)) +
  theme_bw()
ggplotly(p_intensity)

Figure 1 - Year of Most Intense MHW around New Zealand EEZ.

MEOW from Spalding et al. (2007).

meow_NZ_lonlat2 <- st_read('C:/Users/thoralf/OneDrive - NIWA/Documents/InSitu_Data/Marine_Ecoregions_Of_the_World_(MEOW)-shp/Marine_Ecoregions_Of_the_World__MEOW_.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")))
## Reading layer `Marine_Ecoregions_Of_the_World__MEOW_' from data source 
##   `C:\Users\thoralf\OneDrive - NIWA\Documents\InSitu_Data\Marine_Ecoregions_Of_the_World_(MEOW)-shp\Marine_Ecoregions_Of_the_World__MEOW_.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 232 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -20037510 ymin: -30240970 xmax: 20037510 ymax: 23063390
## Projected CRS: WGS 84 / Pseudo-Mercator
map_bioregion <- ggplot() +
  geom_sf(data=meow_NZ_lonlat2,aes(fill=ECOREGION)) +
  scale_fill_viridis(begin = 0, end = .75,option="viridis",discrete = T) +
  borders('world2',xlim=c(160,190),ylim=c(-55,-25),fill = 'grey') + xlim(c(160,190)) + ylim(c(-55,-25)) +
  theme(panel.background = element_rect(fill = 'aliceblue')) +
  xlab('') + ylab('') +
  theme(legend.position = "right",
        legend.title = element_blank(),
        legend.text = element_text(size=16))
map_bioregion
Figure 2 - MEOW ecoregions around New Zealand EEZ.

Figure 2 - MEOW ecoregions around New Zealand EEZ.

Pixel-based case scenario

We use the time series of SST available within the heatwaveR package (sst_WA) and investigate the trend analysis in 5 metrics by grouping metric values per year as previously.

ts <- ts2clm(sst_WA, climatologyPeriod = c("1982-01-01", "2011-12-31"))
mhw <- detect_event(ts)

mhw_metrics_year <- mhw$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)) %>% 
  # 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) %>% 
  ungroup() %>% 
  dplyr::select(-event) %>% 
  pivot_longer(-c(year),names_to='Metrics',values_to='values')
## `summarise()` has grouped output by 'event'. You can override using the `.groups` argument.
#Replace NAs by 0?
mhw_metrics_year$values[is.na(mhw_metrics_year$values)] <- 0

p <- ggplot(mhw_metrics_year,aes(x=year,y=values)) +
  geom_point() +
  geom_line() +
  facet_wrap(~Metrics,scales = 'free') +
  geom_smooth(se=T,method = 'lm') +
  xlab('Year') + 
  theme_bw() + 
  theme(legend.position="none") 
  
ggplotly(p)
## `geom_smooth()` using formula 'y ~ x'

Figure 5. Trend in 5 MHW metrics for the sst_WA time series.

The lm() regression seems to suject some upward trends in the metrics. What does Mann-Kendall have to say about these?

MHW_trends <- mhw_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) %>% #,
         #lm_pvalue = parameters::p_value(aa$lm[[1]])[2,2]) %>% 
  #How to include p-value of lm? summary(aa$lm[[3]])$coefficient[,4] or parameters::p_value(aa$lm[[3]])[2,2]
  # 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)


kable(MHW_trends,caption = "Table 3 - Sens's Slope and p-value.")  %>%
  kable_classic()
Table 3 - Sens’s Slope and p-value.
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
nevents 0.0000000 0.4914680 2006 0.4124316 0.0346136 0.0000000 0.3784001 -0.2916667 0.1710692
nMHWdays 0.0000000 0.3216169 1995 0.3395085 0.9663348 -0.5555556 0.1416319 0.0000000 0.8998992
int_cumulative 0.0000000 0.3152513 1995 0.3204836 2.0934086 -0.8027143 0.1124402 0.0000000 0.9198687
int_mean 0.0000000 0.3152513 1995 0.2847829 0.0097617 -0.0084139 0.1763236 0.0000000 0.7628066
int_max 0.0064262 0.1939960 1995 0.2223932 0.0292780 -0.0277000 0.1416319 0.0000000 0.9598837

For some reason, the Mann-Kendall and Sen Slope analyses don’t seem to be useful here as they return suspiciously too low trends (vs lm_slope, which is the traditional linear regression using lm()) and quite high p-values

Summarizing MHW metrics by year in a given pixel will likely result in some years with no MHWs, hence bringing some 0 values into the time series. I am not too sure how sensitive the MK and Pettitt (breakpoint) analyses are to 0 values, but I can dig this out.

Bibliography

Holbrook, Neil J, Hillary A Scannell, Alexander Sen Gupta, Jessica A Benthuysen, Ming Feng, Eric CJ Oliver, Lisa V Alexander, et al. 2019. “A Global Assessment of Marine Heatwaves and Their Drivers.” Nature Communications 10 (1): 1–13.
Schlegel, Robert W., and Albertus J. Smit. 2018. “heatwaveR: A Central Algorithm for the Detection of Heatwaves and Cold-Spells.” Journal of Open Source Software 3 (27): 821. https://doi.org/10.21105/joss.00821.
Spalding, Mark D, Helen E Fox, Gerald R Allen, Nick Davidson, Zach A Ferdaña, MAX Finlayson, Benjamin S Halpern, et al. 2007. “Marine Ecoregions of the World: A Bioregionalization of Coastal and Shelf Areas.” BioScience 57 (7): 573–83.
Thoral, François, Shinae Montie, Mads S Thomsen, Leigh W Tait, Matthew H Pinkerton, and David R Schiel. 2022. “Unravelling Seasonal Trends in Coastal Marine Heatwave Metrics Across Global Biogeographical Realms.” Scientific Reports 12 (1): 1–13.