WAP-Heatwaves.Rmd

---
title: Detection of extreme heatwaves and low sea ice cover in the Antarctic Peninsula  for
  evaluation of effects over Krill
author: "Lucas Krüger, Maurício Mardones, Lorena Rebolledo"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#---------- Download NOAA erddap optimum interpolation SST----

#this part of the script is a modification of the scripts in the heatwaveR vignette
#available at https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html


```{r, echo =FALSE, message=FALSE}
#---------- Download NOAA erddap optimum interpolation SST----

# The packages we will use
library(dplyr) # A staple for modern data management in R
library(lubridate) # Useful functions for dealing with dates
library(ggplot2) # The preferred library for data visualisation
library(patchwork)
library(tidync) # For easily dealing with NetCDF data
library(rerddap) # For easily downloading subsets of data
library(doParallel) # For parallel processing
library(raster)
library(terra)
library(sf)
library(ncdf4)
library(mapdata)
library(mapproj)
library(fields)
library(sp)
#library(rgdal)

```

# The information for the NOAA OISST data

```{r, echo=FALSE}


rerddap::info(datasetid = "ncdcOisst21Agg_LonPM180", url = "https://coastwatch.pfeg.noaa.gov/erddap/")

```

# Download and prepare data based on user provided start and end dates using three different requests per geographical strata so it can be run on slow PCs.

```{r}
# spatial units to use


shp<-shapefile("Clipped_StrataGCS.shp")


crs=CRS( "+proj=longlat +datum=WGS84 +no_defs") # coordinate reference system CRS

shp2<-spTransform(shp,crs) # be sure it is in the correct coordinate system


#load antarctic peninsula shapefile to plot it
# this shapefile was cut from high resolution file from the SCAR Antarctic Digital Database, 2023
# Gerrish, L., Ireland, L., Fretwell, P., & Cooper, P. (2023). High resolution vector polygons of the Antarctic coastline (7.7) [Data set]. UK Polar Data Centre, Natural Environment Research Council, UK Research & Innovation. https://doi.org/10.5285/0be5339c-9d35-44c9-a10f-da4b5356840b'

wap<-shapefile("peninsula_correct.shp") 


# plot map with shapefile (Fig 1)
shpdf<-data.frame(shp2)




shp2$ID[shp2$ID=="JOIN"]<-"JI"
shp2$ID[shp2$ID=="Extra"]<-"GS"
shp2$ID[shp2$ID=="SSIW"]<-"SS"


shp2$ID[shp2$ID=="OuterSW"]<-"SW"

shp2<-subset(shp2,ID!="OuterNW"|ID!="OuterSE"|ID!="OuterNE")

ggplot()+
  
  xlab("Longitude")+ylab("Latitude")+
  theme_bw()+
    geom_polygon(data=shp2,aes(x = long, y = lat, group = group,label=group),fill="lightblue",colour="black")+
  geom_polygon(data=wap,aes(x = long, y = lat, group = group))+
geom_text(data=shpdf,aes(label=ID, 
            x=Labx,
            y=Laby),colour="red2",
            size = 3, hjust=0, vjust=-1) +
      coord_quickmap(xlim=c(-75,-50),ylim=c(-75, -59), expand = TRUE, clip = "on")


shp.sf<-st_as_sf(shp)
unique(shp.sf$ID)

#Extract coordinates and generates a bounding box 

geom<-st_geometry(shp.sf)

bbox_list <- lapply(st_geometry(shp.sf), st_bbox)
```


```{r}
#min and max coordinates of each strata as data frame

df1<-data.frame(t(as.data.frame((bbox_list[[1]])))) # Elephant
df2<-data.frame(t(as.data.frame((bbox_list[[2]])))) # South Shetlands
df3<-data.frame(t(as.data.frame((bbox_list[[3]])))) # Bransfield
df4<-data.frame(t(as.data.frame((bbox_list[[4]])))) #Joinville
df5<-data.frame(t(as.data.frame((bbox_list[[5]])))) # Outer (gerlache?)
df9<-data.frame(t(as.data.frame((bbox_list[[9]])))) #Extra SW

```

#download and process data for Elephant Island Strata

