library(devtools)
library(knitr)
library(ezknitr)
library(dplyr)
library(gdata)
library(sp)
library(lubridate)
library(rgeos)
library(rgdal)
#Clear environment
rm(list = ls())
Import processed fall16 dataset
fall.df<-read.csv("Processed Data/fall16.df.csv")
ep.fall<-fall.df[fall.df$pop=="EP",]
mcp.fall<-fall.df[fall.df$pop=="MCP",]
ep.sp<-ep.fall
mcp.sp<-mcp.fall
Create SpatialPointsDataFrame
coordinates(ep.sp)<-c("x","y")
coordinates(mcp.sp)<-c("x","y")
Define projection
proj4string(ep.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 ")
proj4string(mcp.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 ")
Buffer each layer of roost site locations
ep.buff<-gBuffer(ep.sp, width=5000) #5km
mcp.buff<-gBuffer(mcp.sp, width=5000) #5km
#5 kilometer buffer used instead of 3km to merge two adjoining areas that are indistinct
Calculate the intersection of the Eastern and Mid-Continent buffered roost site layers
both.pop<-gIntersection(ep.buff, mcp.buff)
Plot results
plot(ep.buff, col="red")
plot(mcp.buff, col="blue", add=T)
plot(both.pop, col="brown", add=T)
Convert from SpatialPolygons to SpatialPolygonsDataFrame so it can be exported as a shapefile
#Create a dataframe and display default rownames
buff.df<-data.frame(ID=1:length(both.pop))
both.pop<- SpatialPolygonsDataFrame(both.pop, buff.df)
#export as shapefile
# writeOGR(obj=both.pop, dsn="overlap_polys16", layer="fall_pop_overlap16_no_2flyway_cranes", driver="ESRI Shapefile")
Turn points layer into spatial first
fall.sp<-fall.df
coordinates(fall.sp)<-c("x", "y")
proj4string(fall.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 ")
Extract points within the overlap polygons
over.pts<-over(fall.sp, both.pop)
ind<-which(!is.na(over.pts))
tot.overlap<-fall.df[ind,] #this is a df of all the points that fall into any of the 3 overlap regions
#split overlap SpatialPolygonsDataFrame into each individual area of overlap
test<-SpatialPolygons(both.pop@polygons, proj4string = both.pop@proj4string) #from spdf to spatialPolygon
ptest<-disaggregate(test) #split apart from list of one into list of many
Create a SpatialPolygonDataframe for each area of overlap
for (i in 1:4){
poly<-ptest@polygons[[i]] #extract each polygon object
poly<-SpatialPolygons(list(poly),proj4string=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 ")) #convert to SpatialPolygon object
poly.id <- sapply(slot(poly, "polygons"), function(x) slot(x, "ID")) # pull info on ID slot
poly.df<-data.frame(ID=1:length(poly), row.names=poly.id) # use that info to create dataframe
poly<-SpatialPolygonsDataFrame(poly, poly.df) # create SpatialPolygonsDataFrame object
assign(paste("stage",i,sep=""),poly) # give unique name for each area of overlap
#export here as a shapefile if wanted
}
Extract crane observations from periods of overlap
st1<-over(fall.sp, stage1) # finds matching geometry between polygons and points
ind1<-which(!is.na(st1$ID)) # creates index for points that were in overlap polygons
st1.pts<-fall.df[ind1,] # subsets points layer for locations that were in area of overlap
st1.pts<-droplevels(st1.pts) # removes id labels from cranes that aren't present
# This gives overlap, but doesn't explicity say when cranes from diff pop's were there at the same time...
# I'll run a loop to determine which days have multiple cranes, and
# also which days have cranes from different populations present simultaneously
Create dataframe of overlap information
#Empty objects to fill in loop
dat.summary.1<-NULL
tempdat<-NULL
pops<-c("EP", "MCP")
st1alt<-st1.pts[order(st1.pts$day, st1.pts$id),] # order points by ascending day and id
st1days<-as.list(unique(st1alt$day)) # pull off days that had cranes present
days.vec<-sort(unlist(st1days)) # convert from list to vector
firstday<-days.vec[1] # pull off value for first day
lastday<-last(days.vec) # pull off value for last day
for(i in firstday:lastday){
tempdat<-st1alt[which(st1alt$day==i),] # pull off info for each day
if(nrow(tempdat)>0){
tempdat$mult<-ifelse(length(unique(tempdat$id))>1,"T","F") # assign if there were multiple cranes present on that day
ifelse(pops%in%tempdat$pop,tempdat$pop.over<-"T",tempdat$pop.over<-"F") # assign if there were cranes from different populations present on that day
dat.summary.1<-rbind(dat.summary.1, tempdat) # store results
}}
#This process can be automated to do all areas at once, but focusing on each area one at a
#time makes it easier to focus on the results of a specific area of overlap
Calculate the number of days visited by each crane
aggregate(dat.summary.1$day, by=list(id=dat.summary.1$id), FUN=function(x){length(unique(x))})
## id x
## 1 0J (Waubun adult #2) 52
## 2 2A (Helliksen July 2015 colt) 4
## 3 4K (Waubun #1) 57
## 4 4M (Ogema adult) 17
Create indicator variable of days when multiple cranes were present on a staging area
st1mult<-dat.summary.1[dat.summary.1$mult=="T",]
length(unique(st1mult$day)) #53 days
## [1] 53
Create indicator variable of days in which multiple cranes from different population were present
st1pov<-dat.summary.1[dat.summary.1$pop.over=="T",]
length(unique(st1pov$day)) #50 days
## [1] 50
Recreate previous steps for the second staging area
st2<-over(fall.sp, stage2) # finds matching geometry between polygons and points
ind2<-which(!is.na(st2$ID)) # creates index for points that were in overlap polygons
st2.pts<-fall.df[ind2,] # subsets points layer for locations that were in area of overlap
st2.pts<-droplevels(st2.pts) # removes id labels from cranes that aren't present
#Empty objects to fill in loop
dat.summary.2<-NULL
tempdat<-NULL
st2alt<-st2.pts[order(st2.pts$day, st2.pts$id),] # order points by ascending day and id
st2days<-as.list(unique(st2alt$day)) # pull off days that had cranes present
days.vec<-sort(unlist(st2days)) # convert from list to vector
firstday<-days.vec[1] # pull off value for first day
lastday<-last(days.vec) # pull off value for last day
for(i in firstday:lastday){
tempdat<-st2alt[which(st2alt$day==i),] # pull off info for each day
if(nrow(tempdat)>0){
tempdat$mult<-ifelse(length(unique(tempdat$id))>1,"T","F") # assign if there were multiple cranes present on that day
ifelse(pops%in%tempdat$pop,tempdat$pop.over<-"T",tempdat$pop.over<-"F") # assign if there were cranes from different populations present on that day
dat.summary.2<-rbind(dat.summary.2, tempdat) # store results
}}
Calculate the number of days visited by each crane
aggregate(dat.summary.2$day, by=list(id=dat.summary.2$id), FUN=function(x){length(unique(x))})
## id x
## 1 3K (Rasmussen adult) 16
## 2 5C (Stockyard colt) 24
## 3 5K (Helliksen July adult 2015) 21
## 4 6M (Larson colt #2) 19
## 5 7E (Nora Township colt) 11
Create indicator variable of days when multiple cranes were present on a staging area
st2mult<-dat.summary.2[dat.summary.2$mult=="T",]
length(unique(st2mult$day)) #28 days
## [1] 28
Create indicator variable of days in which multiple cranes from different population were present
st2pov<-dat.summary.2[dat.summary.2$pop.over=="T",]
length(unique(st2pov$day)) #21 days
## [1] 21
The third area of overlap
st3<-over(fall.sp, stage3) # finds matching geometry between polygons and points
ind3<-which(!is.na(st3$ID)) # creates index for points that were in overlap polygons
st3.pts<-fall.df[ind3,] # subsets points layer for locations that were in area of overlap
st3.pts<-droplevels(st3.pts) # removes id labels from cranes that aren't present
#Empty objects to fill in loop
dat.summary.3<-NULL
tempdat<-NULL
st3alt<-st3.pts[order(st3.pts$day, st3.pts$id),] # order points by ascending day and id
st3days<-as.list(unique(st3alt$day)) # pull off days that had cranes present
days.vec<-sort(unlist(st3days)) # convert from list to vector
firstday<-days.vec[1] # pull off value for first day
lastday<-last(days.vec) # pull off value for last day
for(i in firstday:lastday){
tempdat<-st3alt[which(st3alt$day==i),] # pull off info for each day
if(nrow(tempdat)>0){
tempdat$mult<-ifelse(length(unique(tempdat$id))>1,"T","F") # assign if there were multiple cranes present on that day
ifelse(pops%in%tempdat$pop,tempdat$pop.over<-"T",tempdat$pop.over<-"F") # assign if there were cranes from different populations present on that day
dat.summary.3<-rbind(dat.summary.3, tempdat) # store results
}
}
Calculate the number of days visited by each crane
aggregate(dat.summary.3$day, by=list(id=dat.summary.3$id), FUN=function(x){length(unique(x))})
## id x
## 1 2A (Helliksen July 2015 colt) 13
## 2 2M (Bagley) 33
## 3 3E (Wambach adult) 10
## 4 4E Larson adult #1 35
## 5 5C (Stockyard colt) 20
## 6 7E (Nora Township colt) 3
## 7 7J (Melby colt #1) 18
## 8 7K (Larson colt #1) 2
## 9 8M (Stockyard adult) 13
Create indicator variable of days when multiple cranes were present on a staging area
st3mult<-dat.summary.3[dat.summary.3$mult=="T",]
length(unique(st3mult$day)) #24 days
## [1] 40
Create indicator variable of days in which multiple cranes from different population were present
st3pov<-dat.summary.3[dat.summary.3$pop.over=="T",]
length(unique(st3pov$day)) #13 days
## [1] 13
The fourth area of overlap
st4<-over(fall.sp, stage4) # finds matching geometry between polygons and points
ind4<-which(!is.na(st4$ID)) # creates index for points that were in overlap polygons
st4.pts<-fall.df[ind4,] # subsets points layer for locations that were in area of overlap
st4.pts<-droplevels(st4.pts) # removes id labels from cranes that aren't present
#Empty objects to fill in loop
dat.summary.4<-NULL
tempdat<-NULL
st4alt<-st4.pts[order(st4.pts$day, st4.pts$id),] # order points by ascending day and id
st4days<-as.list(unique(st4alt$day)) # pull off days that had cranes present
days.vec<-sort(unlist(st4days)) # convert from list to vector
firstday<-days.vec[1] # pull off value for first day
lastday<-last(days.vec) # pull off value for last day
for(i in firstday:lastday){
tempdat<-st4alt[which(st4alt$day==i),] # pull off info for each day
if(nrow(tempdat)>0){
tempdat$mult<-ifelse(length(unique(tempdat$id))>1,"T","F") # assign if there were multiple cranes present on that day
ifelse(pops%in%tempdat$pop,tempdat$pop.over<-"T",tempdat$pop.over<-"F") # assign if there were cranes from different populations present on that day
dat.summary.4<-rbind(dat.summary.4, tempdat) # store results
}
}
Calculate the number of days visited each crane
aggregate(dat.summary.4$day, by=list(id=dat.summary.4$id), FUN=function(x){length(unique(x))})
## id x
## 1 5K (Helliksen July adult 2015) 2
## 2 7E (Nora Township colt) 3
Create indicator variable of days when multiple cranes were present on a staging area
st4mult<-dat.summary.4[dat.summary.4$mult=="T",]
length(unique(st4mult$day)) #0 days
## [1] 0
Create indicator variable of days in which multiple cranes from different population were present
st4pov<-dat.summary.4[dat.summary.4$pop.over=="T",]
length(unique(st4pov$day)) #0 days
## [1] 0
Combine all dataframes together for export
dat.summary.1$area<-"1"
dat.summary.2$area<-"2"
dat.summary.3$area<-"3"
dat.summary.4$area<-"4"
overlap.df<-rbind(dat.summary.1, dat.summary.2, dat.summary.3, dat.summary.4)
stage.summary<-aggregate(overlap.df$day, by=list(id=overlap.df$id, pop=overlap.df$pop, age=overlap.df$age, gender=overlap.df$gender),
FUN=function(x){length(unique(x))})
stage.summary[order(stage.summary$pop, stage.summary$x),]
## id pop age gender x
## 13 7K (Larson colt #1) EP colt male 2
## 3 8M (Stockyard adult) EP adult female 13
## 10 3K (Rasmussen adult) EP adult male 16
## 12 7E (Nora Township colt) EP colt male 17
## 8 7J (Melby colt #1) EP colt female 18
## 7 6M (Larson colt #2) EP colt female 19
## 1 2M (Bagley) EP adult female 33
## 2 4E Larson adult #1 EP adult female 35
## 6 5C (Stockyard colt) EP colt female 44
## 9 0J (Waubun adult #2) EP adult male 52
## 11 3E (Wambach adult) MCP adult male 10
## 5 4M (Ogema adult) MCP adult female 17
## 14 2A (Helliksen July 2015 colt) MCP colt male 17
## 15 5K (Helliksen July adult 2015) MCP colt male 22
## 4 4K (Waubun #1) MCP adult female 57
# So 14 cranes were in the overlapped staging areas, and here are the number of days each spend in any one of them.
Export for use in Event Chart script
# write.csv(x=overlap.df, file="overlap16.df.csv",row.names=FALSE)
Session Info
devtools::session_info()
## Session info -------------------------------------------------------------
## setting value
## version R version 3.4.0 (2017-04-21)
## system x86_64, mingw32
## ui RStudio (1.0.143)
## language (EN)
## collate English_United States.1252
## tz America/Denver
## date 2017-06-01
## Packages -----------------------------------------------------------------
## package * version date source
## assertthat 0.2.0 2017-04-11 CRAN (R 3.4.0)
## base * 3.4.0 2017-04-21 local
## compiler 3.4.0 2017-04-21 local
## datasets * 3.4.0 2017-04-21 local
## DBI 0.6-1 2017-04-01 CRAN (R 3.4.0)
## devtools * 1.13.1 2017-05-13 CRAN (R 3.4.0)
## digest 0.6.12 2017-01-27 CRAN (R 3.4.0)
## dplyr * 0.5.0 2016-06-24 CRAN (R 3.4.0)
## evaluate 0.10 2016-10-11 CRAN (R 3.4.0)
## ezknitr * 0.6 2016-09-16 CRAN (R 3.4.0)
## gdata * 2.17.0 2015-07-04 CRAN (R 3.4.0)
## graphics * 3.4.0 2017-04-21 local
## grDevices * 3.4.0 2017-04-21 local
## grid 3.4.0 2017-04-21 local
## gtools 3.5.0 2015-05-29 CRAN (R 3.4.0)
## highr 0.6 2016-05-09 CRAN (R 3.4.0)
## hms 0.3 2016-11-22 CRAN (R 3.4.0)
## knitr * 1.15.1 2016-11-22 CRAN (R 3.4.0)
## lattice 0.20-35 2017-03-25 CRAN (R 3.4.0)
## lazyeval 0.2.0 2016-06-12 CRAN (R 3.4.0)
## lubridate * 1.6.0 2016-09-13 CRAN (R 3.4.0)
## magrittr 1.5 2014-11-22 CRAN (R 3.4.0)
## memoise 1.1.0 2017-04-21 CRAN (R 3.4.0)
## methods * 3.4.0 2017-04-21 local
## R.methodsS3 1.7.1 2016-02-16 CRAN (R 3.4.0)
## R.oo 1.21.0 2016-11-01 CRAN (R 3.4.0)
## R.utils 2.5.0 2016-11-07 CRAN (R 3.4.0)
## R6 2.2.1 2017-05-10 CRAN (R 3.4.0)
## Rcpp 0.12.10 2017-03-19 CRAN (R 3.4.0)
## readr * 1.1.0 2017-03-22 CRAN (R 3.4.0)
## rgdal * 1.2-7 2017-04-25 CRAN (R 3.4.0)
## rgeos * 0.3-23 2017-04-06 CRAN (R 3.4.0)
## sp * 1.2-4 2016-12-22 CRAN (R 3.4.0)
## stats * 3.4.0 2017-04-21 local
## stringi 1.1.5 2017-04-07 CRAN (R 3.4.0)
## stringr 1.2.0 2017-02-18 CRAN (R 3.4.0)
## tibble 1.3.0 2017-04-01 CRAN (R 3.4.0)
## tools 3.4.0 2017-04-21 local
## utils * 3.4.0 2017-04-21 local
## withr 1.0.2 2016-06-20 CRAN (R 3.4.0)
spun with: ezknitr::ezspin(“Overlap Scripts/Fall 2016 Overlap.R”, out_dir=“output”)