Hundley et al temperature analysis

Author

Justin Baumann

Load Packages

#load packages
library(tidyverse)
library(lubridate)
library(tibbletime)
#devtools::install_github("dbca-wa/dbcaDHW")
library(dbcaDHW) #DHW and other CRW metrics
library(heatwaveR)
library(patchwork)
library(rerddap)
library(scales)
library(knitr)

Download SST Data from ERDDAP –> ONLY DO THIS ONCE

sst_info<-rerddap::info(datasetid = "jplMURSST42", url = "https://coastwatch.pfeg.noaa.gov/erddap/") #datasetid should come straight from ERDDAP

#False MUR --> NOTE that this is QUITE a bit south of our site, but the nearest NS pixel we can get that isn't on land!
false_murSST = griddap(sst_info, 
                       latitude = c(16.375, 16.375),
                       longitude = c(-88.375,-88.375),
                       time=c("2002-09-02","2024-05-18"), 
                       fmt = "csv")
write.csv(false_murSST, file='false_MURSST.csv')

#Silk MUR --> THIS IS REALLY close to the actual site!
silk_murSST = griddap(sst_info, 
                       latitude = c(16.451775, 16.451775),
                       longitude = c(-88.045393,-88.045393),
                       time=c("2002-09-02","2024-05-18"), 
                       fmt = "csv")
write.csv(silk_murSST, file='silk_murSST.csv')

Read in Data

false<-read_csv("https://raw.githubusercontent.com/jbaumann3/Hundley_et_al_Belize_RT_5_year/main/foh_false_2018_2023.csv") #insitu

false_MURSST<-read_csv("https://raw.githubusercontent.com/jbaumann3/Hundley_et_al_Belize_RT_5_year/main/false_MURSST.csv") #satellite

silk<-read.csv("https://raw.githubusercontent.com/jbaumann3/Hundley_et_al_Belize_RT_5_year/main/foh_silk_2015_2023.csv") #insitu

silk_murSST<-read_csv("https://raw.githubusercontent.com/jbaumann3/Hundley_et_al_Belize_RT_5_year/main/silk_murSST.csv")#satellite

prepare in-situ data

make silk & False date into a date object

silk$date<- mdy_hm(silk$date)
false$date<- mdy_hm(false$date)

tail(silk)
                     date Temp_C year month day hour min sec
76644 2024-06-10 18:00:00 30.976 2024     6  10   18   0   0
76645 2024-06-10 19:00:00 30.900 2024     6  10   19   0   0
76646 2024-06-10 20:00:00 30.722 2024     6  10   20   0   0
76647 2024-06-10 21:00:00 30.722 2024     6  10   21   0   0
76648 2024-06-10 22:00:00 30.722 2024     6  10   22   0   0
76649 2024-06-10 23:00:00 30.798 2024     6  10   23   0   0
tail(false)
# A tibble: 6 × 8
  date                Temp_C  year month   day  hour   min   sec
  <dttm>               <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2024-03-29 12:00:00   26.5  2024     3    29    12     0     0
2 2024-03-29 13:00:00   30.1  2024     3    29    13     0     0
3 2024-03-29 14:00:00   30.4  2024     3    29    14     0     0
4 2024-03-29 15:00:00   23.8  2024     3    29    15     0     0
5 2024-03-29 16:00:00   23.2  2024     3    29    16     0     0
6 2024-03-29 17:00:00   23.4  2024     3    29    17     0     0

add site column and combine data

false$site<-'False Caye'
silk$site<-'Silk Caye'

bind together and save as .csv

tempdat<-rbind(false,silk)
#write.csv(tempdat,file='bz_rt_5yr_temp.csv')

calculate daily average temps to match 1 measurement per day SST data and preferred MHW and DHW calcs

daily<- tempdat %>%
  mutate(day=as_date(date)) %>%
  group_by(site, day) %>%
  drop_na(date) %>%
  summarize(meandailytemp=mean(Temp_C), sd=sd(Temp_C), n=n(), se=sd/sqrt(n))
`summarise()` has grouped output by 'site'. You can override using the
`.groups` argument.
falseday<-daily %>%
  filter(site=='False Caye')

silkday<-daily %>%
  filter(site=='Silk Caye')

prepare satellite data

make date column a date object

false_MURSST
# A tibble: 7,724 × 9
    ...1 time                latitude longitude analysed_sst analysis_error
   <dbl> <dttm>                 <dbl>     <dbl>        <dbl>          <dbl>
 1     1 2002-09-02 09:00:00     16.4     -88.4         29.4           0.39
 2     2 2002-09-03 09:00:00     16.4     -88.4         29.4           0.39
 3     3 2002-09-04 09:00:00     16.4     -88.4         29.4           0.39
 4     4 2002-09-05 09:00:00     16.4     -88.4         29.8           0.4 
 5     5 2002-09-06 09:00:00     16.4     -88.4         29.9           0.39
 6     6 2002-09-07 09:00:00     16.4     -88.4         30.1           0.39
 7     7 2002-09-08 09:00:00     16.4     -88.4         30.4           0.37
 8     8 2002-09-09 09:00:00     16.4     -88.4         30.5           0.36
 9     9 2002-09-10 09:00:00     16.4     -88.4         30.4           0.36
10    10 2002-09-11 09:00:00     16.4     -88.4         30.2           0.37
# ℹ 7,714 more rows
# ℹ 3 more variables: mask <dbl>, sea_ice_fraction <lgl>, sst_anomaly <dbl>
false_MURSST<- false_MURSST %>%
  select(-c(...1,sea_ice_fraction, mask)) %>%
  mutate(date=ymd_hms(time))

silk_MURSST<- silk_murSST %>%
  select(-c(...1,sea_ice_fraction, mask)) %>%
  mutate(date=ymd_hms(time))
silk_MURSST
# A tibble: 7,724 × 7
   time                latitude longitude analysed_sst analysis_error
   <dttm>                 <dbl>     <dbl>        <dbl>          <dbl>
 1 2002-09-02 09:00:00     16.4     -88.1         29.4           0.38
 2 2002-09-03 09:00:00     16.4     -88.1         29.5           0.37
 3 2002-09-04 09:00:00     16.4     -88.1         29.6           0.37
 4 2002-09-05 09:00:00     16.4     -88.1         29.8           0.38
 5 2002-09-06 09:00:00     16.4     -88.1         29.8           0.37
 6 2002-09-07 09:00:00     16.4     -88.1         29.9           0.37
 7 2002-09-08 09:00:00     16.4     -88.1         30.0           0.36
 8 2002-09-09 09:00:00     16.4     -88.1         29.8           0.36
 9 2002-09-10 09:00:00     16.4     -88.1         29.5           0.36
