-
Notifications
You must be signed in to change notification settings - Fork 0
/
09_setting_migration_dates.R
122 lines (92 loc) · 3.21 KB
/
09_setting_migration_dates.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
## Indicating where geese are in migration
## 2022-08-24
## Author: S. Cunningham
library(tidyverse)
library(sf)
## Read in staging areas shapefiles
# North America - Prairie Pothole Region
ppr <- st_read("data/shapefiles", "gmannppr")
ppr <- st_transform(ppr, "+proj=aeqd +lat_0=55 +lon_0=-60 +datum=WGS84")
ppr$NAME_0 <- "PPR"
ppr <- ppr[,"NAME_0"]
# Greenland - Iceland
isl <- st_read("data/shapefiles", "ISL_adm0")
isl <- st_transform(isl, "+proj=aeqd +lat_0=55 +lon_0=-60 +datum=WGS84")
isl <- isl[,"NAME_0"]
## Read in wintering areas
# North America - Texas
texas <- st_read("data/shapefiles", "NAMCwinter_ne10m")
texas <- st_transform(texas, "+proj=aeqd +lat_0=55 +lon_0=-60 +datum=WGS84")
texas$NAME_0 <- "Texas"
texas <- texas[,"NAME_0"]
# Greenland - Ireland/Scotland
irescot <- st_read("data/shapefiles", "IrelandUK_ne10m")
irescot <- st_transform(irescot, "+proj=aeqd +lat_0=55 +lon_0=-60 +datum=WGS84")
irescot <- irescot[,"ADMIN"]
names(irescot)[1] <- "NAME_0"
## Read goose location data (one point per day
dat <- read.csv("files_for_models/daily_odba_behavior.csv")
dat <- dat[,c(1,3:6,10)]
# List of unique bird IDs
un.id <- unique(dat$animal_id)
# Set up empty data frame
mig_dates <- data.frame()
# Loop over individuals
for (i in 1:length(un.id)) {
# Subset by individual
temp <- dat[dat$animal_id==un.id[i],]
# Determine which population
pop <- temp$pop[1]
# Convert from data frame to sf
projcrs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"
temp_sf <- st_as_sf(x=temp, coords=c("longitude", "latitude"), crs=projcrs)
temp_sf <- st_transform(temp_sf, "+proj=aeqd +lat_0=55 +lon_0=-60 +datum=WGS84")
plot(st_geometry(temp_sf))
## Staging Area
if (pop=="GRLD") {
# Spatial join
pts_stage <- st_join(temp_sf, isl, join=st_within)
} else {
# Spatial join
pts_stage <- st_join(temp_sf, ppr, join=st_within)
}
# Convert to data frame
pts_stage <- st_drop_geometry(pts_stage)
pts_stage <- pts_stage[!is.na(pts_stage$NAME_0),]
# Extract date last at staging area, add 14
end <- max(pts_stage$julian) + 14
## Wintering areas
if (pop=="GRLD") {
# Spatial join
pts_stage <- st_join(temp_sf, irescot, join=st_within)
} else {
# Spatial join
pts_stage <- st_join(temp_sf, texas, join=st_within)
}
# Convert to data frame
pts_stage <- st_drop_geometry(pts_stage)
pts_stage <- pts_stage[!is.na(pts_stage$NAME_0),]
# What is last day on wintering areas?
lwint <- max(pts_stage$julian)
# what is earliest day on wintering areas?
ewint <- min(pts_stage$julian)
# how long on winter areas?
wint <- lwint - ewint
if (wint > 14) {
pts_stage <- pts_stage[pts_stage$julian>=(lwint-14),]
print(un.id[i])
}
# start date
start <- min(pts_stage$julian)
## Combine into data frame
new <- data.frame(animal_id=temp$animal_id[1], start=start, end=end)
# add row to existing data
mig_dates <- rbind(mig_dates, new)
}
mig_dates$start[mig_dates$animal_id=="RP01F"] <- 59
# Calculate duration
mig_dates <- mutate(mig_dates, duration=end-start)
# Add one to be inclusive of start and end dates
mig_dates$duration <- mig_dates$duration + 1
# Write to file
write_csv(mig_dates, "files_for_models/migration_dates.csv")