Load libraries
library(data.table)
library(recurse)
library(tidyverse)
library(lubridate)
library(amt)
library(sp)
library(magrittr)Import dataset.
The csv imported is the original raw data. The only pre-processing is that columns not needed in this analysis have been removed.
crane<-fread("raw_data_but_trimmed_columns/raw_data.csv")
names(crane)[1:2]<-c("long", "lat")
names(crane)[3]<-"ID"
names(crane)[4]<-"Timestamp"
crane$Timestamp<-as.POSIXct(crane$Timestamp, format="%Y-%m-%d %H:%M:%OS")Get rid of missing locations
ids <- crane %>% select("long", "lat", "Timestamp") %>% complete.cases
crane<-crane %>% filter(., ids==TRUE)Two of the individuals used in this analysis had breeding territories in Ontario, though they were captured in northern Minnesota. All other cranes were captured on or near their breeding/natal territory. Because we filter each crane’s locations (to extract just breeding range locations for each year) individually using a cutoff based on net-squared displacement (NSD) with the origin being the capture location, these two cranes need to be pulled out of the dataset and treated individually.
# Remove the two Ontario cranes from the full dataset
cr.sub<-crane%>%filter(ID!="3E (Wambach adult)"&ID!="9C (Helliksen adult April 2015)")
# Pull those individuals out to calculate NSD sepearately for them.
cans<-crane%>%filter(ID=="3E (Wambach adult)"|ID=="9C (Helliksen adult April 2015)")
cans$yr<-year(cans$Timestamp)
wam<-cans%>%filter(ID=="3E (Wambach adult)")
hel<-cans%>%filter(ID=="9C (Helliksen adult April 2015)")Now I’ll calculate the net-squared displacement for these two cranes. I chose an origin based on visual inspection of the territory used each summer. Both cranes had high site fidelity to a small localized area.
# First crane
# manually create origin
hel[1,]<-NA
hel[1,'long']<-(-79.947)
hel[1,'lat']<-51.220
hel[1,'ID']<-"(Helliksen adult April 2015)"
# convert to UTM's
# first convert to spatial object
hel.sp<-SpatialPoints(hel[c('long','lat')], proj4string = CRS("+proj=longlat +ellps=WGS84"))
# then project to Albers Equal Area projection with specified meridians for boundaries
hel.sp<-spTransform(hel.sp, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "))
hel<-cbind(hel, hel.sp@coords)
colnames(hel)[6:7]<-c('utm_x', 'utm_y')
# now calculate net-squared displacement (NSD)
hel$nsd<-((hel$utm_x - hel$utm_x[1])^2 + (hel$utm_y - hel$utm_y[1])^2) #in meters squared
########################################################################
# Second crane
# manually create origin
wam[1,]<-NA
wam[1,'long']<-(-83.357)
wam[1,'lat']<-51.253
wam[1,'ID']<-"3E (Wambach adult)"
# There is an outlier to remove first
wam<-wam%>%filter(long>(-105)&lat<55) # 1 point dropped
# convert to UTM's
# first convert to spatial object
wam.sp<-SpatialPoints(wam[c('long','lat')], proj4string = CRS("+proj=longlat +ellps=WGS84"))
# then project to Albers Equal Area projection with specified meridians for boundaries
wam.sp<-spTransform(wam.sp, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "))
wam<-cbind(wam, wam.sp@coords)
colnames(wam)[6:7]<-c('utm_x', 'utm_y')
# now calculate NSD
wam$nsd<-((wam$utm_x - wam$utm_x[1])^2 + (wam$utm_y - wam$utm_y[1])^2) #in meters squared
###########################################################################
# Now combine both cranes together
# I'll take out the imported origin locations so they don't affect later nsd filtering
hel<-hel[-1,]
wam<-wam[-1,]
cans<-rbind(hel, wam)Now to calculate some metrics on the trajectories. I also add age information to reflect the age of the crane at time of capture as either adult or colt, which is a colloquial term for hatch-year cranes.
# First drop a couple outlier points that are errors
cr.sub<-cr.sub[cr.sub$long<(-80),]
cr.sub<-cr.sub%>%filter(lat<50)
#Switch to UTM's instead of lat long
cr<-cr.sub
# Create spatial object
cr.sp<-SpatialPoints(cr[c('long','lat')], proj4string = CRS("+proj=longlat +ellps=WGS84"))
# Project to Albers Equal Area projection with specified meridians for boundaries
cr.sp<-spTransform(cr.sp, CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "))
cr<-cbind(cr, cr.sp@coords)
colnames(cr)[5:6]<-c('utm_x','utm_y')
cr<-cr[,c('Timestamp', 'long', 'lat', 'ID', 'utm_x', 'utm_y')]
# Add in informaiton on which cranes were adults at time of capture (all others were hatch-year birds)
adult.list<-c("0J (Waubun adult #2)", "0K Deglman #1","1J (Santiago #3)" ,
"2C (Melby adult)" , "2M (Bagley)" , "3C (Larson adult #2)" ,
"3E (Wambach adult)" , "3K (Rasmussen adult)", "3M (Callaway adult)" ,
"4E Larson adult #1" , "4J (Santiago #1)", "4K (Waubun #1)" , "4M (Ogema adult)" ,
"5J (Deglman #3)" , "5K (Helliksen July adult 2015)", "5M (Beaulieu adult)" ,
"6E (Linden-Kologi)", "8K (Deglman #2)","8M (Stockyard adult)",
"9A (Santer)" , "9C (Helliksen adult April 2015)","9K (Santiago #2)" ,"9M (Inman adult)" )
cr<-mutate(cr, age=ifelse(ID%in%adult.list, "adult", "colt"))
cans<-mutate(cans, age=ifelse(ID%in%adult.list, "adult", "colt"))
#' Create a track class object using the amt package
trk <- mk_track(cr, .x=utm_x, .y=utm_y, .t=Timestamp, id = ID, age=age, lat=lat, long=long)
# Make the track "regular" by resampling to a 30 minute interval and dropping points if they aren't within 3 minutes of that frequency.
regtrk<- trk %>% nest(-id) %>%
mutate(regtrk = map(data, function(d) {
d %>%
track_resample(rate = minutes(30), tolerance = minutes(3)) %>%
filter_min_n_burst(min_n = 3)
})) %>% select(id, regtrk) %>% unnest()## Warning: All elements of `...` must be named.
## Did you want `data = c(x_, y_, t_, age, lat, long)`?## Warning: `cols` is now required.
## Please use `cols = c(regtrk)`#Convert back to track_xy object
regtrk1<-mk_track(regtrk, .x=x_, .y=y_, .t=t_, id = id, age=age, lat=lat, long=long)
#get additional track metrics
trk2<-regtrk1 %>% nest(-id) %>%
mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"),
dir_rel = map(data, direction_rel),
sl = map(data, step_lengths),
speed_m_s=map(data, speed),
nsd_=map(data, nsd))%>%unnest()## Warning: All elements of `...` must be named.
## Did you want `data = c(x_, y_, t_, age, lat, long)`?## Warning: `cols` is now required.
## Please use `cols = c(data, dir_abs, dir_rel, sl, speed_m_s, nsd_)`Now I’ll filter out some of the more extreme locations based on speed and step length.
trk.f<-trk2%>%filter(speed_m_s<25)
trk.f<-trk.f%>%filter(sl<45500) # sl is step length (this is the equivalent of 25m/s)Now to put everything back together
trk.f$yr<-year(trk.f$t_)
us.add<-trk.f%>%select(-dir_abs, -dir_rel, -sl, -speed_m_s)
# Get variables in the same order before combining datasets
cans<-cans[,c('ID', 'utm_x', 'utm_y', 'Timestamp', 'age', 'lat', 'long', 'nsd', 'yr')]
# Make variable names the same
colnames(us.add)<-c('ID', 'utm_x', 'utm_y', 'Timestamp', 'age', 'lat', 'long', 'nsd', 'yr')
df<-rbind(us.add, cans)I’ll remove individuals that experience mortality as hatch-year cranes during the summer they were captured
morts<-c("0A (Callaway colt 2015)", "0E (recovered from Ras colt #2)", "2K (Rasmussen colt #1)",
"4C (Star Lake/Tweeton)", "5E (Linden Kologi colt #2)", "7M (E of Mud Lake colt)")
filt<-df[!df$ID%in%morts,]
filt<-droplevels(filt)Now I’m going to filter the dataset so that the locations retained from each individual start from the first time in each calendar year that their NSD values go below the predetermined cutoff threshold of 100km and end at the last time that their NSD values go above the threshold.
nsd_db100<-data.frame(ID=NA, nsd=NA, utm_x=NA, utm_y=NA, Timestamp=NA, age=NA, lat=NA, long=NA, yr=NA, spring_cutoff=NA,fall_cutoff=NA,nsd_cutoff=NA )
nsd_threshold<-10000000000 #this corresponds to 100km from first gps point
# Compile filtered database
unique_ids<-unique(filt$ID)
for(i in 1:length(unique_ids)){
tmp<-filt[filt$ID==unique_ids[i],]
tmp_yrs<-unique(tmp$yr)
for(j in 1:length(tmp_yrs)){
tmp_yr_db<-tmp[tmp$yr==tmp_yrs[j],]
tmp_nsd_filt<-tmp_yr_db[tmp_yr_db$nsd<=nsd_threshold,]
if(nrow(tmp_nsd_filt)>0){
spring_cutoff<-as.POSIXct(min(tmp_nsd_filt$Timestamp))
fall_cutoff<-as.POSIXct(max(tmp_nsd_filt$Timestamp))
# filter by actual dates instead of nsd cutoff to allow cranes to roam further than the cutoff length
tmp_yr_db<-tmp_yr_db%>%filter(Timestamp>spring_cutoff)
tmp_yr_db<-tmp_yr_db%>%filter(Timestamp<fall_cutoff)
# add on year-specific temp info to keep track in the dataset
tmp_yr_db$spring_cutoff<-spring_cutoff
tmp_yr_db$fall_cutoff<-fall_cutoff
tmp_yr_db$nsd_cutoff<-nsd_threshold
#combine results
nsd_db100<-bind_rows(nsd_db100, tmp_yr_db)
}
}
}Next I’ll add specific information to consider the age category of each subadult during each calendar year.
nsd_db100$age_yr<-NA
nsd_db100<-nsd_db100[-1,]
nsd_db100[nsd_db100$yr==2014,]%<>%mutate(age_yr=case_when(
age=='adult'~'adult',
age=='colt'~'hatch-yr'))
# There were 4 colts from 2014. Two of them didn't make it to the next yr.
# 7A Hackensack didn't make it past the first fall before radio silence.
# 7C Helliksen was shot by a hunter in Texas the first winter.
# 2015
nsd_db100[nsd_db100$yr==2015,]%<>%mutate(age_yr=case_when(
age=='adult'~'adult',
age=='colt'&ID=="0C (Callaway colt 2014)"|ID=="6C (Rice Lake colt 2014)"~'second-yr',
age=='colt'&ID!="0C (Callaway colt 2014)"&ID!="6C (Rice Lake colt 2014)"~'hatch-yr'))
# 2016
nsd_db100[nsd_db100$yr==2016,]%<>%mutate(age_yr=case_when(
age=='adult'~'adult',
age=='colt'&ID=="0C (Callaway colt 2014)"|ID=="6C (Rice Lake colt 2014)"~'third-yr',
age=='colt'&ID!="0C (Callaway colt 2014)"&ID!="6C (Rice Lake colt 2014)"~'second-yr'))
# 2017
nsd_db100[nsd_db100$yr==2017,]%<>%mutate(age_yr=case_when(
age=='adult'~'adult',
age=='colt'&ID=="0C (Callaway colt 2014)"|ID=="6C (Rice Lake colt 2014)"~'fourth-yr',
age=='colt'&ID!="0C (Callaway colt 2014)"&ID!="6C (Rice Lake colt 2014)"~'third-yr'))
# 2018
nsd_db100[nsd_db100$yr==2018,]%<>%mutate(age_yr=case_when(
age=='adult'~'adult',
age=='colt'&ID=="0C (Callaway colt 2014)"|ID=="6C (Rice Lake colt 2014)"~'fifth-yr', #no fifth-yr nsd_db100anes left in 2018
age=='colt'&ID!="0C (Callaway colt 2014)"&ID!="6C (Rice Lake colt 2014)"~'fourth-yr'))The dataset at this point has been saved under the name “nsd100.csv” to be saved for future analysis.
Now I’ll get recursion information for each crane’s trajectory each year. Note that the radius used is 150 meters.
nsd_db100$ID<-as.factor(nsd_db100$ID)
nsd_db100.r = nsd_db100[, c("utm_x", "utm_y", 'Timestamp', 'ID')]
nsd_db100_recurse<-NULL
unique.ids<-unique(nsd_db100$ID)
# have to have verbose=F to not run out of memory
for(i in 1:length(unique.ids)){
tmpdat<-nsd_db100.r%>%filter(ID==unique.ids[i])
tmpvars<-nsd_db100%>%filter(ID==unique.ids[i])
tmprecurse<-getRecursions(tmpdat, 150)
nsd_db100_recurse<-rbind(nsd_db100_recurse, data.frame(id=tmpdat$ID,timestamp=tmpdat$Timestamp,
revisits=tmprecurse$revisits,
res_time=tmprecurse$residenceTime,
timeunits=tmprecurse$timeunits,
radius='150_meters',
age_yr=tmpvars$age_yr))} The dataset at this point has been saved under the name “nsd100_with_recursion.csv” to be saved for future analysis. Note that this csv is much smaller than nsd100 because only the recursion information for every buffered area is retained (i.e. not all locations).
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] magrittr_1.5 sp_1.3-1 amt_0.0.7
## [4] lubridate_1.7.4 forcats_0.4.0 stringr_1.4.0
## [7] dplyr_0.8.3 purrr_0.3.3 readr_1.3.1
## [10] tidyr_1.0.0 tibble_2.1.3 ggplot2_3.2.1
## [13] tidyverse_1.2.1 recurse_1.1.0 data.table_1.12.6
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 xfun_0.10 splines_3.6.0
## [4] haven_2.1.1 lattice_0.20-38 colorspace_1.4-1
## [7] vctrs_0.2.0 generics_0.0.2 htmltools_0.4.0
## [10] yaml_2.2.0 survival_2.44-1.1 rlang_0.4.0
## [13] pillar_1.4.2 glue_1.3.1 withr_2.1.2
## [16] modelr_0.1.5 readxl_1.3.1 lifecycle_0.1.0
## [19] munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0
## [22] raster_3.0-7 rvest_0.3.4 codetools_0.2-16
## [25] evaluate_0.14 knitr_1.25 broom_0.5.2
## [28] Rcpp_1.0.2 scales_1.0.0 backports_1.1.5
## [31] jsonlite_1.6 hms_0.5.1 digest_0.6.22
## [34] stringi_1.4.3 grid_3.6.0 rgdal_1.4-6
## [37] cli_1.1.0 tools_3.6.0 lazyeval_0.2.2
## [40] crayon_1.3.4 pkgconfig_2.0.3 zeallot_0.1.0
## [43] Matrix_1.2-17 xml2_1.2.2 assertthat_0.2.1
## [46] rmarkdown_1.16 httr_1.4.1 rstudioapi_0.10
## [49] R6_2.4.0 nlme_3.1-139 compiler_3.6.0