10 2002-09-11 09:00:00     16.4     -88.1         29.7           0.38
# ℹ 7,714 more rows
# ℹ 2 more variables: sst_anomaly <dbl>, date <dttm>

add site column and combine data

false_MURSST$site<-'False Caye'
silk_MURSST$site<-'Silk Caye'

##combine datasets into 1

satdat<-rbind(false_MURSST,silk_MURSST) %>%
  drop_na(analysed_sst) %>%
  mutate(analysed_sst = as.numeric(analysed_sst)) %>%
  mutate(analysis_error=as.numeric(analysis_error)) %>%
  mutate(sst_anomaly=as.numeric(sst_anomaly))

Comparing Sat and in-situ data

#trim and combine
satdat$measurement <- "Satellite JPL MUR"
daily$measurement <- "in-situ"

lilsat<- satdat %>%
  mutate(date=as.Date(date)) %>%
  select(date, analysed_sst, measurement, site) %>%
  mutate(temp=analysed_sst) %>%
  select(-analysed_sst)

liltemp<- daily %>%
  select(day, meandailytemp, measurement, site) %>%
  mutate(temp=meandailytemp, date=day) %>%
  select(-c(meandailytemp, day))

alltemp<- rbind(lilsat, liltemp)


satandinsitu<-ggplot(data=alltemp, aes(x=date, y=temp, color=measurement))+
  geom_line(linewidth=0.75)+
  theme_classic()+
  scale_color_manual(values=c('firebrick2','black'))+
  scale_x_date(date_breaks = "6 months", limits=ymd(c('2017-12-01', '2024-01-01')), expand=c(0,0))+
  theme(axis.text.x=element_text(angle=90), text=element_text(size=12))+
  facet_wrap(~site, ncol=1)

satandinsitu
Warning: Removed 6623 rows containing missing values or values outside the scale range
(`geom_line()`).

Difference between sat and in-situ daily values

widetemp<- alltemp %>%
  pivot_wider(names_from = measurement, values_from = temp) %>%
  drop_na(c(`Satellite JPL MUR`, `in-situ`)) %>%
  mutate(diff= `in-situ`-`Satellite JPL MUR`)
widetemp
# A tibble: 5,115 × 5
   date       site       `Satellite JPL MUR` `in-situ`    diff
   <date>     <chr>                    <dbl>     <dbl>   <dbl>
 1 2018-12-05 False Caye                28.2      28.8  0.622 
 2 2018-12-06 False Caye                28.5      28.6  0.105 
 3 2018-12-07 False Caye                28.4      28.5  0.106 
 4 2018-12-08 False Caye                28.1      28.5  0.380 
 5 2018-12-09 False Caye                28.0      28.6  0.535 
 6 2018-12-10 False Caye                28.0      28.3  0.357 
 7 2018-12-11 False Caye                27.9      27.6 -0.268 
 8 2018-12-12 False Caye                27.8      27.1 -0.636 
 9 2018-12-13 False Caye                27.4      27.5  0.0770
10 2018-12-14 False Caye                27.5      27.8  0.364 
# ℹ 5,105 more rows
meandiff<- widetemp %>%
  summarize(mean= mean(diff), sd=sd(diff), n=n(), se=sd/sqrt(n))

meandiff
# A tibble: 1 × 4
   mean    sd     n      se
  <dbl> <dbl> <int>   <dbl>
1 0.392 0.491  5115 0.00687
#In-situ is ~ 0.391 +/- 0.00728 (SEM) higher than satellite across this time period

satvsinsitu<-ggplot(widetemp, aes(x=date, y=diff, color=site))+
  geom_line()+
  geom_hline(yintercept=0, linetype='dashed')+
  theme_classic()+
  facet_wrap(~site, ncol=1)+
  scale_color_manual(values=c('firebrick','dodgerblue3'))+
  labs(x=NULL, y='In-situ temp - SST')

satvsinsitu

DHW using dbcaDHW

#Trying max monthly means using our satellite data as past dataset (pre-2020)

head(satdat)
# A tibble: 6 × 9
  time                latitude longitude analysed_sst analysis_error sst_anomaly
  <dttm>                 <dbl>     <dbl>        <dbl>          <dbl>       <dbl>
1 2002-09-02 09:00:00     16.4     -88.4         29.4           0.39      -0.5  
2 2002-09-03 09:00:00     16.4     -88.4         29.4           0.39      -0.518
3 2002-09-04 09:00:00     16.4     -88.4         29.4           0.39      -0.439
4 2002-09-05 09:00:00     16.4     -88.4         29.8           0.4       -0.025
5 2002-09-06 09:00:00     16.4     -88.4         29.9           0.39      -0.109
6 2002-09-07 09:00:00     16.4     -88.4         30.1           0.39       0.088
# ℹ 3 more variables: date <dttm>, site <chr>, measurement <chr>
mmtdata <- satdat %>%
  #tidyr::gather("site", "Temp_C", 2:all_of(g_num)) %>%
  dplyr::filter(date < "2017-12-01") %>% #define past (pre-experiment!)
  dplyr::mutate(month = month(date),
                year = year(date)) %>%
  dplyr::group_by(site, year, month) %>%
  drop_na(date, analysed_sst) %>%
  dplyr::summarise(avg = mean(analysed_sst)) %>% #obtain mthly avgs
  dplyr::group_by(site) %>%
  dplyr::summarise(mmmt = max(avg)) #obtain max mthly avg
`summarise()` has grouped output by 'site', 'year'. You can override using the
`.groups` argument.
mmtdata #shows average MMM for each site from 2002-2020 (18 years-- length of dataset)
# A tibble: 2 × 2
  site        mmmt
  <chr>      <dbl>
1 False Caye  30.4
2 Silk Caye   30.3
#False=30.36
#Silk=30.33

DHW calcs

## function for rolling cusum of 12 weeks - 84 days)
dhw_calc_12 <- tibbletime::rollify(sum, window = 84)

observation period, hotspot & dhw calc USING TEMPDAT (in-situ) with baselines (MMMT) from satellite

## observation period, hotspot & dhw calc
dhw <- tempdat %>%
  #tidyr::gather("site", "sst", 2:all_of(g_num)) %>%
  dplyr::filter(date >= "2017-12-01") %>% #define present
  dplyr::mutate(month = month(date)) %>%
  dplyr::full_join(mmtdata, by = "site") %>%
  dplyr::mutate(hspt = Temp_C - mmmt,
                hsptm = ifelse(hspt >= 1, hspt, 0),
                dhw = dhw_calc_12(hsptm)*(1/7)) %>% #hotspot based on 1 degree
  dplyr::select(-month, -hsptm)