```{r, echo=FALSE}


OISST_sub_dl <- function(time_df){
  OISST_DATA1 <- rerddap::griddap(datasetx = "ncdcOisst21Agg_LonPM180",
                                url = "https://coastwatch.pfeg.noaa.gov/erddap/", 
                                time = c(time_df$start, time_df$end), 
                                zlev = c(0, 0),
                                latitude = c(min(df1$ymin),max(df1$ymax)), 
                                longitude = c(min(df1$xmin),max(df1$xmax)),
                                fields = c("sst","ice"))$data %>% 
    dplyr::mutate(time = base::as.Date(stringr::str_remove(time, "T12:00:00Z"))) %>% 
    dplyr::rename(t = time, temp = sst, sic=ice,lon = longitude, lat = latitude) %>% 
    dplyr::select(lon, lat, t, temp,sic) %>% 
    stats::na.omit()
}




dl_years1 <- data.frame(date_index = 1:3,
                        start = as.character(c("1985-01-01","1990-01-01","1995-01-01")),
                        
                        end = as.character(c("1989-12-31","1994-12-31","1999-12-31")))


dl_years2 <- data.frame(date_index = 1:3,
                        start = as.character(c("2000-01-01", "2005-01-01", "2010-01-01")),
                        
                        end = as.character(c("2004-12-31", "2009-12-31","2014-12-31")))


dl_years3 <- data.frame(date_index = 1:3,
                        start = as.character(c("2015-01-01","2020-01-01","2021-01-01")),
                        
                        end = as.character(c("2019-12-31","2020-12-31","2022-12-31")))


###download and save data


base::system.time(
  OISST_EL1 <- dl_years1 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
) 


base::system.time(
  OISST_EL2 <- dl_years2 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)


base::system.time(
  OISST_EL3 <- dl_years3 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)

base::saveRDS(OISST_EL1, file = "OISST_EL1.Rds")
base::saveRDS(OISST_EL2, file = "OISST_EL2.Rds")
base::saveRDS(OISST_EL3, file = "OISST_EL3.Rds")

### summarize data per day

EL1<-plyr::ddply(OISST_EL1, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

EL2<-plyr::ddply(OISST_EL2, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

EL3<-plyr::ddply(OISST_EL3, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

### merge data

EL12<-rbind(EL1,EL2)
EL123<-rbind(EL12,EL3)

base::saveRDS(EL123, file = "Elephant_Strata.Rds")

#remove data and release unused memory

rm(list=c("OISST_EL1","OISST_EL2","OISST_EL3","OISST_sub_dl","shp","bbox_list",
          "EL1","EL2","EL3","EL12","EL123"))

gc()



```
#download and process data for South Shetland West Strata

```{r, echo=FALSE}


OISST_sub_dl <- function(time_df){
  OISST_DATA1 <- rerddap::griddap(datasetx = "ncdcOisst21Agg_LonPM180",
                                url = "https://coastwatch.pfeg.noaa.gov/erddap/", 
                                time = c(time_df$start, time_df$end), 
                                zlev = c(0, 0),
                                latitude = c(min(df2$ymin),max(df2$ymax)),
                                longitude = c(min(df2$xmin),max(df2$xmax)),
                                fields = c("sst","ice"))$data %>% 
    dplyr::mutate(time = base::as.Date(stringr::str_remove(time, "T12:00:00Z"))) %>% 
    dplyr::rename(t = time, temp = sst, sic=ice,lon = longitude, lat = latitude) %>% 
    dplyr::select(lon, lat, t, temp,sic) %>% 
    stats::na.omit()
}




dl_years1 <- data.frame(date_index = 1:3,
                        start = as.character(c("1985-01-01","1990-01-01","1995-01-01")),
                        
                        end = as.character(c("1989-12-31","1994-12-31","1999-12-31")))


dl_years2 <- data.frame(date_index = 1:3,
                        start = as.character(c("2000-01-01", "2005-01-01", "2010-01-01")),
                        
                        end = as.character(c("2004-12-31", "2009-12-31","2014-12-31")))


dl_years3 <- data.frame(date_index = 1:3,
                        start = as.character(c("2015-01-01","2020-01-01","2021-01-01")),
                        
                        end = as.character(c("2019-12-31","2020-12-31","2022-12-31")))


###download and save data


base::system.time(
  OISST_SS1 <- dl_years1 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
) 


base::system.time(
  OISST_SS2 <- dl_years2 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)


base::system.time(
  OISST_SS3 <- dl_years3 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)

base::saveRDS(OISST_SS1, file = "OISST_SS1.Rds")
base::saveRDS(OISST_SS2, file = "OISST_SS2.Rds")
base::saveRDS(OISST_SS3, file = "OISST_SS3.Rds")

### summarize data per day

SS1<-plyr::ddply(OISST_SS1, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

SS2<-plyr::ddply(OISST_SS2, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

SS3<-plyr::ddply(OISST_SS3, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

### merge data

SS12<-rbind(SS1,SS2)
SS123<-rbind(SS12,SS3)

base::saveRDS(SS123, file = "West_SSI_Strata.Rds")

#remove data and release unused memory

rm(list=c("OISST_SS1","OISST_SS2","OISST_SS3","OISST_sub_dl",
          "SS1","SS2","SS3","SS12","SS123"))

gc()



```


