Relationship of Avian Influenza Outbreaks to BWTE Residence Time
Joint modeling of duck telemetry and A.I. outbreaks across North America
Contact:
John: jmh09r@my.fsu.edu
date() ## [1] "Sun Aug 19 21:24:07 2018"library(knitr)
library(rgl)
opts_knit$set(verbose = FALSE)
knit_hooks$set(webgl = hook_webgl)suppressMessages(library(INLA))
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(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))
suppressMessages(library(changepoint))
suppressMessages(library(ecp))
suppressMessages(library(geosphere))
suppressMessages(library(moveHMM))
suppressMessages(library(adegenet))
suppressMessages(library(deldir))
suppressMessages(library(sp))
suppressMessages(library(spatstat))
suppressMessages(library(classInt))
suppressMessages(library(viridis))
suppressMessages(library(kableExtra))
source("RScripts.R")Load and Prep Data
Removing “dead” locations.
Selecting “plausible” records from the DAF_Filter.
Formating date, adding local time (America/Chicago) and getting Hour, Day, Week, season.
Argos_sp = read.csv("C:/Users/humph173/Documents/LabWork/Patuxent/Douglas_080218/BWTeal_2013-2017_Ramey_argos_tabular_diag_gis_filtered.csv",
stringsAsFactors=FALSE,
header = TRUE, sep=",")[,1:35] %>%
filter(Fate == "alive",
DAF_Filter == 0) %>%
mutate(ID = Animal_ID,
ObsDate = as.POSIXct(Location_Datetime_GMT,
tz = "GMT",
format = "%Y/%m/%d %H:%M:%S"),
ObsDateD = as.Date(Location_Datetime_GMT),
LocT = ymd_hms(format(ObsDate, tz="America/Chicago", usetz=TRUE)),
LocHR = hour(LocT),
LocDY = day(LocT),
LocWK = week(LocT),
LocMN = month(LocT))
Argos_sp$Season = ifelse(Argos_sp$LocMN == 12 | Argos_sp$LocMN == 1, "Winter",
ifelse(Argos_sp$LocMN == 3, "Spring",
ifelse(Argos_sp$LocMN >= 4 & Argos_sp$LocMN <= 8, "Breeding",
ifelse(Argos_sp$LocMN >= 9 & Argos_sp$LocMN <= 11, "Fall",
ifelse(Argos_sp$LocMN == 2 & Argos_sp$Day <= 15, "Winter",
ifelse(Argos_sp$LocMN == 2 & Argos_sp$Day >= 16, "Spring",
NA))))))Get Boundaries of North America
NAmer.reg = readOGR(dsn = "./Arc/NAmer/States",
layer = "NA_regions",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\LabWork\Patuxent\Patux_level1\Arc\NAmer\States", layer: "NA_regions"
## with 127 features
## It has 36 fields
## Integer64 fields read as strings: FID_NA_reg FID_Cuba FID_ FID_NA_r_1 FID_0 MEXST0_ MEXST0_ID TOTPOP MALE FEMALE FID_NA_sta FID_Bahama OBJECTID TOTPOP_CYLL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
NAmer = gUnaryUnion(NAmer.reg)
tmp.ras = raster(res=0.5) #raster verion
extent(tmp.ras) = extent(NAmer)
Domain.r = rasterize(NAmer,
tmp.ras,
field = 0,
background = NA)
Grd.pnt = rasterToPoints(Domain.r, sp = TRUE)
Grd.pnt@data = Grd.pnt@data %>%
mutate(Long = Grd.pnt@coords[,1],
Lat = Grd.pnt@coords[,2]) %>%
select(-layer)
GLakes = readOGR(dsn = "./Arc/RedRast/GLake",
layer = "GLakes",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Documents\LabWork\Patuxent\Patux_level1\Arc\RedRast\GLake", layer: "GLakes"
## with 5 features
## It has 6 fields
## Integer64 fields read as strings: OBJECTID_1 ObjectIDConvert to Points
Locs = cbind(Argos_sp$Longitude, Argos_sp$Latitude)
Sp.pnts = SpatialPointsDataFrame(Locs, Argos_sp)
proj4string(Sp.pnts) = proj4string(NAmer.reg)Grid Residence Time
Estimating duration of time (hours) spent in each grid cell. Individual birds are returned as layers in a raster stack. Birds are then combined to get the maximum duration any one bird spent in each cell. Time spent in grid cells lacking satellite fixes are interpolated assuming uniform movement between successful satellite captures.
Animal.IDs = levels(as.factor(Sp.pnts$Animal_ID))
for(i in 1:length(Animal.IDs)){
Ani.tmp = subset(Sp.pnts, Animal_ID == Animal.IDs[i])
tmp.res = ResidenceRaster(Ani.tmp, Domain.r)
if(i == 1){Res.stk = tmp.res}
else{Res.stk = stack(Res.stk, tmp.res)}
}
names(Res.stk) = Animal.IDs
Res.stk
ZZ = max(Res.stk)
ZZView Results
rng = seq(0, 105, 1)
mCols = brewer.pal(9, "YlOrRd")[-c(1:2)]
cr0 = colorRampPalette(mCols)(n = 2000)
cr = colorRampPalette(c(cr0),
bias = 4)
#Removing background zeros for plotting
ZZ2 = ZZ
ZZ2[ZZ2 == 0] = NA
Plt.a = levelplot(ZZ2,
margin = FALSE,
xlab =NULL,
ylab = NULL,
main = "Residence Time (days)",
maxpixels = 1e5,
col.regions = cr, at =rng,
colorkey = list(labels=list(at=c(0, 20, 40, 60, 80, 103),
labels=c("0", "20", "40", "60", "80", "100 Days"), cex=1),
labels=list(cex=12),
space = "bottom"),
par.strip.text = list(fontface='bold', cex=1),
par.settings = list(axis.line = list(col = "black"),
strip.background = list(col = 'transparent'),
strip.border = list(col = 'transparent')),
scales = list(cex = 1)) +
latticeExtra::layer(sp.polygons(NAmer.reg, col = "darkgray", lwd = 0.5))
Plt.a
Construct Triangulated Mesh
NA.bf = gBuffer(NAmer,
width = 1,
byid = FALSE)
NA.bnds = inla.sp2segment(NA.bf)
MaxEdge = 0.5
mesh0 = inla.mesh.2d(boundary = NA.bnds,
cutoff = 1,
max.edge = MaxEdge,
min.angle = 21)
MeshLocs = cbind(mesh0$loc[,1], mesh0$loc[,2])
xyz = as.data.frame(
inla.mesh.map(MeshLocs,
projection = "longlat",
inverse = TRUE))
true.radius.of.earth = 6371 #km
radius.of.earth = 1
mesh.dom = inla.mesh.create(loc = xyz,
cutoff = 35/true.radius.of.earth,
refine=list(max.edge = 1000/true.radius.of.earth,
min.angle = 26))Click, drag, and zoom to view mesh
mesh.dom$n ## [1] 9555plot(mesh.dom, rgl = TRUE, main = " ")