dhw
# A tibble: 100,338 × 11
   date                Temp_C  year   day  hour   min   sec site      mmmt  hspt
   <dttm>               <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>    <dbl> <dbl>
 1 2018-12-05 05:00:00   28.5  2018     5     4    59    59 False C…  30.4 -1.84
 2 2018-12-05 06:00:00   28.5  2018     5     6     0     0 False C…  30.4 -1.89
 3 2018-12-05 07:00:00   28.5  2018     5     7     0     0 False C…  30.4 -1.82
 4 2018-12-05 08:00:00   28.5  2018     5     7    59    59 False C…  30.4 -1.87
 5 2018-12-05 09:00:00   28.5  2018     5     9     0     0 False C…  30.4 -1.82
 6 2018-12-05 10:00:00   28.5  2018     5    10     0     0 False C…  30.4 -1.84
 7 2018-12-05 11:00:00   28.5  2018     5    10    59    59 False C…  30.4 -1.84
 8 2018-12-05 12:00:00   28.5  2018     5    12     0     0 False C…  30.4 -1.87
 9 2018-12-05 13:00:00   28.4  2018     5    13     0     0 False C…  30.4 -1.92
10 2018-12-05 14:00:00   28.6  2018     5    13    59    59 False C…  30.4 -1.74
# ℹ 100,328 more rows
# ℹ 1 more variable: dhw <dbl>
tail(dhw)
# A tibble: 6 × 11
  date                Temp_C  year   day  hour   min   sec site       mmmt  hspt
  <dttm>               <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>     <dbl> <dbl>
1 2024-06-10 18:00:00   31.0  2024    10    18     0     0 Silk Caye  30.3 0.644
2 2024-06-10 19:00:00   30.9  2024    10    19     0     0 Silk Caye  30.3 0.568
3 2024-06-10 20:00:00   30.7  2024    10    20     0     0 Silk Caye  30.3 0.390
4 2024-06-10 21:00:00   30.7  2024    10    21     0     0 Silk Caye  30.3 0.390
5 2024-06-10 22:00:00   30.7  2024    10    22     0     0 Silk Caye  30.3 0.390
6 2024-06-10 23:00:00   30.8  2024    10    23     0     0 Silk Caye  30.3 0.466
# ℹ 1 more variable: dhw <dbl>

#max DHW

Table of DHW by year of the experiment

dhwm <- dhw %>%
  dplyr::group_by(site, year) %>%
  dplyr::summarise(dhwm = max(dhw, na.rm = TRUE)) %>%
  kable()
`summarise()` has grouped output by 'site'. You can override using the
`.groups` argument.
dhwm
site year dhwm
False Caye 2018 0.000000
False Caye 2019 11.256848
False Caye 2020 18.358924
False Caye 2021 9.566876
False Caye 2022 4.597010
False Caye 2023 24.806829
False Caye 2024 6.214305
Silk Caye 2017 6.214305
Silk Caye 2018 0.000000
Silk Caye 2019 6.207591
Silk Caye 2020 7.242629
Silk Caye 2021 1.228390
Silk Caye 2022 0.376581
Silk Caye 2023 9.643905
Silk Caye 2024 0.000000

Graph of DHW by year of the experiment

dhwgraph<-ggplot(dhw, aes(x=date, y=dhw, color=site))+
  geom_line(size=0.5)+
  theme_classic()+
  scale_color_manual(values=c('firebrick','dodgerblue3'))+
  scale_x_datetime(date_breaks = "3 months")+
  theme(axis.text.x=element_text(angle=90), text=element_text(size=12))+
  facet_wrap(~site, ncol=1)+
  geom_hline(yintercept=4, linetype='dotdash', color='black')+
  geom_hline(yintercept=8, linetype='dashed', color='red')+
  labs(x=NULL, y='Degree Heating Weeks')
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
dhwgraph
Warning: Removed 83 rows containing missing values or values outside the scale range
(`geom_line()`).

heatwaveR MHW and climatology calcs

Part 1: Climatologies

-data should include 2 columns (datetime (t) and temp (temp)) -one site at a time for this

-Starting with the satellite data!

reshape the data as required

### SATELLITE
head(false_MURSST)
# A tibble: 6 × 8
  time                latitude longitude analysed_sst analysis_error sst_anomaly
  <dttm>                 <dbl>     <dbl>        <dbl>          <dbl>       <dbl>
1 2002-09-02 09:00:00     16.4     -88.4         29.4           0.39      -0.5  
2 2002-09-03 09:00:00     16.4     -88.4         29.4           0.39      -0.518
3 2002-09-04 09:00:00     16.4     -88.4         29.4           0.39      -0.439
4 2002-09-05 09:00:00     16.4     -88.4         29.8           0.4       -0.025
5 2002-09-06 09:00:00     16.4     -88.4         29.9           0.39      -0.109
6 2002-09-07 09:00:00     16.4     -88.4         30.1           0.39       0.088
# ℹ 2 more variables: date <dttm>, site <chr>
head(silk_MURSST)
# A tibble: 6 × 8
  time                latitude longitude analysed_sst analysis_error sst_anomaly
  <dttm>                 <dbl>     <dbl>        <dbl>          <dbl>       <dbl>
1 2002-09-02 09:00:00     16.4     -88.1         29.4           0.38      -0.377
2 2002-09-03 09:00:00     16.4     -88.1         29.5           0.37      -0.179
3 2002-09-04 09:00:00     16.4     -88.1         29.6           0.37      -0.027
4 2002-09-05 09:00:00     16.4     -88.1         29.8           0.38       0.038
5 2002-09-06 09:00:00     16.4     -88.1         29.8           0.37       0.011
6 2002-09-07 09:00:00     16.4     -88.1         29.9           0.37       0.046
# ℹ 2 more variables: date <dttm>, site <chr>
false_clim<- false_MURSST[-1,] %>%
  mutate(date=as.Date(date)) %>%
  mutate(temp=as.numeric(analysed_sst)) %>%
  select(date, temp)
false_clim
# A tibble: 7,723 × 2
   date        temp
   <date>     <dbl>
 1 2002-09-03  29.4
 2 2002-09-04  29.4
 3 2002-09-05  29.8
 4 2002-09-06  29.9
 5 2002-09-07  30.1
 6 2002-09-08  30.4
 7 2002-09-09  30.5
 8 2002-09-10  30.4
 9 2002-09-11  30.2
10 2002-09-12  29.9
# ℹ 7,713 more rows
silk_clim<-silk_MURSST %>%
  mutate(date=as.Date(date)) %>%
  mutate(temp=analysed_sst) %>%
  select(date, temp)