#download and process data for Bransfield Strait Strata

```{r, echo=FALSE}

df3$ID[3]

OISST_sub_dl <- function(time_df){
  OISST_DATA1 <- rerddap::griddap(datasetx = "ncdcOisst21Agg_LonPM180",
                                url = "https://coastwatch.pfeg.noaa.gov/erddap/", 
                                time = c(time_df$start, time_df$end), 
                                zlev = c(0, 0),
                                latitude = c(min(df3$ymin),max(df3$ymax)),
                                longitude = c(min(df3$xmin),max(df3$xmax)),
                                fields = c("sst","ice"))$data %>% 
    dplyr::mutate(time = base::as.Date(stringr::str_remove(time, "T12:00:00Z"))) %>% 
    dplyr::rename(t = time, temp = sst, sic=ice,lon = longitude, lat = latitude) %>% 
    dplyr::select(lon, lat, t, temp,sic) %>% 
    stats::na.omit()
}




dl_years1 <- data.frame(date_index = 1:3,
                        start = as.character(c("1985-01-01","1990-01-01","1995-01-01")),
                        
                        end = as.character(c("1989-12-31","1994-12-31","1999-12-31")))


dl_years2 <- data.frame(date_index = 1:3,
                        start = as.character(c("2000-01-01", "2005-01-01", "2010-01-01")),
                        
                        end = as.character(c("2004-12-31", "2009-12-31","2014-12-31")))


dl_years3 <- data.frame(date_index = 1:3,
                        start = as.character(c("2015-01-01","2020-01-01","2021-01-01")),
                        
                        end = as.character(c("2019-12-31","2020-12-31","2022-12-31")))


###download and save data


base::system.time(
  OISST_BS1 <- dl_years1 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
) 


base::system.time(
  OISST_BS2 <- dl_years2 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)


base::system.time(
  OISST_BS3 <- dl_years3 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)

base::saveRDS(OISST_BS1, file = "OISST_BS1.Rds")
base::saveRDS(OISST_BS2, file = "OISST_BS2.Rds")
base::saveRDS(OISST_BS3, file = "OISST_BS3.Rds")

OISST_BS1<-readRDS("OISST_BS1.Rds")
OISST_BS2<-readRDS("OISST_BS2.Rds")
OISST_BS3<-readRDS("OISST_BS3.Rds")
### summarize data per day

BS1<-plyr::ddply(OISST_BS1, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

BS2<-plyr::ddply(OISST_BS2, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

BS3<-plyr::ddply(OISST_BS3, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

### merge data



BS12<-rbind(BS1,BS2)
BS123<-rbind(BS12,BS3)

base::saveRDS(BS123, file = "Bransfield_Strata.Rds")

#remove data and release unused memory

rm(list=c("OISST_BS1","OISST_BS2","OISST_BS3","OISST_sub_dl",
          "BS1","BS2","BS3","BS12","BS123"))

gc()



```




#download and process data for Joinville Island Strata

