Weekly Risk
1 Load Boundaries
2 Example
ISSUE: Here, the raw abundance estimates for Northern Shoveler (NORS) during the Spring Season (February 16 through March 31) are loaded. Because some of the eBird observations used in estimating abundance were from the southern U.S., predicted abundances include birds just starting migration (early February) and those already at breeding grounds (late March). The map coverage is very broad, making it challenging to identify the most probable locations of bird occurrence.
NORS.spring = raster("E:/Dabbler_Pre/Dabbler_project/Raster_outputs/abun_nors_Spring.tif")
breaks = seq(0, 101, 0.1)
mCols = brewer.pal(9, "YlOrRd")
cr = colorRampPalette(mCols)(length(breaks)-1)
cr[1] = "lightgray"
cr = colorRampPalette(cr,
bias = 2.5, space = "rgb")
M0 = levelplot(NORS.spring,
margin = FALSE,
xlab = "NORS ABundance",
ylab = NULL,
maxpixels = 1e5,
col.regions = cr, at = breaks,
colorkey = list(labels=list(at=c(0, 20, 40, 60, 80, 100),
labels=c("0", "20", "40", "60", "80", "100 +"), cex=1.5),
labels=list(cex=12),
space = "right"),
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.5)) +
latticeExtra::layer(sp.polygons(States, col = "darkgray", lwd = 0.5))
M0
Load and Clean the Orginal Observation Data (eBird) used to model NORS abundance. Also, using the reported eBird observation date to determine the week of year that each observation was made.
CleanBird("E:/Dabbler_Pre/Dabbler_project/eBird_data/All_Observations/ebd_norsho_relFeb-2018.txt", "NORS", 2016, 2017)## Parsed with column specification:
## cols(
## .default = col_character(),
## `LAST EDITED DATE` = col_datetime(format = ""),
## `TAXONOMIC ORDER` = col_double(),
## `SUBSPECIES COMMON NAME` = col_logical(),
## `SUBSPECIES SCIENTIFIC NAME` = col_logical(),
## `BREEDING BIRD ATLAS CODE` = col_logical(),
## `BREEDING BIRD ATLAS CATEGORY` = col_logical(),
## COUNTY = col_logical(),
## SUBNATIONAL2_CODE = col_logical(),
## `IBA CODE` = col_logical(),
## `BCR CODE` = col_logical(),
## `USFWS CODE` = col_logical(),
## `ATLAS BLOCK` = col_logical(),
## LATITUDE = col_double(),
## LONGITUDE = col_double(),
## `OBSERVATION DATE` = col_date(format = ""),
## `TIME OBSERVATIONS STARTED` = col_time(format = ""),
## `DURATION MINUTES` = col_double(),
## `EFFORT DISTANCE KM` = col_double(),
## `EFFORT AREA HA` = col_double(),
## `NUMBER OBSERVERS` = col_double()
## # ... with 6 more columns
## )## See spec(...) for full column specifications.NORS$Week = week(NORS$Date)
head(NORS)## Spp Long Lat Count nObs dMins EDis EArea
## 1 Spatula clypeata 55.55655 25.43041 16 1 60 3.000 NA
## 2 Spatula clypeata 55.62782 24.08784 20 1 90 2.200 NA
## 3 Spatula clypeata 55.24670 25.18562 5 3 60 2.000 NA
## 4 Spatula clypeata 55.25211 24.83147 8 1 180 3.000 NA
## 5 Spatula clypeata 55.32127 24.84597 1 2 180 64.372 NA
## 6 Spatula clypeata 55.30002 24.83536 1 2 180 64.372 NA
## Date Year Month Day Season uStep Week
## 1 2016-04-01 2016 4 1 Breeding 2016Breeding 14
## 2 2016-04-02 2016 4 2 Breeding 2016Breeding 14
## 3 2016-04-04 2016 4 4 Breeding 2016Breeding 14
## 4 2016-04-18 2016 4 18 Breeding 2016Breeding 16
## 5 2016-05-14 2016 5 14 Breeding 2016Breeding 20
## 6 2016-05-14 2016 5 14 Breeding 2016Breeding 20Although the abundance estimates shown above cover an entire season, migration timing can be approximated more precisely. For example, here we calculate and then graph the average weekly latitude of birds from the orginal eBird data. The vertical dotted lines shwon below highlight weeks 7 through 13 which coincide with the spring season. As indicated by the graph, NORS average latitude varies within the spring season; it steadily increases. The variabilty shown by this data can be used to estimate the most probable latitude for any week of the year.
Avg.lat = as.data.frame(
NORS %>%
group_by(Week) %>%
summarise(Average = mean(Lat)))
ggplot(Avg.lat,
aes(Week, Average)) +
geom_line(col = "black",
lwd=1) +
geom_vline(xintercept = c(7,13),
linetype = "dotted",
col = "red",
size = 0.5) +
xlab("Week of Year") +
ylab(expression("Average Latitude")) +
scale_x_continuous(breaks = seq(0, 50, 10),
labels = paste(seq(0, 50, 10)),
limits = c(0,53)) +
theme_classic() +
theme(axis.text.x = element_text(face="bold", size=12),
axis.text.y = element_text(face="bold", size=14),
strip.text = element_text(face="bold", size = 20),
axis.title.y = element_text(size = 20),
axis.title.x = element_text(size = 20, vjust=-2))
Example (Week 7) Before running the full analysis, let’s look at Week 7 as an example of how the process would work.
First, we create a copy of the spring abindnace raster and make a dense point grid version of the data. Essentially, we identify the lat and long of each cell in the abundance raster.
NORS.grd = rasterToPoints(NORS.spring, sp = TRUE)
NORS.grd@data = NORS.grd@data %>%
mutate(Long = NORS.grd@coords[,1],
Lat = NORS.grd@coords[,2])Next, we use the eBird data to find the probability of each possible latitude in the point grid above (representing the whole U.S.). We find that probability through comparison with the average and standard deviation (SD) of eBird observations made in a given week (i.e., the latitudes of observation locations).
mySD = 0.25 #Precision threshold. The number of SD's to approximate. Higher values will produce broader coverage, lower values more narrow.
NORS.spr.df = NORS %>% filter(Week >= 7 & Week <= 13) #Extract the spring season from the full data set
Week.lvls = levels(factor(NORS.spr.df$Week)) #List of weeks to model
weekly.obs = NORS.spr.df %>% filter(Week == Week.lvls[1]) #take the first week in the set, Week 7
avg.rng = mean(weekly.obs$Lat) #Find the average Latitude for Week 7
sd.rng = sd(weekly.obs$Lat) #find the SD of the Week 7 latitudes
Top.lim = avg.rng+(sd.rng*mySD) #Add the precision threhsold entered above from the average (the new "max" latitude)
Bot.lim = avg.rng-(sd.rng*mySD) #As above, but now subtracting from the average to get the "minimum" latitude
#Lower probability, Finding probabilty that actual is the given latitude or LOWER
NORS.grd$ProbL = pnorm(NORS.grd$Lat, #NORS.grd includes all latitudes in the U.S.; pnorm assumes a normal distribution.
mean = Bot.lim, #Minimum lat calculated from eBird data with precision
sd = mySD*sd.rng, #Updated sd
lower.tail = TRUE) #Finding probabilty that actual is the given latitude or lower
NORS.grd$ProbL[NORS.grd$ProbL > 0.5] = 0 #Set less than 50% to zero, because it's calculated seperately below based on our precision.
#As above, but now find probability that the actual is equal or HIGHER than the mean
NORS.grd$ProbH = pnorm(NORS.grd$Lat,
mean = Top.lim,
sd = mySD*sd.rng,
lower.tail = FALSE) #Don't want the lower tail
NORS.grd$ProbH[NORS.grd$ProbH > 0.5] = 0
#Finally, values are mapped to the point grid (combining lower and higher probabilties to one map)
NORS.grd$ProbC = Scale(ifelse(NORS.grd$Lat <= Top.lim & NORS.grd$Lat >= Bot.lim, 0.5,
ifelse(NORS.grd$Lat < Bot.lim, NORS.grd$ProbL,
ifelse(NORS.grd$Lat > Top.lim, NORS.grd$ProbH, 0))))
#Create a raster of the result
tmp.r = rasterize(NORS.grd,
NORS.spring,
"ProbC",
background = NA)View the example result. Recall, this map shows the probability that locations (latitudes) in the U.S. are equal to the average Latitude of NORS during Week 7 of the year within 0.5 of a SD (2 x the mySD threshold specified above).
plot(tmp.r)
Now, the probabilty of a bird being at a latitude is multiplied times the orginally estimated abundance.
tmp2.r = tmp.r*NORS.spring #Abundance at a location times the probabilty of being at that latitudeCompare the Week 7 estimate to the whole spring season
Exampl.stk =stack(NORS.spring, tmp2.r)
names(Exampl.stk) = c("Whole.Season", "Week.7")
breaks = seq(0, 101, 0.1)
mCols = brewer.pal(9, "YlOrRd")
cr = colorRampPalette(mCols)(length(breaks)-1)
cr[1] = "lightgray"
cr = colorRampPalette(cr,
bias = 2.5, space = "rgb")
M0 = levelplot(Exampl.stk,
margin = FALSE,
layout=c(1, 2),
xlab = "NORS ABundance",
ylab = NULL,
maxpixels = 1e5,
col.regions = cr, at = breaks,
colorkey = list(labels=list(at=c(0, 20, 40, 60, 80, 100),
labels=c("0", "20", "40", "60", "80", "100 +"), cex=1.5),
labels=list(cex=12),
space = "right"),
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.5)) +
latticeExtra::layer(sp.polygons(States, col = "darkgray", lwd = 0.5))
M0
3 Model Full Year
In two week intervals.
Load Poultry Data Commercial Poultry Proximity (USDA Hybrid Model). Just using broilers in this example.
Hybrid.df = read.csv("C:/Users/humph173/Documents/LabWork/Patuxent/HybridModel/broilers_hybrid_FINAL.csv",
stringsAsFactors=FALSE) %>%
select(POINT_X, POINT_Y, population, commodityType, category)
Hybrid.pnt = SpatialPointsDataFrame(Hybrid.df[,c("POINT_X", "POINT_Y")], Hybrid.df)
proj4string(Hybrid.pnt) = proj4string(States)Create raster Version
Blank.r = NORS.spring
Blank.r[Blank.r >= 0] = 0
Broilers.r = rasterize(Hybrid.pnt,
Blank.r,
"population",
background = 0)
Broilers.r = Broilers.r + Blank.r #Outside US to NA
plot(Broilers.r) 
Loop through Full Year In two week intervals.
#Load additional seasons
NORS.brd = raster("E:/Dabbler_Pre/Dabbler_project/Raster_outputs/abun_nors_Breeding.tif")
NORS.fal = raster("E:/Dabbler_Pre/Dabbler_project/Raster_outputs/abun_nors_Fall.tif")
NORS.wint = raster("E:/Dabbler_Pre/Dabbler_project/Raster_outputs/abun_nors_Winter.tif")
NORS.spr.df = NORS
Week.lvls = seq(1, 52, 2)
for(i in 1:length(Week.lvls)){
weekly.obs = NORS.spr.df %>% filter(Week == Week.lvls[i])
avg.rng = mean(weekly.obs$Lat)
sd.rng = sd(weekly.obs$Lat)
Top.lim = avg.rng+(sd.rng*mySD)
Bot.lim = avg.rng-(sd.rng*mySD)
NORS.grd$ProbL = pnorm(NORS.grd$Lat,
mean = Bot.lim,
sd = mySD*sd.rng,
lower.tail = TRUE)
NORS.grd$ProbL[NORS.grd$ProbL > 0.5] = 0
NORS.grd$ProbH = pnorm(NORS.grd$Lat,
mean = Top.lim,
sd = mySD*sd.rng,
lower.tail = FALSE)
NORS.grd$ProbH[NORS.grd$ProbH > 0.5] = 0
NORS.grd$ProbC = Scale(ifelse(NORS.grd$Lat <= Top.lim & NORS.grd$Lat >= Bot.lim, 0.5,
ifelse(NORS.grd$Lat < Bot.lim, NORS.grd$ProbL,
ifelse(NORS.grd$Lat > Top.lim, NORS.grd$ProbH, 0))))
tmp.r = rasterize(NORS.grd,
NORS.spring,
"ProbC",
background = NA)
if(i > 48 | i <= 6){tmp2.r = tmp.r*NORS.wint*0.09*Broilers.r} #Nonsense risk formula
if(i > 6 & i <= 13){tmp2.r = tmp.r*NORS.spring*0.1*Broilers.r}
if(i > 13 & i <= 35){tmp2.r = tmp.r*NORS.brd*0.15*Broilers.r}
if(i > 35 & i <= 48){tmp2.r = tmp.r*NORS.fal*0.2*Broilers.r}
values(tmp2.r) = Scale(values(tmp2.r)) #scaling
if(i == 1){Full.Year = tmp2.r
} else{Full.Year = stack(Full.Year, tmp2.r)}
}
names(Full.Year) = paste("Week ", Week.lvls, sep="")4 Relative Risk Maps
Showing every two weeks.
breaks = seq(0, 1, 0.01)
mCols = brewer.pal(9, "YlOrRd")
cr = colorRampPalette(mCols)(length(breaks)-1)
cr[1] = "lightgray"
cr = colorRampPalette(cr,
bias = 2.5, space = "rgb")
Full.Year[Full.Year == 0] = NA
for(i in 1:dim(Full.Year)[3]){
myPlt = levelplot(Full.Year[[i]],
#layout=c(1, dim(Full.Year)[3]),
margin = TRUE,
xlab = paste(names(Full.Year))[i],
ylab = NULL,
#names.attr= paste(names(Full.Year)),
maxpixels = 1e5,
col.regions = cr, at = breaks,
colorkey = list(labels=list(at=seq(0, 1, 0.1),
labels=paste(seq(0, 1, 0.1)), cex=1.5),
labels=list(cex=12),
space = "right"),
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.5)) +
latticeExtra::layer(sp.polygons(States, col = "darkgray", lwd = 0.5))
plot(myPlt)
}
