-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathextract_reefs.Rd
147 lines (106 loc) · 4.93 KB
/
extract_reefs.Rd
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/extract_reefs.R
\name{extract_reefs}
\alias{extract_reefs}
\title{Extract reef}
\usage{
extract_reefs(
input,
shpfile,
output = "sf",
fun = "mean",
weights = "area",
timeseries = "daily",
varname = "sst"
)
}
\arguments{
\item{input}{Input raster data.}
\item{shpfile}{Location of shapefile mask.}
\item{output}{Output format, "sf" (default) or "df".}
\item{fun}{see ?exact_extract for details}
\item{weights}{see ?exact_extract for details}
\item{timeseries}{time series}
\item{varname}{varname}
}
\value{
Shapefile \code{sf} object with SST details or \code{data.frame} (see above for details).
}
\description{
Function to extract SST data for shapefile overlay.
Processing large daily datasets output in sf may result in significant lag due to tidyverse
pivot_longer approach (i.e. 365 days results in an sf with 1357104 rows). Alternative approach
implemented earlier (looping over time periods) is still slow - to be refined later.
See vignette for further details
}
\examples{
\dontrun{
#### annual data
ecmwfr_combined <- process_ERA5("/Users/rof011/GBR-dhw/datasets/era5/")
gbr_era5_annual_max <- summarise_raster(ecmwfr_combined, index = "years", fun = max, na.rm = TRUE)
gbr_era5_sf_annual_max <- extract_reefs(gbr_era5_annual_max, gbr_files, output="sf", weights="area", fun="mean")
### daily data for 2024 only
ecmwfr_2024 <- ecmwfr_combined[[format(terra::time(ecmwfr_combined), "\%Y") == "2024"]]
gbr_era5_sf_2024 <- extract_reefs(ecmwfr_2024, gbr_files, output="sf", weights="area", fun="mean")
### ERA5 climatology
ecmwfr_combined <- process_ERA5("/Users/rof011/GBR-dhw/datasets/era5/")
ecmwfr_climatology <- create_climatology(ecmwfr_combined, baa=FALSE)
###### DHW
### maxDHW per year for era5 in rast
gbr_era5_dhw_annual_max <- summarise_raster(ecmwfr_climatology$dhw, index = "years", fun = max, na.rm = TRUE)
ggplot() + theme_bw() +
facet_wrap(~lyr, ncol=10) +
tidyterra::geom_spatraster(data=gbr_era5_dhw_annual_max) +
scale_fill_distiller(palette="Reds", direction=1)
### annual mean reef-average maxDHW
gbr_era5_dhw_mean_annual_max <- extract_reefs(gbr_era5_dhw_annual_max, gbr_files, output="df", varname="dhw", weights="area", fun="mean")
gbr_era5_dhw_mean_annual_max_summarised <- gbr_era5_dhw_mean_annual_max |>
group_by(date) |>
summarise(dhw = mean(dhw, na.rm=TRUE))
ggplot() + theme_bw() +
geom_vline(xintercept=c("1998", "2002", "2016", "2017", "2020", "2022", "2024"), linewidth=2, alpha=0.2) +
geom_col(data=gbr_era5_dhw_mean_annual_max_summarised, aes(x=date, y=dhw, fill=dhw), color="black", show.legend=FALSE) +
scale_fill_distiller(palette="Reds", direction=1)
###### SST
### maxSST per year for era5 in rast
gbr_era5_sst_annual_max <- summarise_raster(ecmwfr_climatology$sst, index = "years", fun = max, na.rm = TRUE)
ggplot() + theme_bw() +
facet_wrap(~lyr, ncol=10) +
tidyterra::geom_spatraster(data=gbr_era5_sst_annual_max) +
scale_fill_distiller(palette="RdBu")
### annual mean reef-average SST
gbr_era5_SST_mean_annual_max <- extract_reefs(gbr_era5_sst_annual_max, gbr_files, output="df", varname="sst", weights="area", fun="mean")
gbr_era5_sst_mean_annual_max_summarised <- gbr_era5_SST_mean_annual_max |>
group_by(date) |>
summarise(sst = mean(sst, na.rm=TRUE)) |>
mutate(date=as.numeric(date))
ggplot() + theme_bw() +
geom_vline(xintercept=c(1998, 2002, 2016, 2017, 2020, 2022, 2024), linewidth=2, alpha=0.2) +
geom_point(data=gbr_era5_sst_mean_annual_max_summarised, aes(x=date, y=sst, fill=sst),
color="black", show.legend=FALSE, shape=21) +
geom_smooth(data=gbr_era5_sst_mean_annual_max_summarised, aes(x=date, y=sst), method = "lm") +
scale_fill_distiller(palette="RdBu", direction=-1) +
coord_cartesian(ylim = c(28,30))
###### SSTanom
### meanSSTanom per year for era5 in rast
gbr_era5_sstanom_annual_max <- summarise_raster(ecmwfr_climatology$anomaly, index = "years", fun = mean, na.rm = TRUE)
ggplot() + theme_bw() +
facet_wrap(~lyr, ncol=10) +
tidyterra::geom_spatraster(data=gbr_era5_sstanom_annual_max) +
scale_fill_distiller(palette="RdBu")
### annual mean reef-average SST
gbr_era5_SSTanom_mean_annual_max <- extract_reefs(gbr_era5_sstanom_annual_max, gbr_files, output="df", varname="sstanom", weights="area", fun="mean")
gbr_era5_sstanom_mean_annual_max_summarised <- gbr_era5_SSTanom_mean_annual_max |>
group_by(date) |>
summarise(sstanom = mean(sstanom, na.rm=TRUE)) |>
mutate(date=as.numeric(date))
ggplot() + theme_bw() +
geom_vline(xintercept=c(1998, 2002, 2016, 2017, 2020, 2022, 2024), linewidth=2, alpha=0.2) +
geom_line(data=gbr_era5_sstanom_mean_annual_max_summarised, aes(x=date, y=sstanom),
color="black", show.legend=FALSE) +
geom_smooth(data=gbr_era5_sstanom_mean_annual_max_summarised, aes(x=date, y=sstanom), method = "lm", alpha=0.2) +
scale_fill_distiller(palette="RdBu", direction=-1) +
geom_hline(yintercept=0, lwd=0.4)
#coord_cartesian(ylim = c(28,30))
}
}