```{r, echo=FALSE}



OISST_sub_dl <- function(time_df){
  OISST_DATA1 <- rerddap::griddap(datasetx = "ncdcOisst21Agg_LonPM180",
                                url = "https://coastwatch.pfeg.noaa.gov/erddap/", 
                                time = c(time_df$start, time_df$end), 
                                zlev = c(0, 0),
                                latitude = c(min(df4$ymin),max(df4$ymax)),
                                longitude = c(min(df4$xmin),max(df4$xmax)),
                                fields = c("sst","ice"))$data %>% 
    dplyr::mutate(time = base::as.Date(stringr::str_remove(time, "T12:00:00Z"))) %>% 
    dplyr::rename(t = time, temp = sst, sic=ice,lon = longitude, lat = latitude) %>% 
    dplyr::select(lon, lat, t, temp,sic) %>% 
    stats::na.omit()
}




dl_years1 <- data.frame(date_index = 1:3,
                        start = as.character(c("1985-01-01","1990-01-01","1995-01-01")),
                        
                        end = as.character(c("1989-12-31","1994-12-31","1999-12-31")))


dl_years2 <- data.frame(date_index = 1:3,
                        start = as.character(c("2000-01-01", "2005-01-01", "2010-01-01")),
                        
                        end = as.character(c("2004-12-31", "2009-12-31","2014-12-31")))


dl_years3 <- data.frame(date_index = 1:3,
                        start = as.character(c("2015-01-01","2020-01-01","2021-01-01")),
                        
                        end = as.character(c("2019-12-31","2020-12-31","2022-12-31")))


###download and save data


base::system.time(
  OISST_JI1 <- dl_years1 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
) 


base::system.time(
  OISST_JI2 <- dl_years2 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)


base::system.time(
  OISST_JI3 <- dl_years3 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)

base::saveRDS(OISST_JI1, file = "OISST_JI1.Rds")
base::saveRDS(OISST_JI2, file = "OISST_JI2.Rds")
base::saveRDS(OISST_JI3, file = "OISST_JI3.Rds")


#OISST_JI1<-base::readRDS("OISST_JI1.Rds")
#OISST_JI2<-base::readRDS("OISST_JI2.Rds")
#OISST_JI3<-base::readRDS("OISST_JI3.Rds")


### summarize data per day

JI1<-plyr::ddply(OISST_JI1, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

JI2<-plyr::ddply(OISST_JI2, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

JI3<-plyr::ddply(OISST_JI3, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

### merge data

JI12<-rbind(JI1,JI2)
JI123<-rbind(JI12,JI3)

base::saveRDS(JI123, file = "Joinvile_Strata.Rds")

#remove data and release unused memory

rm(list=c("OISST_JI1","OISST_JI2","OISST_JI3","OISST_sub_dl",
          "JI1","JI2","JI3","JI12","JI123"))

gc()



```


#download and process data for Extra Strata (Gerlache)

```{r, echo=FALSE}


OISST_sub_dl <- function(time_df){
  OISST_DATA1 <- rerddap::griddap(datasetx = "ncdcOisst21Agg_LonPM180",
                                url = "https://coastwatch.pfeg.noaa.gov/erddap/", 
                                time = c(time_df$start, time_df$end), 
                                zlev = c(0, 0),
                                latitude = c(min(df5$ymin),max(df5$ymax)),
                                longitude = c(min(df5$xmin),max(df5$xmax)),
                                fields = c("sst","ice"))$data %>% 
    dplyr::mutate(time = base::as.Date(stringr::str_remove(time, "T12:00:00Z"))) %>% 
    dplyr::rename(t = time, temp = sst, sic=ice,lon = longitude, lat = latitude) %>% 
    dplyr::select(lon, lat, t, temp,sic) %>% 
    stats::na.omit()
}




dl_years1 <- data.frame(date_index = 1:3,
                        start = as.character(c("1985-01-01","1990-01-01","1995-01-01")),
                        
                        end = as.character(c("1989-12-31","1994-12-31","1999-12-31")))


dl_years2 <- data.frame(date_index = 1:3,
                        start = as.character(c("2000-01-01", "2005-01-01", "2010-01-01")),
                        
                        end = as.character(c("2004-12-31", "2009-12-31","2014-12-31")))


dl_years3 <- data.frame(date_index = 1:3,
                        start = as.character(c("2015-01-01","2020-01-01","2021-01-01")),
                        
                        end = as.character(c("2019-12-31","2020-12-31","2022-12-31")))


###download and save data


base::system.time(
  OISST_EX1 <- dl_years1 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
) 


base::system.time(
  OISST_EX2 <- dl_years2 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)


base::system.time(
  OISST_EX3 <- dl_years3 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)

base::saveRDS(OISST_EX1, file = "OISST_EX1.Rds")
base::saveRDS(OISST_EX2, file = "OISST_EX2.Rds")
base::saveRDS(OISST_EX3, file = "OISST_EX3.Rds")

### summarize data per day

EX1<-plyr::ddply(OISST_EX1, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

EX2<-plyr::ddply(OISST_EX2, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

EX3<-plyr::ddply(OISST_EX3, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

### merge data

EX12<-rbind(EX1,EX2)
EX123<-rbind(EX12,EX3)

base::saveRDS(EX123, file = "Extra_Strata.Rds")

#remove data and release unused memory

rm(list=c("OISST_EX1","OISST_EX2","OISST_EX3","OISST_sub_dl",
          "EX1","EX2","EX3","EX12","EX123"))

gc()



```


