-
Notifications
You must be signed in to change notification settings - Fork 0
/
2_SSLAvailAndCovarExtraction.Rmd
198 lines (146 loc) · 7.23 KB
/
2_SSLAvailAndCovarExtraction.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
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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
---
title: "SSLAvailAndCovarExtraction"
author: "Kelly Kapsar"
date: "10/1/2021"
output: html_document
---
```{r load libraries}
library(tidyr)
library(dplyr)
library(sf)
library(raster)
library(ggplot2)
library(scales)
library(ggmap)
library(amt)
```
```{r read in data}
# Set seed
set.seed(1015)
# Projection information for WGS84/UTM Zone 5N (EPSG:32605)
prj <- 32605
# Create study area polygon
coords <- data.frame(lat=c(56, 62, 62, 56, 56), lon=c(-155, -155, -143, -143, -155), id="study")
study <- coords %>%
st_as_sf(coords = c("lon", "lat"), crs=4326) %>%
group_by(id) %>%
summarize(geometry = st_combine(geometry)) %>%
st_cast("POLYGON") %>%
st_transform(prj)
# Read in clean, non-land ssl used locations
ssl <- readRDS("../Data_Processed/Telemetry/watersealis.rds")
ssl$used <- 1
# Create an ID column that identifies each unique individual/year/week combination
ssl$weeklyhr_id <- factor(paste0(ssl$deploy_id, substr(ssl$weekofyear,1,4), substr(ssl$weekofyear,7,8)))
# st_write(ssl, "../Data_Processed/Telemetry/USedLocs_Clean.shp")
# Read in covariates
landmask <- raster("../Data_Processed/Landmask_GEBCO.tif")
dist_500m <- raster("../Data_Processed/Dist500m.tif") %>% raster::mask(landmask, maskvalue=1)
depth <- raster("../Data_Processed/Bathymetry.tif") %>% raster::mask(landmask, maskvalue=1)
dist_land <- raster("../Data_Processed/DistLand.tif") %>% raster::mask(landmask, maskvalue=1)
slope <- raster("../Data_Processed/slope.tif") %>% raster::mask(landmask, maskvalue=1)
ship <- readRDS("../Data_Processed/AIS_AllOther.rds")
fish <- readRDS("../Data_Processed/AIS_Fishing.rds")
sst <- readRDS("../Data_Processed/sst_weekly.rds")
wind <- readRDS("../Data_Processed/wind_weekly.rds")
# Create a list of all covariate files and check resolution
raslist <- list(depth, dist_land, dist_500m, slope, ship, fish, sst, wind)
rasres <- lapply(raslist, function(x) res(x)/1000)
# Extract dynamic covariate data at used locations
### Modified from FROM AMT PACKAGE
# https://github.com/jmsigner/amt/blob/master/R/extract_covariates.R
# Also helpful: https://cran.r-project.org/web/packages/amt/vignettes/p3_rsf.html
# I modified it so that it works with weekly timescale data
# DOES NOT work with holes in the data set any more
extract_covar_var_time_custom <- function(
xy, t, covariates) {
t_covar <- raster::getZ(covariates)
t_obs <- format(as.POSIXct(t), "%G-W%V") # Convert timestamp to week
# Same week of year as lubridate::isoweek(ssl_dates$date)
wr <- sapply(t_obs, function(x) which(x == t_covar)) # Identify which slice to select
ev <- terra::extract(covariates, xy) # Extract covariate values for all time slices at all used locations
cov_val <- ev[cbind(seq_along(wr), wr)] # select only the relevant time slice for each location
return(cov_val)
}
```
```{r avail method 4: weekly kde home ranges}
# Isolate out variables of interest
ssl_simple <- ssl %>% select(weeklyhr_id, date, lon, lat)
# Turn into an amt track
trk <- amt::make_track(st_drop_geometry(ssl_simple), .x=lon, .y=lat, .t=date, id = weeklyhr_id, date=date,
crs = CRS("+init=epsg:4326"))
trk.class<-class(trk)
# Calculate time of day based on lat/lon and timestamp
trk <- trk %>% time_of_day()
class(trk) <- trk.class
#' Now, we can transform back to geographic coordinates
trk <- amt::transform_coords(trk, CRS("+init=epsg:32605"))
# Summarize sampling rate
summarize_sampling_rate(trk)
summarize_sampling_rate_many(trk, c("id"))
# Create a list of static covariates
staticcovars <- raster::stack(dist_land, dist_500m, depth, slope)
pointcounts <- ssl %>% st_drop_geometry() %>% group_by(weeklyhr_id) %>% summarize(n=n()) %>% ungroup()
ssl_hr <- data.frame()
avail_pts <- data.frame()
# Create 5 random non-land points within each weekly MCP for each used point
for(i in 1:length(unique(trk$id))){
if(i %% 10 == 0){ print(i)}
t <- trk[which(trk$id == unique(trk$id)[i]),]
npts <- pointcounts$n[which(pointcounts$weeklyhr_id == unique(trk$id)[i])]*5
hr_kde <- hr_kde(t, levels=c(0.95))
pts <- random_points(hr_kde, n = npts) %>% mutate(weeklyhr_id = unique(trk$id)[i])
landpts <- st_as_sf(pts, coords=c("x_", "y_"), crs=prj, remove = FALSE) %>%
cbind(., raster::extract(staticcovars, .)) %>%
st_drop_geometry() %>%
filter(!is.na(.data$Bathymetry))
while(length(landpts$case_) < npts){
newpts <- random_points(hr_kde, n = npts-length(landpts$case_)) %>% mutate(weeklyhr_id = unique(trk$id)[i])
newlandpts <- st_as_sf(newpts, coords=c("x_", "y_"), crs=prj, remove = FALSE) %>%
cbind(., raster::extract(staticcovars, .)) %>%
st_drop_geometry() %>%
filter(!is.na(.data$Bathymetry))
landpts <- rbind(landpts, newlandpts)
}
hr_kde <- hr_isopleths(hr_kde) %>% mutate(weeklyhr_id = unique(trk$id)[i])
ssl_hr <- rbind(ssl_hr, hr_kde)
avail_pts <- rbind(avail_pts, landpts)
}
# Save output polygons
st_write(ssl_hr, "../Data_Processed/Telemetry/Homerange_KDE_weekly_20220705.shp")
ssl_new <- ssl %>%
select(deploy_id, weeklyhr_id, date, year, weekofyear, northing, easting) %>%
cbind(., raster::extract(staticcovars, .)) %>%
rename(dist_land = DistLand, dist_500m = Dist500m) %>%
mutate(used = 1)
# Extract week of year as date from original data
ssl_dates <- ssl %>% st_drop_geometry() %>% group_by(weeklyhr_id) %>% summarize(date=min(date))
avail_pts <- left_join(avail_pts, ssl_dates, by=c("weeklyhr_id" = "weeklyhr_id"))
# Check that week of year re-calculation worked okay
# It works!
# ssl_dates$weekofyear <- lubridate::isoweek(ssl_dates$date)
# ssl_dates$weekofyear2 <- format(ssl_dates$date, "%G-W%V")
# ssl_datetest <- left_join(ssl, ssl_dates, by= "weeklyhr_id")
# length(which(ssl_datetest$weekofyear.x != ssl_datetest$weekofyear.y))
# Clean up column names in available points
avail_pts <- avail_pts %>% mutate(deploy_id = substr(avail_pts$weeklyhr_id, 1, 13),
year = lubridate::year(avail_pts$date),
weekofyear = lubridate::isoweek(avail_pts$date),
used = 0) %>%
rename(northing = y_, easting = x_, dist_land = DistLand, dist_500m = Dist500m) %>%
select(-case_)
# Merge used with available
all_pts <- avail_pts %>% st_as_sf(., coords=c("easting", "northing"), crs=prj, remove=FALSE) %>% rbind(., ssl_new)
# Rename to match other naming schemes
ssl4 <- all_pts
# Extract dynamic weekly covariate values at used and available locations
# Extract covariate data at used locations
ssl4$sst <- extract_covar_var_time_custom(xy= st_coordinates(ssl4), t = ssl4$date, covariates=sst)
ssl4$wind <- extract_covar_var_time_custom(xy= st_coordinates(ssl4), t = ssl4$date, covariates=wind)
ssl4$ship <- extract_covar_var_time_custom(xy= st_coordinates(ssl4), t = ssl4$date, covariates=ship)
ssl4$fish <- extract_covar_var_time_custom(xy= st_coordinates(ssl4), t = ssl4$date, covariates=fish)
# convert column names to lowercase
colnames(ssl4) <- tolower(colnames(ssl4))
save(ssl4, file="../Data_Processed/Telemetry/TEMP_20220706.rda")
# browseURL("https://www.youtube.com/watch?v=K1b8AhIsSYQ&list=RDK1b8AhIsSYQ&start_radio=1")
```