Argos Telemetry Data
Spatial Analysis (BWTE 2013-2017)
Contact: jmh09r@my.fsu.edu
date() ## [1] "Mon Mar 19 23:38:32 2018"Options
library(knitr)
library(rgl)
opts_knit$set(verbose = TRUE)
knit_hooks$set(webgl = hook_webgl)Set Directory
setwd("D:/Patuxent/BWTE_ARamey/BWTE_ARamey_2013-2017_Argos_data-20180307T004644Z-001/BWTE_ARamey_2013-2017_Argos_data/Argos_data")Load Packages:
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(raster))
suppressMessages(library(factoextra))
suppressMessages(library(corrplot))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(gganimate))
suppressMessages(library(ggthemes))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(gridExtra))
suppressMessages(library(dplyr))
suppressMessages(library(data.table))
suppressMessages(library(readr))
suppressMessages(library(stringr))
suppressMessages(library(lubridate))
suppressMessages(library(RcppRoll))## Warning: package 'RcppRoll' was built under R version 3.4.4suppressMessages(library(changepoint))## Warning: package 'changepoint' was built under R version 3.4.4suppressMessages(library(ecp))## Warning: package 'ecp' was built under R version 3.4.4suppressMessages(library(geosphere))
source("RScripts.R")Load Data
Second dataset inlcudes “dead dates” for each bird as derived here: http://rpubs.com/JMHumphreys/BWTE_Clean
Argos_sp = read.csv("./BWTeal_2013-2017_Ramey_argos_tabular_ds.csv",
stringsAsFactors=FALSE,
header = TRUE, sep=",")
DeadShed.df = read.csv("./DeadShed031918.csv",
stringsAsFactors=FALSE,
header = TRUE, sep=",")Remove “Dead” Locations
Also, removing locations with missing Lon/Lat solutions.
Animal.IDs = levels(as.factor(Argos_sp$Animal_ID))
for(i in 1:length(Animal.IDs)){
DeathPoint = (DeadShed.df %>%
filter(Animal == Animal.IDs[i]))[,18]
Argos_sp_cln = Argos_sp %>%
filter(Animal_ID == Animal.IDs[i]) %>%
mutate(ObsDate = ymd_hms(Message_Datetime_GMT),
ObsDateD = as.Date(Message_Datetime_GMT))
Argos_sp_cln = Argos_sp_cln %>%
mutate(Idx = 1:nrow(Argos_sp_cln)) %>%
filter(Idx <= DeathPoint,
is.na(Lon_Solution_1) == FALSE,
is.na(Lat_Solution_1) == FALSE)
if(i == 1){SpOut.df = Argos_sp_cln}
else{SpOut.df = rbind(SpOut.df, Argos_sp_cln)}
}
Argos_sp_cln = SpOut.dfCheck Accuracy Codes
Also, removing codes: 0, A, B
Code descriptions: http://www.argos-system.org/manual/3-location/34_location_classes.htm
Argos_sp_cln %>%
group_by(Location_Class) %>%
summarise(Count = length(Location_Class))## # A tibble: 7 x 2
## Location_Class Count
## <chr> <int>
## 1 0 8499
## 2 1 21325
## 3 2 36804
## 4 3 66862
## 5 A 6661
## 6 B 6547
## 7 Z 142Argos_sp_cln = Argos_sp_cln %>%
filter(Location_Class == 1 |
Location_Class == 2 |
Location_Class == 3)Get Boundaries of North America
Used to define model boundaries and to visualize results.
NAmer = readOGR(dsn = "D:/Patuxent/Patux_level1/Arc/NAmer",
layer = "NAmer",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "D:/Patuxent/Patux_level1/Arc/NAmer", layer: "NAmer"
## with 4 features
## It has 5 fields
## Integer64 fields read as strings: FID_LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
tmp.ras = raster(res=1) #raster verion
extent(tmp.ras) = extent(NAmer)
Domain.r = rasterize(NAmer,
tmp.ras,
field = 0,
background = NA)
#State boundaries
States = map("state",
fill = TRUE,
plot = FALSE)
IDs = sapply(strsplit(States$names, ":"),
function(x) x[1])
States = map2SpatialPolygons(States, IDs = IDs,
proj4string = CRS(LL84))
#Add a dataframe
pid = sapply(slot(States, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame(ID=1:length(States),
row.names = pid)
States = SpatialPolygonsDataFrame(States, p.df)Convert Argos to Spatial Points
Locs = cbind(Argos_sp_cln$Lon_Solution_1, Argos_sp_cln$Lat_Solution_1)
Sp.pnts = SpatialPointsDataFrame(Locs , Argos_sp_cln)
proj4string(Sp.pnts) = proj4string(States)Quick Plot
Arbitrarily selecting a couple birds to check points.
Obs.pnt = subset(Sp.pnts, Animal_ID == "BWTE13S_131011" | Animal_ID == "BWTE15S_135860")
wmap_df = fortify(NAmer)## Regions defined for each Polygonsstates_df = fortify(States)## Regions defined for each Polygonstmp.df = Obs.pnt@data %>%
select(Lon_Solution_1, Lat_Solution_1, Animal_ID)
ggplot(wmap_df, aes(long, lat, group=group)) +
geom_polygon(fill = "tan", col="black") +
geom_polygon(data = states_df, aes(long,lat, group=group),
fill = "tan", col="black") +
geom_point(data=tmp.df,
aes(Lon_Solution_1, Lat_Solution_1,
group = tmp.df$Animal_ID, fill=NULL),
col = "red", size=1.5) +
facet_wrap(~Animal_ID, ncol = 1) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("BWTE") +
scale_x_longitude(xmin=-175, xmax=50, step=20) +
scale_y_latitude(ymin=10, ymax=80, step=20) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_text(size=20, face="bold", hjust = 0.5))
Create Tracks
nProj = CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5
+lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=km +no_defs")
Sp.pntsp = spTransform(Sp.pnts, nProj)
BWTE.P = list(0)
TagID = sort(unique(Animal.IDs)) # Creates vector of the tag ids in alphabetical order
for(i in 1:length(TagID))
BWTE.P[[i]] = Lines(list(Line(coordinates(Sp.pntsp[Sp.pntsp@data$Animal_ID==TagID[i],]))), ID=TagID[i])
BWTE.trak = SpatialLines(BWTE.P)
proj4string(BWTE.trak) = nProj
Distances = SpatialLinesLengths(BWTE.trak, longlat=FALSE)
DistTravelled = data.frame(Animal_ID = TagID, Distances)
BWTE.trak = SpatialLinesDataFrame(BWTE.trak,
data = DistTravelled,
match.ID=FALSE)Quick Plot
Arbitrarily selecting a couple birds to check points.
Obs.trk = subset(spTransform(BWTE.trak, proj4string(States)),
Animal_ID == "BWTE13S_131011" )
Obs.trk = fortify(Obs.trk)
wmap_df = fortify(NAmer)## Regions defined for each Polygonsstates_df = fortify(States)## Regions defined for each Polygonsggplot(wmap_df, aes(long, lat, group=group)) +
geom_polygon(fill = "tan", col="black") +
geom_polygon(data = states_df, aes(long,lat, group=group),
fill = "tan", col="black") +
geom_path(data=Obs.trk,
aes(long, lat, group=group),
col = "red", size=0.75) +
xlab("Longitude") +
ylab("Latitude") +
ggtitle("BWTE13S_131011") +
scale_x_longitude(xmin=-175, xmax=50, step=20) +
scale_y_latitude(ymin=10, ymax=80, step=20) +
coord_equal() +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size=12, face="bold"),
axis.title.x = element_text(size=16, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
plot.title = element_text(size=20, face="bold", hjust = 0.5))
Total Distance
ggplot(DistTravelled, aes(Animal_ID, Distances)) +
geom_bar(colour="black", stat="identity") +
ylab("Total Distance (km)") +
theme_classic() +
guides(col = guide_legend(ncol = 16)) +
theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=12, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(face="bold", size=12, vjust=1, hjust=1, angle=90),
legend.position="bottom")
Daily Distance
Sp.pntsp$Dist = c(0, spDists(Sp.pntsp,
longlat = F,
segments = T,
diagonal = FALSE))
Obs.pnt = subset(Sp.pntsp, Animal_ID == "BWTE15S_135861" | Animal_ID == "BWTE15S_135855")
Obs.df = as.data.frame(Obs.pnt@data %>%
group_by(ObsDateD, Animal_ID) %>%
summarise(Ddist = sum(Dist))) %>%
filter(Ddist <= 1000)
ggplot(Obs.df, aes(ObsDateD, Ddist)) +
geom_bar(colour="black", stat="identity") +
ylab("Daily Distance (km)") +
theme_classic() +
facet_wrap(~Animal_ID, ncol = 1) +
guides(col = guide_legend(ncol = 16)) +
theme(panel.border = element_rect(color = "black", fill = NA, size = 1),
strip.text.x = element_text(size=12, face="bold"),
axis.title.y = element_text(size=16, face="bold"),
axis.text.x = element_text(face="bold", size=12, vjust=1, hjust=1, angle=90),
legend.position="bottom")