#download and process data for Outer SW Strata (Grandidier Channel, Margueritte Bay)

```{r, echo=FALSE}


OISST_sub_dl <- function(time_df){
  OISST_DATA1 <- rerddap::griddap(datasetx = "ncdcOisst21Agg_LonPM180",
                                url = "https://coastwatch.pfeg.noaa.gov/erddap/", 
                                time = c(time_df$start, time_df$end), 
                                zlev = c(0, 0),
                                latitude = c(min(df9$ymin),max(df9$ymax)),
                                longitude = c(min(df9$xmin),max(df9$xmax)),
                                fields = c("sst","ice"))$data %>% 
    dplyr::mutate(time = base::as.Date(stringr::str_remove(time, "T12:00:00Z"))) %>% 
    dplyr::rename(t = time, temp = sst, sic=ice,lon = longitude, lat = latitude) %>% 
    dplyr::select(lon, lat, t, temp,sic) %>% 
    stats::na.omit()
}




dl_years1 <- data.frame(date_index = 1:3,
                        start = as.character(c("1985-01-01","1990-01-01","1995-01-01")),
                        
                        end = as.character(c("1989-12-31","1994-12-31","1999-12-31")))


dl_years2 <- data.frame(date_index = 1:3,
                        start = as.character(c("2000-01-01", "2005-01-01", "2010-01-01")),
                        
                        end = as.character(c("2004-12-31", "2009-12-31","2014-12-31")))


dl_years3 <- data.frame(date_index = 1:3,
                        start = as.character(c("2015-01-01","2020-01-01","2021-01-01")),
                        
                        end = as.character(c("2019-12-31","2020-12-31","2022-12-31")))


###download and save data


base::system.time(
  OISST_SW1 <- dl_years1 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
) 


base::system.time(
  OISST_SW2 <- dl_years2 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)


base::system.time(
  OISST_SW3 <- dl_years3 %>% 
    dplyr::group_by(date_index) %>% 
    dplyr::group_modify(~OISST_sub_dl(.x)) %>% 
    dplyr::ungroup() %>% 
    dplyr::select(lon, lat, t, temp,sic)
)

base::saveRDS(OISST_SW1, file = "OISST_SW1.Rds")
base::saveRDS(OISST_SW2, file = "OISST_SW2.Rds")
base::saveRDS(OISST_SW3, file = "OISST_SW3.Rds")

### summarize data per day

SW1<-plyr::ddply(OISST_SW1, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

SW2<-plyr::ddply(OISST_SW2, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

SW3<-plyr::ddply(OISST_SW3, c("t"), summarise,
            tSD=sd(temp),tMax=max(temp),
            tMin=min(temp),temp=mean(temp),
            sicMean=mean(sic),sicMax=max(sic),
            sicMin=min(sic),sicSD=sd(sic))

### merge data

SW12<-rbind(SW1,SW2)
SW123<-rbind(SW12,SW3)

base::saveRDS(SW123, file = "OuterSW_Strata.Rds")

#remove data and release unused memory

rm(list=c("OISST_SW1","OISST_SW2","OISST_SW3","OISST_sub_dl",
          "SW1","SW2","SW3","SW12","SW123"))

gc()



```


###----------- calculating marine heatwaves-------------


