-
Notifications
You must be signed in to change notification settings - Fork 0
/
temperature_mackerel.R
126 lines (97 loc) · 3.98 KB
/
temperature_mackerel.R
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
`%>%` <- magrittr::`%>%`
# shapefile setup
crs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
# equal area crs
new_crs <- "+proj=utm +zone=12 +datum=NAD83 +no_defs +ellps=GRS80"
# different geom fix
sf::sf_use_s2(FALSE)
mab <- NEesp::shape %>%
dplyr::select(STRATA, geometry) %>%
sf::st_transform(proj4string = new_crs) %>%
# try geom fix?
# dplyr::mutate(geometry = geometry %>%
# s2::s2_rebuild() %>%
# sf::st_as_sfc()) %>%
dplyr::summarise(geometry = sf::st_union(geometry)) %>%
sf::st_crop(y = c(xmin = -80, xmax = -69,
ymax = 41.5, ymin = 35.8327))
# gulf of st lawrence shapefile downloaded from:
# https://www.marineregions.org/gazetteer.php?p=details&id=4290
gosl <- sf::st_read(here::here("data-raw/MackerelShapefiles/MackerelNAFO.shp"))
# southern GOSL
ggplot2::ggplot() +
ggplot2::geom_sf(data = gosl %>%
dplyr::filter(Label == "4T"),
ggplot2::aes(fill = Label)) +
viridis::scale_fill_viridis(discrete = TRUE) +
ggplot2::theme_minimal()
gosl <- gosl %>%
dplyr::filter(Label == "4T")
years <- 1982:2021
prop_spring <- c()
area_days_spring <- c()
for(j in years) {
message(paste("starting", j))
# download data ----
dir.create(here::here("data-raw","gridded", "sst_data"), recursive = TRUE)
url <- paste0("https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.", j, ".v2.nc")
download.file(url, destfile = "test.nc")
#
# R can't open the file (will have to do this in a gh action...)
# download file manually for testing on desktop
# name <- paste0(j, ".nc")
name <- "test.nc"
data <- ecopull::nc_to_raster(nc = name, varname = 'sst') # converts to NAD83
data <- raster::rotate(data)
message("converted to raster...")
# make sure all days are there ----
if(raster::nlayers(data) < 365) {
message(j, " does not have a full year of data! skipping!")
} else {
# crop to MAB ----
# filter to just may and june
# months <- c(paste0("X", j, ".05.0", 1:9),
# paste0("X", j, ".05.", 10:31),
# paste0("X", j, ".06.0", 1:9),
# paste0("X", j, ".06.", 10:30))
# filter to just june and july
months <- c(paste0("X", j, ".07.0", 1:9),
paste0("X", j, ".07.", 10:31),
paste0("X", j, ".06.0", 1:9),
paste0("X", j, ".06.", 10:30))
# mab_temp <- raster::mask(x = data[[months]],
# mask = mab)
# message("cropped to MAB...")
mab_temp <- raster::mask(x = data[[months]],
mask = gosl)
message("cropped to Gulf of St Lawrence...")
# create dataframes ----
rast_mab_df <- raster::as.data.frame(mab_temp, xy = TRUE)
message("created data frames...")
# calculate proportion between 12-16C by day ----
# calculate total area ----
total_area <- raster::area(mab_temp, na.rm = TRUE) %>%
raster::as.data.frame(xy = TRUE) %>%
dplyr::select(-c(x, y)) %>%
colSums(na.rm = TRUE)
# calculate area between 12-16C ----
opt_mab_temp <- mab_temp
opt_mab_temp@data@values[which(mab_temp@data@values > 16)] <- NA
opt_mab_temp@data@values[which(mab_temp@data@values < 12)] <- NA
optimal_area <- raster::area(opt_mab_temp, na.rm = TRUE) %>%
raster::as.data.frame(xy = TRUE) %>%
dplyr::select(-c(x, y)) %>%
colSums(na.rm = TRUE)
this_dat <- tibble::tibble(Year = j,
mean_total_area = mean(total_area),
mean_optimal_area = mean(optimal_area),
sum_optimal_area = sum(optimal_area),
proportion_optimal = mean_optimal_area/mean_total_area)
message("calculated indicators...")
area_days_spring <- rbind(area_days_spring, this_dat)
}
message(paste("done with", j))
}
print(area_days_spring)
out_data <- tibble::tibble(area_days_spring)
write.csv(out_data, here::here("data-raw/temperature_mackerel_gosl.csv"))