silk_clim
# A tibble: 7,724 × 2
   date        temp
   <date>     <dbl>
 1 2002-09-02  29.4
 2 2002-09-03  29.5
 3 2002-09-04  29.6
 4 2002-09-05  29.8
 5 2002-09-06  29.8
 6 2002-09-07  29.9
 7 2002-09-08  30.0
 8 2002-09-09  29.8
 9 2002-09-10  29.5
10 2002-09-11  29.7
# ℹ 7,714 more rows
###IN SITU
head(falseday)
# A tibble: 6 × 6
# Groups:   site [1]
  site       day        meandailytemp    sd     n     se
  <chr>      <date>             <dbl> <dbl> <int>  <dbl>
1 False Caye 2018-12-05          28.8 0.374    19 0.0858
2 False Caye 2018-12-06          28.6 0.284    24 0.0580
3 False Caye 2018-12-07          28.5 0.345    24 0.0704
4 False Caye 2018-12-08          28.5 0.440    24 0.0897
5 False Caye 2018-12-09          28.6 0.447    24 0.0913
6 False Caye 2018-12-10          28.3 0.179    24 0.0366
head(silkday)
# A tibble: 6 × 6
# Groups:   site [1]
  site      day        meandailytemp    sd     n     se
  <chr>     <date>             <dbl> <dbl> <int>  <dbl>
1 Silk Caye 2015-04-18          29.3 0.591     6 0.241 
2 Silk Caye 2015-04-19          29.0 0.580    24 0.118 
3 Silk Caye 2015-04-20          29.1 0.567    24 0.116 
4 Silk Caye 2015-04-21          29.3 0.580    24 0.118 
5 Silk Caye 2015-04-22          29.5 0.489    24 0.0999
6 Silk Caye 2015-04-23          29.5 0.565    24 0.115 
false_isclim<- falseday[,-1] %>%
  mutate(date=as.Date(day)) %>%
  mutate(temp=meandailytemp) %>%
  select(date, temp)
false_isclim
# A tibble: 1,942 × 2
   date        temp
   <date>     <dbl>
 1 2018-12-05  28.8
 2 2018-12-06  28.6
 3 2018-12-07  28.5
 4 2018-12-08  28.5
 5 2018-12-09  28.6
 6 2018-12-10  28.3
 7 2018-12-11  27.6
 8 2018-12-12  27.1
 9 2018-12-13  27.5
10 2018-12-14  27.8
# ℹ 1,932 more rows
silk_isclim<- silkday[,-1] %>%
  mutate(date=as.Date(day)) %>%
  mutate(temp=meandailytemp) %>%
  select(date, temp)
silk_isclim
# A tibble: 3,203 × 2
   date        temp
   <date>     <dbl>
 1 2015-04-18  29.3
 2 2015-04-19  29.0
 3 2015-04-20  29.1
 4 2015-04-21  29.3
 5 2015-04-22  29.5
 6 2015-04-23  29.5
 7 2015-04-24  29.4
 8 2015-04-25  29.4
 9 2015-04-26  29.5
10 2015-04-27  29.6
# ℹ 3,193 more rows

Build climatolgies

With sat data (2002-present) the climatology is 20+ years!  -In-situ data are limited to only last 5 years (more for Silk)

#use the ts2clm function

tssilk<-ts2clm(silk_clim, x=date, y=temp, climatologyPeriod = c("2002-09-02", "2024-05-18"),var=TRUE)
tssilk
# A tibble: 7,930 × 6
     doy date        temp  seas thresh   var
   <int> <date>     <dbl> <dbl>  <dbl> <dbl>
 1   246 2002-09-02  29.4  29.8   30.4 0.475
 2   247 2002-09-03  29.5  29.8   30.4 0.476
 3   248 2002-09-04  29.6  29.9   30.4 0.477
 4   249 2002-09-05  29.8  29.9   30.5 0.478
 5   250 2002-09-06  29.8  29.9   30.5 0.480
 6   251 2002-09-07  29.9  29.9   30.5 0.482
 7   252 2002-09-08  30.0  29.9   30.5 0.485
 8   253 2002-09-09  29.8  29.9   30.5 0.487
 9   254 2002-09-10  29.5  29.9   30.5 0.488
10   255 2002-09-11  29.7  29.9   30.5 0.490
# ℹ 7,920 more rows
tsfalse<-ts2clm(false_clim, x=date, y=temp, climatologyPeriod = c("2018-12-05", "2024-05-18"),var=TRUE)
tsfalse
# A tibble: 7,929 × 6
     doy date        temp  seas thresh   var
   <int> <date>     <dbl> <dbl>  <dbl> <dbl>
 1   247 2002-09-03  29.4  30.2   30.6 0.348
 2   248 2002-09-04  29.4  30.2   30.6 0.346
 3   249 2002-09-05  29.8  30.2   30.6 0.345
 4   250 2002-09-06  29.9  30.2   30.6 0.346
 5   251 2002-09-07  30.1  30.2   30.6 0.347
 6   252 2002-09-08  30.4  30.2   30.6 0.350
 7   253 2002-09-09  30.5  30.2   30.6 0.354
 8   254 2002-09-10  30.4  30.2   30.6 0.358
 9   255 2002-09-11  30.2  30.2   30.6 0.362
10   256 2002-09-12  29.9  30.2   30.6 0.367
# ℹ 7,919 more rows

ts2clim returns the following columns:
-temp (as per original data)
-seas (daily climatological cycle in degree C)
-thresh (daily varying threshold (90th percentile is default) that is used in detect_event to ID a MHW)
-var: daily variance (default is var=FALSE, which does not report var)

Graphs of the variables for each site

silkclimgraph<- ggplot(data=tssilk, aes(x=date))+
  #geom_line(aes(y=temp),color='black', alpha=0.5, linetype='dashed')+
  geom_line(aes(y=seas), color='dodgerblue3')+
  geom_line(aes(y=thresh), color='firebrick', size=1)+
  #geom_line(aes(y=var), color='orange3')+
  theme_bw()
silkclimgraph

falseclimgraph<- ggplot(data=tsfalse, aes(x=date))+
  #geom_line(aes(y=temp),color='black', alpha=0.5, linetype='dashed')+
  geom_line(aes(y=seas), color='dodgerblue3')+
  geom_line(aes(y=thresh), color='firebrick', size=1)+
  #geom_line(aes(y=var), color='orange3')+
  theme_bw()
falseclimgraph

falseclimgraph+silkclimgraph

detect MHW using detect_event

Note: It is preferred to use 30 year climatologies for this, we only have data form 2002-present!

#S = FALSE tells the function we are in the northern hemisphere