```{r,echo=FALSE}

library(heatwaveR)


#load and merge data

BS<-base::readRDS("Bransfield_Strata.Rds")
EI<-base::readRDS("Elephant_Strata.Rds")
EX<-base::readRDS("Extra_Strata.Rds")
JI<-base::readRDS("Joinvile_Strata.Rds")
SS<-base::readRDS("West_SSI_Strata.Rds")
SW<-base::readRDS("OuterSW_Strata.Rds")

### calculate thresholds from climatology

# tMax thresholds
tMax_clim.BS <- ts2clm(data = BS, 
                     y = tMax, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 90)

tMax_clim.EI <- ts2clm(data = EI, 
                     y = tMax, climatologyPeriod = c("1986-01-01", "2022-12-31"), pctile = 90)

tMax_clim.EX <- ts2clm(data = EX, 
                     y = tMax, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 90)

tMax_clim.JI <- ts2clm(data = JI, 
                     y = tMax, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 90)

tMax_clim.SS <- ts2clm(data = SS, 
                     y = tMax, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 90)

tMax_clim.SW <- ts2clm(data = SW, 
                     y = tMax, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 90)



# iMin thresholds

sicMean_clim.BS <- ts2clm(data = BS, 
                     y = sicMean, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 20)
sicMean_clim.EI <- ts2clm(data = EI, 
                     y = sicMean, climatologyPeriod = c("1986-01-01", "2022-12-31"), pctile = 20)
sicMean_clim.EX <- ts2clm(data = EX, 
                     y = sicMean, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 20)

sicMean_clim.JI <- ts2clm(data = JI, 
                     y = sicMean, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 20)
sicMean_clim.SS <- ts2clm(data = SS, 
                     y = sicMean, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 20)
sicMean_clim.SW <- ts2clm(data = SW, 
                     y = sicMean, climatologyPeriod = c("1985-01-01", "2022-12-31"), pctile = 20)



#DETECT alldfS
eBS.t <- detect_event(data = tMax_clim.BS, y = tMax,minDuration = 5,maxGap = 7) 
eEI.t <- detect_event(data = tMax_clim.EI, y = tMax,minDuration = 5,maxGap = 7) 
eEX.t <- detect_event(data = tMax_clim.EX, y = tMax,minDuration = 5,maxGap = 7) 
eJI.t <- detect_event(data = tMax_clim.JI, y = tMax,minDuration = 5,maxGap = 7) 
eSS.t <- detect_event(data = tMax_clim.SS, y = tMax,minDuration = 5,maxGap = 7) 
eSW.t <- detect_event(data = tMax_clim.SW, y = tMax,minDuration = 5,maxGap = 7) 


eBS.sic <- detect_event(data = sicMean_clim.BS, y = sicMean,minDuration = 3,maxGap = 7,coldSpells = T) 
eEI.sic <- detect_event(data = sicMean_clim.EI, y = sicMean,minDuration = 3,maxGap = 7,coldSpells = T) 
eEX.sic <- detect_event(data = sicMean_clim.EX, y = sicMean,minDuration = 3,maxGap = 7,coldSpells = T) 
eJI.sic <- detect_event(data = sicMean_clim.JI, y = sicMean,minDuration = 3,maxGap = 7,coldSpells = T) 
eSS.sic <- detect_event(data = sicMean_clim.SS, y = sicMean,minDuration = 3,maxGap = 7,coldSpells = T) 
eSW.sic <- detect_event(data = sicMean_clim.SW, y = sicMean,minDuration = 3,maxGap = 7,coldSpells = T) 


#export climatologies and create dataframes

bsclimT<-data.frame(eBS.t$climatology,sector=c("BS")) 
elclimT<-data.frame(eEI.t$climatology,sector=c("EI")) 
exclimT<-data.frame(eEX.t$climatology,sector=c("EX")) 
jiclimT<-data.frame(eJI.t$climatology,sector=c("JI")) 
ssclimT<-data.frame(eSS.t$climatology,sector=c("SS")) 
swclimT<-data.frame(eSW.t$climatology,sector=c("SW")) 

bsclimS<-data.frame(eBS.sic$climatology,sector=c("BS")) 
elclimS<-data.frame(eEI.sic$climatology,sector=c("EI")) 
exclimS<-data.frame(eEX.sic$climatology,sector=c("EX")) 
jiclimS<-data.frame(eJI.sic$climatology,sector=c("JI")) 
ssclimS<-data.frame(eSS.sic$climatology,sector=c("SS")) 
swclimS<-data.frame(eSW.sic$climatology,sector=c("SW")) 

climT<-rbind(bsclimT,elclimT,exclimT,jiclimT,ssclimT,swclimT) # join data
climS<-rbind(bsclimS,elclimS,exclimS,jiclimS,ssclimS,swclimS) # join data


#climT$Day<-(decimal_date(climT$date_peak)-year(climT$date_peak))*365 # day of the year
climT$quart<-quarter(climT$t) # tri-mensual classification

#climS$Day<-(decimal_date(climS$date_peak)-year(climS$date_peak))*365 # day of the year
climS$quart<-quarter(climS$t) # tri-mensual classification


# change names
climT$Season[climT$quart=="1"]<-"JFM"
climT$Season[climT$quart=="2"]<-"AMJ"
climT$Season[climT$quart=="3"]<-"JAS"
climT$Season[climT$quart=="4"]<-"OND"

climS$Season[climS$quart=="1"]<-"JFM"
climS$Season[climS$quart=="2"]<-"AMJ"
climS$Season[climS$quart=="3"]<-"JAS"
climS$Season[climS$quart=="4"]<-"OND"





# plots for temperature and sea ice

# Create category breaks and select slice of data.frame


# Set line colours
lineColCat <- c(
  "Temperature" = "black",
  "Climatology" = "gray20",
  "Threshold" = "darkgreen",
  "2x Threshold" = "darkgreen",
  "3x Threshold" = "darkgreen",
  "4x Threshold" = "darkgreen"
  )

# Set category fill colours
fillColCat <- c(
  "Moderate" = "#ffc866",
  "Strong" = "#ff6900",
  "Severe" = "#9e0000",
  "Extreme" = "#2d0000"
  )

# Set category fill colours
fillColCat2 <- c(
  "Moderate" = "#ffc866",
  "Normal" = "#85B7CC"
  )

lineColCat2 <- c(
  "Sea Ice" = "black",
  "Climatology" = "blue3",
  "Threshold" = "red3"
  )
# Create category breaks 
climT
climT$diff = climT$thresh - climT$seas
climT$thresh_2x = climT$thresh + climT$diff
climT$thresh_3x = climT$thresh_2x + climT$diff
climT$thresh_4x = climT$thresh_3x + climT$diff



saveRDS(climS,"climatologyICE.Rds")
saveRDS(climT,"climatologyTEMP.Rds")


#figure 
(ggplot(data = subset(climT,t>'2022-03-31'), aes(x = t, y = tMax)) +
  geom_flame(aes(y2 = thresh, fill = "Moderate")) +
  geom_flame(aes(y2 = thresh_2x, fill = "Strong")) +
  geom_flame(aes(y2 = thresh_3x, fill = "Severe")) +
  geom_flame(aes(y2 = thresh_4x, fill = "Extreme")) +
  geom_line(aes(y = thresh_2x, col = "2x Threshold"), size = 0.7, linetype = "dashed") +
  geom_line(aes(y = thresh_3x, col = "3x Threshold"), size = 0.7, linetype = "dotdash") +
  geom_line(aes(y = thresh_4x, col = "4x Threshold"), size = 0.7, linetype = "dotted") +
  geom_line(aes(y = seas, col = "Climatology"), size = 0.7) +
  geom_line(aes(y = thresh, col = "Threshold"), size = 0.7) +
  geom_line(aes(y = tMax, col = "Temperature"), size = 0.6) +
  scale_colour_manual(name = NULL, values = lineColCat,
                      breaks = c("Temperature", "Climatology", "Threshold",
                                 "2x Threshold", "3x Threshold", "4x Threshold")) +
  scale_fill_manual(name = NULL, values = fillColCat, guide = FALSE) +
  scale_x_date(date_labels = "%b %Y") +
  guides(colour = guide_legend(override.aes = list(linetype = c("solid", "solid", "solid",
                                                                "dashed", "dotdash", "dotted"),
                                                   size = c(0.6, 0.7, 0.7, 0.7, 0.7, 0.7)))) +
  labs(y = "Temperature [°C]", x = NULL)+
  facet_wrap(sector~.)+theme_bw()+theme(panel.spacing.x = unit(2, "lines"))+
  ggtitle(label="a. High temperature extremes"))
  
  
  

#export events and creat data frame

bsevenT<-data.frame(eBS.t$event,sector=c("BS")) 
elevenT<-data.frame(eEI.t$event,sector=c("EI")) 
exevenT<-data.frame(eEX.t$event,sector=c("EX")) 
jievenT<-data.frame(eJI.t$event,sector=c("JI")) 
ssevenT<-data.frame(eSS.t$event,sector=c("SS")) 
swevenT<-data.frame(eSW.t$event,sector=c("SW")) 


evenT<-rbind(bsevenT,elevenT,exevenT,jievenT,ssevenT,swevenT) # join events data

evenT$Day<-(decimal_date(evenT$date_peak)-year(evenT$date_peak))*365 # day of the year
evenT$quart<-quarter(evenT$date_peak) # tri-mensual classification


# change names
evenT$Season[evenT$quart=="1"]<-"JFM"
evenT$Season[evenT$quart=="2"]<-"AMJ"
evenT$Season[evenT$quart=="3"]<-"JAS"
evenT$Season[evenT$quart=="4"]<-"OND"

evenT$Season<-factor(evenT$Season,levels=c("JFM","AMJ","JAS","OND")) # order levels


#PLOT figure 2



(ggplot(evenT,aes(date_peak,duration,colour=Season,linetype=Season))+
    geom_lolli()+
    theme_bw()+ylab("Duration (days)")+
    scale_colour_manual(values=c("red","green4","blue","orange"))+
    facet_wrap(sector~.)+
    xlab("Date peak")+
    ggtitle(label="a. Duration of high temperature events"))+
  
  (ggplot(evenT,aes(date_peak,intensity_max,colour=Season,linetype=Season))+
     geom_lolli()+
     theme_bw()+ylab("Intensity (°C)")+
     scale_colour_manual(values=c("red","green4","blue","orange"))+
     facet_wrap(sector~.)+
     xlab("Date peak")+
     ggtitle(label="b. Intensity of temperature events"))
  
head(evenT)
evenT$Year<-year(evenT$date_start)




#summarize per season and year

evesum<-plyr::ddply(evenT, c("Year","Season","sector"), summarise,
            TN=length(duration),
            TmeanDur=mean(duration),
            TcumDur=sum(duration),
            TmeanInt=mean(intensity_mean),
            TmaxInt=max(intensity_max))



saveRDS(evesum,"events_summary.Rds")

summary(evesum)

evesum$TN[is.na(evesum$TN)]<-0
evesum$TmeanDur[is.na(evesum$TmeanDur)]<-0
evesum$TcumDur[is.na(evesum$TcumDur)]<-0
evesum$TmeanInt[is.na(evesum$TmeanInt)]<-0
evesum$TmaxInt[is.na(evesum$TmaxInt)]<-0

head(evesum)


meanT<-plyr::ddply(evenT, c("Year","sector"), summarise,
                   TNt=length(duration),
                   TDur=sum(duration),
                   Tint=mean(intensity_mean),
                   Tintx=mean(intensity_max))
meanT$sector[meanT$sector=="EX"]<-"GS"

summary(meanT$TNt)

lm1<-lm(TDur~Year,data=meanT)

summary(lm1)

lm2<-lm(Tintx~Year,data=meanT)

summary(lm2)



meanS<-plyr::ddply(climS, c("Year","sector"), summarise,
                  sicm=mean(sicMean))



meanTS<-merge(meanT,meanS)

meanTS$lat[meanTS$sector=="BS"|meanTS$sector=="SS"|meanTS$sector=="EI"]<-"N"
meanTS$lat[is.na(meanTS$lat)]<-"S"
  




lmDurN<-lm(sicm~log(TDur),data=subset(meanTS,lat=="N"))

summary(lmDurN)


lmDurS<-lm(sicm~log(TDur),data=subset(meanTS,lat=="S"))

summary(lmDurS)


lmIntN<-lm(sicm~Tint,data=subset(meanTS,lat=="N"))

summary(lmIntN)


lmIntS<-lm(sicm~Tint,data=subset(meanTS,lat=="S"))

summary(lmIntS)


(ggplot()+
    geom_smooth(method="lm",aes(x=Year,y=TDur),se=F,data=meanT)+
    geom_point(data=meanT,aes(x=Year,y=TDur,colour=sector,shape=sector))+
    theme_bw()+ylab("days above threshold")+#scale_y_log10()+
    scale_colour_manual(values=c("red4","red2","blue2","blue","red","blue4"))+
    
    #facet_wrap(sector~.)+
    xlab("Year")+
    ggtitle(label="a."))+
  
  
  (ggplot()+
     geom_smooth(method="lm",aes(x=Year,y=Tintx),se=F,data=meanT)+
     geom_point(data=meanT,aes(x=Year,y=Tintx,colour=sector,shape=sector))+
     theme_bw()+ylab("°C above threshold")+
     scale_colour_manual(values=c("red4","red2","blue2","blue","red","blue4"))+
     #facet_wrap(sector~.)+
     xlab("Year")+
     ggtitle(label="b."))+

ggplot(meanTS,aes(log(TDur),sicm,colour=lat,shape=lat,linetpe=lat))+
  stat_smooth(method="glm", method.args = list(family = "binomial"),se=F)+
  geom_point()+theme_bw()+scale_colour_manual(values=c("red2","blue2"))+
  xlab("Duration of MHWs (log days)")+
  ylab("Sea ice cover (proportion)")+
  ggtitle(label="c.")+
  
  ggplot(meanTS,aes(Tint,sicm,colour=lat,shape=lat,linetpe=lat))+
  stat_smooth(method="glm", method.args = list(family = "binomial"),se=F)+
  geom_point()+theme_bw()+scale_colour_manual(values=c("red2","blue2"))+
  xlab("Intensity of MHWs (°C)")+
  ylab("")+
  ggtitle(label="d.")





```