#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)Hundley et al temperature analysis
Load Packages
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")#satelliteprepare 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)
satandinsituWarning: 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')
satvsinsituDHW 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.33DHW 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.
dhwgraphWarning: 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()
silkclimgraphfalseclimgraph<- 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()
falseclimgraphfalseclimgraph+silkclimgraphdetect 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_exceedanceWarning: 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_graphsWarning: 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_exceedanceWarning: 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_graphsWarning: 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_climWarning: 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_climWarning: 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