silk_event<-heatwaveR::detect_event(tssilk, x=date, y=temp, S=FALSE)
silk_event
$climatology
# A tibble: 7,930 × 10
     doy date        temp  seas thresh   var threshCriterion durationCriterion
   <int> <date>     <dbl> <dbl>  <dbl> <dbl> <lgl>           <lgl>            
 1   246 2002-09-02  29.4  29.8   30.4 0.475 FALSE           FALSE            
 2   247 2002-09-03  29.5  29.8   30.4 0.476 FALSE           FALSE            
 3   248 2002-09-04  29.6  29.9   30.4 0.477 FALSE           FALSE            
 4   249 2002-09-05  29.8  29.9   30.5 0.478 FALSE           FALSE            
 5   250 2002-09-06  29.8  29.9   30.5 0.480 FALSE           FALSE            
 6   251 2002-09-07  29.9  29.9   30.5 0.482 FALSE           FALSE            
 7   252 2002-09-08  30.0  29.9   30.5 0.485 FALSE           FALSE            
 8   253 2002-09-09  29.8  29.9   30.5 0.487 FALSE           FALSE            
 9   254 2002-09-10  29.5  29.9   30.5 0.488 FALSE           FALSE            
10   255 2002-09-11  29.7  29.9   30.5 0.490 FALSE           FALSE            
# ℹ 7,920 more rows
# ℹ 2 more variables: event <lgl>, event_no <int>

$event
# A tibble: 64 × 22
   event_no index_start index_peak index_end duration date_start date_peak 
      <int>       <int>      <int>     <int>    <dbl> <date>     <date>    
 1        1        1012       1014      1017        6 2005-06-09 2005-06-11
 2        2        1040       1042      1055       16 2005-07-07 2005-07-09
 3        3        1091       1096      1097        7 2005-08-27 2005-09-01
 4        4        1604       1606      1608        5 2007-01-22 2007-01-24
 5        5        1613       1617      1620        8 2007-01-31 2007-02-04
 6        6        1624       1627      1630        7 2007-02-11 2007-02-14
 7        7        1642       1644      1647        6 2007-03-01 2007-03-03
 8        8        1685       1687      1689        5 2007-04-13 2007-04-15
 9        9        2000       2004      2004        5 2008-02-22 2008-02-26
10       10        2199       2201      2204        6 2008-09-08 2008-09-10
# ℹ 54 more rows
# ℹ 15 more variables: date_end <date>, intensity_mean <dbl>,
#   intensity_max <dbl>, intensity_var <dbl>, intensity_cumulative <dbl>,
#   intensity_mean_relThresh <dbl>, intensity_max_relThresh <dbl>,
#   intensity_var_relThresh <dbl>, intensity_cumulative_relThresh <dbl>,
#   intensity_mean_abs <dbl>, intensity_max_abs <dbl>, intensity_var_abs <dbl>,
#   intensity_cumulative_abs <dbl>, rate_onset <dbl>, rate_decline <dbl>
false_event<-heatwaveR::detect_event(tsfalse, x=date, y=temp, S=FALSE)
false_event
$climatology
# A tibble: 7,929 × 10
     doy date        temp  seas thresh   var threshCriterion durationCriterion
   <int> <date>     <dbl> <dbl>  <dbl> <dbl> <lgl>           <lgl>            
 1   247 2002-09-03  29.4  30.2   30.6 0.348 FALSE           FALSE            
 2   248 2002-09-04  29.4  30.2   30.6 0.346 FALSE           FALSE            
 3   249 2002-09-05  29.8  30.2   30.6 0.345 FALSE           FALSE            
 4   250 2002-09-06  29.9  30.2   30.6 0.346 FALSE           FALSE            
 5   251 2002-09-07  30.1  30.2   30.6 0.347 FALSE           FALSE            
 6   252 2002-09-08  30.4  30.2   30.6 0.350 FALSE           FALSE            
 7   253 2002-09-09  30.5  30.2   30.6 0.354 FALSE           FALSE            
 8   254 2002-09-10  30.4  30.2   30.6 0.358 FALSE           FALSE            
 9   255 2002-09-11  30.2  30.2   30.6 0.362 FALSE           FALSE            
10   256 2002-09-12  29.9  30.2   30.6 0.367 FALSE           FALSE            
# ℹ 7,919 more rows
# ℹ 2 more variables: event <lgl>, event_no <int>

$event
# A tibble: 39 × 22
   event_no index_start index_peak index_end duration date_start date_peak 
      <int>       <int>      <int>     <int>    <dbl> <date>     <date>    
 1        1        1011       1013      1015        5 2005-06-09 2005-06-11
 2        2        1039       1043      1055       17 2005-07-07 2005-07-11
 3        3        1091       1095      1096        6 2005-08-28 2005-09-01
 4        4        1613       1616      1617        5 2007-02-01 2007-02-04
 5        5        2199       2202      2205        7 2008-09-09 2008-09-12
 6        6        2491       2492      2495        5 2009-06-28 2009-06-29
 7        7        2567       2569      2573        7 2009-09-12 2009-09-14
 8        8        2655       2658      2662        8 2009-12-09 2009-12-12
 9        9        2796       2796      2807       12 2010-04-29 2010-04-29
10       10        2857       2858      2861        5 2010-06-29 2010-06-30
# ℹ 29 more rows
# ℹ 15 more variables: date_end <date>, intensity_mean <dbl>,
#   intensity_max <dbl>, intensity_var <dbl>, intensity_cumulative <dbl>,
#   intensity_mean_relThresh <dbl>, intensity_max_relThresh <dbl>,
#   intensity_var_relThresh <dbl>, intensity_cumulative_relThresh <dbl>,
#   intensity_mean_abs <dbl>, intensity_max_abs <dbl>, intensity_var_abs <dbl>,
#   intensity_cumulative_abs <dbl>, rate_onset <dbl>, rate_decline <dbl>
#MHW categories
silkcat<-category(silk_event, season="range", S=FALSE)
falsecat<-category(false_event, season="range", S=FALSE)

detect_event() returns the following:
-event_no: sequential number inficating ID and order of events
-index_start: start index of event -index end: end index of event
-duration: duration of event in days
-start, peak, and end date of events -intensity mean, max, and var (this is standard deviation): intensity is temp-seas. -intensity_cumulative = deg C * days
-Many more vars defined here
-category (I, II, III, IV– associated with moderate, strong, severe, extreme) and associated vars including % of time spent at each category level


calculate exceedance of MMMT from previous calculations (30.33 at Silk, 30.36 at False)

Exceedance using SST

