Poultry Mortality Risk due to Avian Influenza
Integrating animal movement and historical disease occurence to improve disease risk models.
Contact:
John: jmh09r@my.fsu.edu
date() ## [1] "Sun Aug 19 20:28:12 2018"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, and 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. The custom function ResidenceRaster() is a wrapper combining the raster and trip packages.
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
#Identifying areas where more than 1day was spent.
ZZ3 = ZZ2
ZZ3[ZZ3 < 1] = NA
REs.stkX = stack(ZZ2, ZZ3)
names(REs.stkX) = c("Duration", "Stopovers")
Plt.a = levelplot(REs.stkX,
margin = FALSE,
layout = c(2,1),
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 = " ")