-
Notifications
You must be signed in to change notification settings - Fork 0
/
heatwave.Rmd
81 lines (61 loc) · 2.78 KB
/
heatwave.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
---
title: "Marine heatwaves with heatwaveR"
output: html_notebook
---
This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code.
Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Cmd+Option+I*.
Basic detection and visualization of events
```{r}
library(dplyr)
library(ggplot2)
library(heatwaveR)
## read data: the csv file is generated by GET_SATDATA.m
# data = read.csv(file="sst_PATAG.csv")
data = read.csv(file="sst_CR.csv")
t = as.Date(data[ ,1])
temp = data[ , 2]
sst_df = data.frame(t, temp)
# Detect the events in a time series
# ts <- ts2clm(sst_df, climatologyPeriod = c("2002-07-27", "2011-07-26"))
ts <- ts2clm(sst_df, climatologyPeriod = c("2003-01-01", "2012-12-31"))
mhw <- detect_event(ts)
# View just a few metrics
mhw$event %>%
dplyr::ungroup() %>%
dplyr::select(event_no, duration, date_start, date_peak, intensity_max, intensity_cumulative) %>%
dplyr::arrange(-intensity_max) %>%
head(5)
event_line(mhw, spread = 180, metric = "intensity_max",
start_date = "2003-01-01", end_date = "2023-07-18") # start_date = "2002-07-27", end_date = "2019-12-31"
# event_line(mhw, spread = 180, start_date = "2002-07-27", end_date = "2019-12-31", category = TRUE)
event_line(mhw, spread = 180, start_date = "2003-01-01", end_date = "2023-07-18", category = TRUE)
lolli_plot(mhw, metric = "intensity_max")
ggplot(mhw$event, aes(x = date_start, y = intensity_max)) +
geom_lolli(colour = "salmon", colour_n = "red", n = 3) +
geom_text(colour = "black", aes(x = as.Date("2003-01-01"), y = 5,
label = "The marine heatwaves\nTend to be left skewed in a\nGiven time series")) +
labs(y = expression(paste("Max. intensity [", degree, "C]")), x = NULL)
```
Calculating and Visualising Exceedances
```{r}
# Calculate exceedence
exc_val <- exceedance(sst_df, threshold = 30)
# Look at a few metrics
exc_val$exceedance %>%
ungroup() %>%
select(exceedance_no, duration, date_start, date_peak, intensity_mean, intensity_max, intensity_cumulative) %>%
dplyr::arrange(-intensity_cumulative) %>%
head(5)
exc_val_thresh <- exc_val$threshold %>%
slice(5300:7000)
ggplot(data = exc_val_thresh, aes(x = t)) +
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" = "salmon")) +
guides(colour = guide_legend(override.aes = list(fill = NA))) +
scale_x_date(date_labels = "%b %Y") +
labs(y = expression(paste("Temperature [", degree, "C]")), x = NULL)
```