silk_exc<- exceedance(silk_clim, x=date, y=temp,threshold=30.33)
silk_exc
$threshold
# A tibble: 7,930 × 7
   date        temp thresh threshCriterion durationCriterion exceedance
   <date>     <dbl>  <dbl> <lgl>           <lgl>             <lgl>     
 1 2002-09-02  29.4   30.3 FALSE           FALSE             FALSE     
 2 2002-09-03  29.5   30.3 FALSE           FALSE             FALSE     
 3 2002-09-04  29.6   30.3 FALSE           FALSE             FALSE     
 4 2002-09-05  29.8   30.3 FALSE           FALSE             FALSE     
 5 2002-09-06  29.8   30.3 FALSE           FALSE             FALSE     
 6 2002-09-07  29.9   30.3 FALSE           FALSE             FALSE     
 7 2002-09-08  30.0   30.3 FALSE           FALSE             FALSE     
 8 2002-09-09  29.8   30.3 FALSE           FALSE             FALSE     
 9 2002-09-10  29.5   30.3 FALSE           FALSE             FALSE     
10 2002-09-11  29.7   30.3 FALSE           FALSE             FALSE     
# ℹ 7,920 more rows
# ℹ 1 more variable: exceedance_no <int>

$exceedance
# A tibble: 19 × 18
   exceedance_no index_start index_peak index_end duration date_start date_peak 
           <dbl>       <dbl>      <dbl>     <dbl>    <dbl> <date>     <date>    
 1             1        1041       1042      1045        5 2005-07-08 2005-07-09
 2             2        1091       1096      1097        7 2005-08-27 2005-09-01
 3             3        2199       2201      2206        8 2008-09-08 2008-09-10
 4             4        2565       2570      2574       10 2009-09-09 2009-09-14
 5             5        2858       2859      2862        5 2010-06-29 2010-06-30
 6             6        2925       2929      2930        6 2010-09-04 2010-09-08
 7             7        2938       2939      2943        6 2010-09-17 2010-09-18
 8             8        3288       3289      3294        7 2011-09-02 2011-09-03
 9             9        3667       3670      3671        5 2012-09-15 2012-09-18
10            10        4040       4043      4044        5 2013-09-23 2013-09-26
11            11        6206       6211      6213        8 2019-08-29 2019-09-03
12            12        6230       6237      6240       11 2019-09-22 2019-09-29
13            13        6257       6259      6261        5 2019-10-19 2019-10-21
14            14        6567       6570      6574        8 2020-08-24 2020-08-27
15            15        6593       6604      6608       16 2020-09-19 2020-09-30
16            16        7585       7588      7593        9 2023-06-08 2023-06-11
17            17        7672       7675      7676        5 2023-09-03 2023-09-06
18            18        7683       7697      7697       15 2023-09-14 2023-09-28
19            19        7705       7712      7717       13 2023-10-06 2023-10-13
# ℹ 11 more variables: date_end <date>, intensity_mean <dbl>,
#   intensity_max <dbl>, intensity_var <dbl>, intensity_cumulative <dbl>,
#   intensity_mean_abs <dbl>, intensity_max_abs <dbl>, intensity_var_abs <dbl>,
#   intensity_cum_abs <dbl>, rate_onset <dbl>, rate_decline <dbl>
false_exc<-exceedance(false_clim, x=date,y=temp,threshold=30.36)
false_exc
$threshold
# A tibble: 7,929 × 7
   date        temp thresh threshCriterion durationCriterion exceedance
   <date>     <dbl>  <dbl> <lgl>           <lgl>             <lgl>     
 1 2002-09-03  29.4   30.4 FALSE           FALSE             FALSE     
 2 2002-09-04  29.4   30.4 FALSE           FALSE             FALSE     
 3 2002-09-05  29.8   30.4 FALSE           FALSE             FALSE     
 4 2002-09-06  29.9   30.4 FALSE           FALSE             FALSE     
 5 2002-09-07  30.1   30.4 FALSE           FALSE             FALSE     
 6 2002-09-08  30.4   30.4 TRUE            FALSE             FALSE     
 7 2002-09-09  30.5   30.4 TRUE            FALSE             FALSE     
 8 2002-09-10  30.4   30.4 TRUE            FALSE             FALSE     
 9 2002-09-11  30.2   30.4 FALSE           FALSE             FALSE     
10 2002-09-12  29.9   30.4 FALSE           FALSE             FALSE     
# ℹ 7,919 more rows
# ℹ 1 more variable: exceedance_no <int>

$exceedance
# A tibble: 33 × 18
   exceedance_no index_start index_peak index_end duration date_start date_peak 
           <dbl>       <dbl>      <dbl>     <dbl>    <dbl> <date>     <date>    
 1             1        1011       1013      1015        5 2005-06-09 2005-06-11
 2             2        1040       1043      1045        6 2005-07-08 2005-07-11
 3             3        1089       1096      1100       12 2005-08-26 2005-09-02
 4             4        2199       2202      2205        7 2008-09-09 2008-09-12
 5             5        2563       2569      2574       12 2009-09-08 2009-09-14
 6             6        2832       2833      2836        5 2010-06-04 2010-06-05
 7             7        2857       2858      2861        5 2010-06-29 2010-06-30
 8             8        2925       2928      2931        7 2010-09-05 2010-09-08
 9             9        2936       2938      2943        8 2010-09-16 2010-09-18
10            10        3265       3268      3270        6 2011-08-11 2011-08-14
# ℹ 23 more rows
# ℹ 11 more variables: date_end <date>, intensity_mean <dbl>,
#   intensity_max <dbl>, intensity_var <dbl>, intensity_cumulative <dbl>,
#   intensity_mean_abs <dbl>, intensity_max_abs <dbl>, intensity_var_abs <dbl>,
#   intensity_cum_abs <dbl>, rate_onset <dbl>, rate_decline <dbl>
false_exceedance<-ggplot(data = false_exc$threshold, aes(x = date)) +
  geom_flame(aes(y = temp, y2 = thresh, fill = "all"), show.legend = F) +
  geom_line(aes(y = temp, colour = "temp")) +
  geom_line(aes(y = thresh, colour = "thresh"), size = 1.0) +
  scale_colour_manual(name = "Line Colour",
                      values = c("temp" = "black", "thresh" =  "forestgreen")) +
  scale_fill_manual(name = "Event Colour", values = c("all" = "firebrick")) +
  guides(colour = guide_legend(override.aes = list(fill = NA))) +
  scale_x_date(date_labels = "%b %Y", limits = as.Date(c('2017-12-01','2023-12-31'))) +
  labs(title="False Caye DHW threshold (30.3°C) exceedance",y = "Temperature (°C)", x = NULL)+
  theme_bw()

silk_exceedance<-ggplot(data = silk_exc$threshold, aes(x = date)) +
  geom_flame(aes(y = temp, y2 = thresh, fill = "all"), show.legend = F) +
  geom_line(aes(y = temp, color='temp')) +
  geom_line(aes(y = thresh, colour = "thresh"), size = 1.0) +
  scale_colour_manual(name = "Line Colour",
                      values = c("temp" = "black", "thresh" =  "forestgreen")) +
  scale_fill_manual(name = "Event Colour", values = c("all" = "firebrick")) +
  guides(colour = guide_legend(override.aes = list(fill = NA))) +
  scale_x_date(date_labels = "%b %Y", limits = as.Date(c('2017-12-01','2023-12-31'))) +
  labs(title= "Silk Caye DHW threshold (30.3°C) exceedance", y = "Temperature (°C)", x = NULL)+
  theme_bw()
  
silk_exceedance
Warning: Removed 5711 rows containing missing values or values outside the scale range
(`geom_flame()`).
Warning: Removed 5708 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 5708 rows containing missing values or values outside the scale range
(`geom_line()`).

exc_graphs<-false_exceedance/silk_exceedance+
  plot_layout(guides='collect')
exc_graphs
Warning: Removed 5710 rows containing missing values or values outside the scale range
(`geom_flame()`).
Warning: Removed 5707 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 5707 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 5711 rows containing missing values or values outside the scale range
(`geom_flame()`).
Warning: Removed 5708 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 5708 rows containing missing values or values outside the scale range
(`geom_line()`).

Exceedance using in-situ data

silk_is_exc<- exceedance(silk_isclim, x=date, y=temp,threshold=30.33)
silk_is_exc
$threshold
# A tibble: 3,342 × 7
   date        temp thresh threshCriterion durationCriterion exceedance
   <date>     <dbl>  <dbl> <lgl>           <lgl>             <lgl>     
 1 2015-04-18  29.3   30.3 FALSE           FALSE             FALSE     
 2 2015-04-19  29.0   30.3 FALSE           FALSE             FALSE     
 3 2015-04-20  29.1   30.3 FALSE           FALSE             FALSE     
 4 2015-04-21  29.3   30.3 FALSE           FALSE             FALSE     
 5 2015-04-22  29.5   30.3 FALSE           FALSE             FALSE     
 6 2015-04-23  29.5   30.3 FALSE           FALSE             FALSE     
 7 2015-04-24  29.4   30.3 FALSE           FALSE             FALSE     
 8 2015-04-25  29.4   30.3 FALSE           FALSE             FALSE     
 9 2015-04-26  29.5   30.3 FALSE           FALSE             FALSE     
10 2015-04-27  29.6   30.3 FALSE           FALSE             FALSE     
# ℹ 3,332 more rows
# ℹ 1 more variable: exceedance_no <int>

$exceedance
# A tibble: 23 × 18
   exceedance_no index_start index_peak index_end duration date_start date_peak 
           <dbl>       <dbl>      <dbl>     <dbl>    <dbl> <date>     <date>    
 1             1         130        132       141       12 2015-08-25 2015-08-27
 2             2         149        155       156        8 2015-09-13 2015-09-19
 3             3         167        175       179       13 2015-10-01 2015-10-09
 4             4         194        198       199        6 2015-10-28 2015-11-01
 5             5         410        420       421       12 2016-05-31 2016-06-10
 6             6         501        505       507        7 2016-08-30 2016-09-03
 7             7         513        515       519        7 2016-09-11 2016-09-13
 8             8         537        540       542        6 2016-10-05 2016-10-08
 9             9         860        866       869       10 2017-08-24 2017-08-30
10            10         873        896       900       28 2017-09-06 2017-09-29
# ℹ 13 more rows
# ℹ 11 more variables: date_end <date>, intensity_mean <dbl>,
#   intensity_max <dbl>, intensity_var <dbl>, intensity_cumulative <dbl>,
#   intensity_mean_abs <dbl>, intensity_max_abs <dbl>, intensity_var_abs <dbl>,
#   intensity_cum_abs <dbl>, rate_onset <dbl>, rate_decline <dbl>
false_is_exc<-exceedance(false_isclim, x=date,y=temp,threshold=30.36)
false_is_exc
$threshold
# A tibble: 1,942 × 7
   date        temp thresh threshCriterion durationCriterion exceedance
   <date>     <dbl>  <dbl> <lgl>           <lgl>             <lgl>     
 1 2018-12-05  28.8   30.4 FALSE           FALSE             FALSE     
 2 2018-12-06  28.6   30.4 FALSE           FALSE             FALSE     
 3 2018-12-07  28.5   30.4 FALSE           FALSE             FALSE     
 4 2018-12-08  28.5   30.4 FALSE           FALSE             FALSE     
 5 2018-12-09  28.6   30.4 FALSE           FALSE             FALSE     
 6 2018-12-10  28.3   30.4 FALSE           FALSE             FALSE     
 7 2018-12-11  27.6   30.4 FALSE           FALSE             FALSE     
 8 2018-12-12  27.1   30.4 FALSE           FALSE             FALSE     
 9 2018-12-13  27.5   30.4 FALSE           FALSE             FALSE     
10 2018-12-14  27.8   30.4 FALSE           FALSE             FALSE     
# ℹ 1,932 more rows
# ℹ 1 more variable: exceedance_no <int>

$exceedance
# A tibble: 20 × 18
   exceedance_no index_start index_peak index_end duration date_start date_peak 
           <dbl>       <dbl>      <dbl>     <dbl>    <dbl> <date>     <date>    
 1             1         159        164       169       11 2019-05-12 2019-05-17
 2             2         176        188       203       28 2019-05-29 2019-06-10
 3             3         219        221       223        5 2019-07-11 2019-07-13
 4             4         240        274       323       84 2019-08-01 2019-09-04
 5             5         503        508       510        8 2020-04-20 2020-04-25
 6             6         525        553       560       36 2020-05-12 2020-06-09
 7             7         583        584       587        5 2020-07-09 2020-07-10
 8             8         592        597       614       23 2020-07-18 2020-07-23
 9             9         625        663       672       48 2020-08-20 2020-09-27
10            10         913        917       918        6 2021-06-04 2021-06-08
11            11         923        931       949       27 2021-06-14 2021-06-22
12            12         964        966       969        6 2021-07-25 2021-07-27
13            13         982        988      1007       26 2021-08-12 2021-08-18
14            14        1014       1018      1050       37 2021-09-13 2021-09-17
15            15        1251       1264      1265       15 2022-05-08 2022-05-21
16            16        1280       1282      1284        5 2022-06-06 2022-06-08
17            17        1357       1359      1362        6 2022-08-22 2022-08-24
18            18        1372       1374      1382       11 2022-09-06 2022-09-08
19            19        1626       1651      1690       65 2023-05-18 2023-06-12
20            20        1696       1774      1780       85 2023-07-27 2023-10-13
# ℹ 11 more variables: date_end <date>, intensity_mean <dbl>,
#   intensity_max <dbl>, intensity_var <dbl>, intensity_cumulative <dbl>,
#   intensity_mean_abs <dbl>, intensity_max_abs <dbl>, intensity_var_abs <dbl>,
#   intensity_cum_abs <dbl>, rate_onset <dbl>, rate_decline <dbl>
false_is_exceedance<-ggplot(data = false_is_exc$threshold, aes(x = date)) +
  geom_flame(aes(y = temp, y2 = thresh, fill = "all"), show.legend = F) +
  geom_line(aes(y = temp, colour = "temp")) +
  geom_line(aes(y = thresh, colour = "thresh"), size = 1.0) +
  scale_colour_manual(name = "Line Colour",
                      values = c("temp" = "black", "thresh" =  "forestgreen")) +
  scale_fill_manual(name = "Event Colour", values = c("all" = "firebrick")) +
  guides(colour = guide_legend(override.aes = list(fill = NA))) +
  scale_x_date(date_labels = "%b %Y", limits = as.Date(c('2017-12-01','2023-12-31'))) +
  labs(title="False Caye DHW threshold (30.3°C) exceedance, in-situ temperature",y = "Temperature (°C)", x = NULL)+
  theme_bw()

silk_is_exceedance<-ggplot(data = silk_is_exc$threshold, aes(x = date)) +
  geom_flame(aes(y = temp, y2 = thresh, fill = "all"), show.legend = F) +
  geom_line(aes(y = temp, color='temp')) +
  geom_line(aes(y = thresh, colour = "thresh"), size = 1.0) +
  scale_colour_manual(name = "Line Colour",
                      values = c("temp" = "black", "thresh" =  "forestgreen")) +
  scale_fill_manual(name = "Event Colour", values = c("all" = "firebrick")) +
  guides(colour = guide_legend(override.aes = list(fill = NA))) +
  scale_x_date(date_labels = "%b %Y", limits = as.Date(c('2017-12-01','2023-12-31'))) +
  labs(title= "Silk Caye DHW threshold (30.3°C) exceedance, in-situ temperature", y = "Temperature (°C)", x = NULL)+
  theme_bw()
  
silk_is_exceedance
Warning: Removed 1259 rows containing missing values or values outside the scale range
(`geom_flame()`).
Warning: Removed 1120 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 1120 rows containing missing values or values outside the scale range
(`geom_line()`).

is_exc_graphs<-false_is_exceedance/silk_is_exceedance+
  plot_layout(guides='collect')
is_exc_graphs
Warning: Removed 89 rows containing missing values or values outside the scale range
(`geom_flame()`).
Warning: Removed 89 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 89 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 1259 rows containing missing values or values outside the scale range
(`geom_flame()`).
Warning: Removed 1120 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 1120 rows containing missing values or values outside the scale range
(`geom_line()`).

joinging the two MHW datasets (false and silk)

silk_event$climatology$site<-'Silk Caye'
silk_event$event$site<-'Silk Caye'

false_event$climatology$site<-'False Caye'
false_event$event$site<-'False Caye'

#rbind each df
satclim<-rbind(silk_event$climatology, false_event$climatology)
satevent<-rbind(silk_event$event, false_event$event)

Some visualizations

#lollipops
ggplot(satevent, aes(x = date_start, y = intensity_max)) +
  geom_lolli(colour = "firebrick", colour_n = "red", n = 3) +
  labs(y = expression(paste("Max. intensity [", degree, "C]")), x = NULL)+
  theme_classic()+
  facet_wrap(~site, ncol=1)+
    scale_x_date(date_labels = "%b %Y", limits = as.Date(c('2017-12-01','2023-12-31')))
Warning: Removed 59 rows containing missing values or values outside the scale range
(`geom_lolli()`).
Warning: Using the `size` aesthetic with geom_segment was deprecated in ggplot2 3.4.0.
ℹ Please use the `linewidth` aesthetic instead.

satclim$dhw_thresh <- ifelse(satclim$site=='Silk Caye', 30.33, 30.36)

#climatologies
all_clim<-ggplot(data = satclim, aes(x = date)) +
  geom_line(aes(y = temp), size=0.1, color='darkgray') +
  geom_line(aes(y = thresh), size = 1.0, alpha=0.5, color='black') +
  geom_line(aes(y = seas), size = 1.2, alpha=0.7,color='forestgreen') +
  geom_flame(aes(y = temp, y2 = thresh), color='firebrick',fill=NA,size=1) +
  geom_hline(aes(yintercept=dhw_thresh), color='black', linetype='dashed')+
  theme_classic()+
  labs(title="MHW Climatology")+
  facet_wrap(~site, ncol=1)
all_clim
Warning: Removed 412 rows containing missing values or values outside the scale range
(`geom_flame()`).

#just exp time
exp_clim<-ggplot(data = satclim, aes(x = date)) +
  geom_line(aes(y = temp), size=0.1, color='darkgray') +
  geom_line(aes(y = thresh), size = 1.0, alpha=0.5, color='black') +
  geom_line(aes(y = seas), size = 1.2, alpha=0.7,color='forestgreen') +
  geom_flame(aes(y = temp, y2 = thresh), color='firebrick',fill=NA,size=1) +
  geom_hline(aes(yintercept=dhw_thresh), color='black', linetype='dashed')+
  theme_classic()+
  labs(title="MHW Climatology")+
  facet_wrap(~site, ncol=1)+
  scale_x_date(date_labels = "%b %Y", limits = as.Date(c('2017-12-01','2023-12-31')))
exp_clim
Warning: Removed 5708 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 5708 rows containing missing values or values outside the scale range
(`geom_line()`).
Removed 5708 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 11421 rows containing missing values or values outside the scale range
(`geom_flame()`).

calling some useful graphs

Satellite (MUR SST) temperature compared to in-situ temperature (HOBO Pro V2)

Difference between in-situ temperature (HOBO Pro V2) and SST (JPL MUR SST) over time

Table of Max DHW by year of experiment

site year dhwm
False Caye 2018 0.000000
False Caye 2019 11.256848
False Caye 2020 18.358924
False Caye 2021 9.566876
False Caye 2022 4.597010
False Caye 2023 24.806829
False Caye 2024 6.214305
Silk Caye 2017 6.214305
Silk Caye 2018 0.000000
Silk Caye 2019 6.207591
Silk Caye 2020 7.242629
Silk Caye 2021 1.228390
Silk Caye 2022 0.376581
Silk Caye 2023 9.643905
Silk Caye 2024 0.000000

Graph of DHW by year of experiment

Climatology of entire SST dataset (2002-2024)

Climatology graph for experimental time (2018-2023)

Exceedance Graph for DHW threshold (30.3 degrees C) - satellite

Exceedance Graph for DHW threshold (30.3 degrees C) - in-situ