Spatial analyses
Gunnar R. Kramer
2023-01-13
Load libraries:
require(rgdal)
require(rgeos)
require(sp)
require(ggplot2)
require(raster)
library(maps)
library(maptools)
require(devtools)
library(spatialEco)
library(OneR)
library(pals)
library(scico)
library(viridis)Load functions
img <- function(obj, nam) {
image(1:length(obj), 1, as.matrix(1:length(obj)), col=obj,
main = nam, ylab = "", xaxt = "n", yaxt = "n", bty = "n")
}
color.bar <- function(lut, min, max=-min, nticks=11, ticks=seq(min, max, len=nticks), title='') {
scale = (length(lut)-1)/(max-min)
dev.new(width=1.75, height=5)
plot(c(0,10), c(min,max), type='n', bty='n', xaxt='n', xlab='', yaxt='n', ylab='', main=title)
axis(2, ticks, las=1)
for (i in 1:(length(lut)-1)) {
y = (i-1)/scale + min
rect(0,y,10,y+1/scale, col=lut[i], border=NA)
}
}Load relevant files
#Load risk-factor rasters (standardized extent, resolution, etc.)
setwd("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/")
#set color ramp for raster plotting
cramp1<-colorRampPalette(c("white","snow2","lightgray","snow3","gold","goldenrod1","hotpink","darkorchid1","darkorchid4"))
#view
color.bar(colorRampPalette(c("white","snow2","lightgray","snow3","gold","goldenrod1","hotpink","darkorchid1","darkorchid4"))(50), min=0, max=1)
#load deployment site shapefile
sites<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/DeploymentSites_4R.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/DeploymentSites_4R.shp", layer: "DeploymentSites_4R"
## with 28 features
## It has 11 fields
## Integer64 fields read as strings: n ngwFall nbwFall nhy ngwSpr nbwSpr#housekeeping of site layer
sites<-spTransform(sites, CRS("+proj=longlat +datum=WGS84"))
breaks<-as.numeric(c(1, 2, 5, 17))
tags<-as.character(c(1, 2, 3))
sites@data$sitesizegwFall<-as.numeric(as.character(sites@data$ngwFall))
sites@data$sitesizebwFall<-as.numeric(as.character(sites@data$nbwFall))
sites@data$sitesizebwBINFall<-cut(sites@data$sitesizebwFall, breaks=breaks, label=tags, include.lowest=TRUE)
sites@data$sitesizegwBINFall<-cut(sites@data$sitesizegwFall, breaks=breaks, labels=tags, include.lowest=TRUE)
sites@data$sitesizegwSpr<-as.numeric(as.character(sites@data$ngwSpr))
sites@data$sitesizebwSpr<-as.numeric(as.character(sites@data$nbwSpr))
sites@data$sitesizebwBINSpr<-cut(sites@data$sitesizebwSpr, breaks=breaks, label=tags, include.lowest=TRUE)
sites@data$sitesizegwBINSpr<-cut(sites@data$sitesizegwSpr, breaks=breaks, labels=tags, include.lowest=TRUE)
trend.breaks<-as.numeric(c(-10, -5, -1.5, 1.5, 5, 100))
trend.tags<-as.character(c(-5, -1.5, 1.5, 5, 10))
sites@data$sitetrendgw<-as.numeric(as.character(sites@data$gwtr))
sites@data$sitetrendbw<-as.numeric(as.character(sites@data$bwtr))
sites@data$gwTrendBin<-cut(sites@data$sitetrendgw, breaks=trend.breaks, label=trend.tags, include.lowest=TRUE)
sites@data$bwTrendBin<-cut(sites@data$sitetrendbw, breaks=trend.breaks, label=trend.tags, include.lowest=TRUE)
#set colors based on population trends
sites@data$sitetrendgwCOL<-ifelse(sites@data$sitetrendgw< -5, "#7E1900", ifelse(sites@data$sitetrendgw >= -5 & sites@data$sitetrendgw < -1.5, "#B48A2C", ifelse(sites@data$sitetrendgw >= -1.5 & sites@data$sitetrendgw < 1.5, "#E5E598", ifelse(sites@data$sitetrendgw >= 1.5 & sites@data$sitetrendgw < 5, "#8CDED9", ifelse(sites@data$sitetrendgw >5, "#3F85BB", 0)))))
sites@data$sitetrendbwCOL<-ifelse(sites@data$sitetrendbw< -5, "#7E1900", ifelse(sites@data$sitetrendbw >= -5 & sites@data$sitetrendbw < -1.5, "#B48A2C", ifelse(sites@data$sitetrendbw >= -1.5 & sites@data$sitetrendbw < 1.5, "#E5E598", ifelse(sites@data$sitetrendbw >= 1.5 & sites@data$sitetrendbw < 5, "#8CDED9", ifelse(sites@data$sitetrendbw >5, "#3F85BB", 0)))))
bwtrendcolList<-sites@data$sitetrendbwCOL[!is.na(sites@data$sitetrendbwCOL)]
gwtrendcolList<-sites@data$sitetrendgwCOL[!is.na(sites@data$sitetrendgwCOL)]
#assign sites to BCRs, assign colors based on BCR
sites@data$BCRbw<-c("BHT", "CH", "CH", "AM", "AM", "PHT", "PHT", "BHT", "AM", "AM", "AM", "BHT", "AM", "BHT", "BHT", "BHT", "BHT", "AM", "CH", "AM", "AM", "PHT", "AM", "OTHER","OTHER","OTHER","OTHER","OTHER")
sites@data$BCRgw<-c("BHT", "CH", "CH", "AM", "AM", "BHT", "BHT", "BHT", "AM", "AM", "BHT", "BHT", "AM", "BHT", "BHT", "BHT","BHT", "AM", "CH", "AM", "AM", "BHT", "AM", "OTHER","OTHER","OTHER","OTHER","OTHER")
sites@data$BCRcolBW<-ifelse(sites@data$BCRbw=="CH", rgb(255/255, 9/255, 231/255, 1), ifelse(sites@data$BCRbw=="PHT", rgb(0/255, 146/255, 146/255, 1), ifelse(sites@data$BCRbw=="AM", rgb(255/255, 157/255, 9/255), "NA")))
sites@data$BCRcolGW<-ifelse(sites@data$BCRgw=="BHT", rgb(146/255,0/255,0/255, 1), ifelse(sites@data$BCRgw=="AM", rgb(255/255, 157/255, 9/255), "NA"))
#load in wintering grounds point shapefile data
bennetpts<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BennetBirdCoords4R.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BennetBirdCoords4R.shp", layer: "BennetBirdCoords4R"
## with 20 features
## It has 3 fieldsbennetpts<-spTransform(bennetpts, CRS("+proj=longlat +datum=WGS84"))
bwcexlist<-subset(sites, as.numeric(as.character(nbwFall))>0)
gwcexlist<-subset(sites, as.numeric(as.character(ngwFall))>0)
summary(gwcexlist)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 -100.58100 -75.180
## coords.x2 10.21512 51.535
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Number of points: 19
## Data attributes:
## lat lon sp n
## Min. :10.22 Min. :-100.58 Length:19 Length:19
## 1st Qu.:25.83 1st Qu.: -93.55 Class :character Class :character
## Median :43.95 Median : -86.03 Mode :character Mode :character
## Mean :36.54 Mean : -87.33
## 3rd Qu.:46.02 3rd Qu.: -81.95
## Max. :51.53 Max. : -75.18
##
## ngwFall nbwFall nhy ngwSpr
## Length:19 Length:19 Length:19 Length:19
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## nbwSpr gwtr bwtr sitesizegwFall
## Length:19 Min. :-8.710 Min. :-1.320 Min. : 1.000
## Class :character 1st Qu.:-7.620 1st Qu.:-1.050 1st Qu.: 1.000
## Mode :character Median : 0.370 Median : 1.055 Median : 2.000
## Mean : 1.873 Mean : 1.817 Mean : 3.316
## 3rd Qu.: 1.740 3rd Qu.: 3.922 3rd Qu.: 3.500
## Max. :20.100 Max. : 6.480 Max. :17.000
## NA's :5 NA's :15
## sitesizebwFall sitesizebwBINFall sitesizegwBINFall sitesizegwSpr
## Min. :0.0000 1 : 2 1:13 Min. : 0.000
## 1st Qu.:0.0000 2 : 2 2: 3 1st Qu.: 1.000
## Median :0.0000 3 : 0 3: 3 Median : 2.000
## Mean :0.5263 NA's:15 Mean : 3.158
## 3rd Qu.:0.0000 3rd Qu.: 3.500
## Max. :4.0000 Max. :13.000
##
## sitesizebwSpr sitesizebwBINSpr sitesizegwBINSpr sitetrendgw
## Min. :0.0000 1 : 2 1 :12 Min. :-8.710
## 1st Qu.:0.0000 2 : 2 2 : 3 1st Qu.:-7.620
## Median :0.0000 3 : 0 3 : 3 Median : 0.370
## Mean :0.5263 NA's:15 NA's: 1 Mean : 1.873
## 3rd Qu.:0.0000 3rd Qu.: 1.740
## Max. :4.0000 Max. :20.100
## NA's :5
## sitetrendbw gwTrendBin bwTrendBin sitetrendgwCOL sitetrendbwCOL
## Min. :-1.320 -5 :5 -5 : 0 Length:19 Length:19
## 1st Qu.:-1.050 -1.5:0 -1.5: 0 Class :character Class :character
## Median : 1.055 1.5 :3 1.5 : 2 Mode :character Mode :character
## Mean : 1.817 5 :3 5 : 1
## 3rd Qu.: 3.922 10 :3 10 : 1
## Max. : 6.480 NA's:5 NA's:15
## NA's :15
## BCRbw BCRgw BCRcolBW BCRcolGW
## Length:19 Length:19 Length:19 Length:19
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## #shapefile layers are not critical for anlaysis. they are used for visualization purposes. geographic layers can be downloaded from many sources, bird range maps can be downloaded from BirdConservation International
setwd("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/")
canada<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Canada.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/Canada.shp", layer: "Canada"
## with 13 features
## It has 2 fieldscanada<-spTransform(canada, CRS("+proj=longlat +datum=WGS84"))
#plot(canada,add=T)
bwwaBR<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BWWA_range.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BWWA_range.shp", layer: "BWWA_range"
## with 2 features
## It has 17 fields
## Integer64 fields read as strings: OBJECTID_1 SISID PRESENCE ORIGIN SEASONALbwwaBR<-spTransform(bwwaBR, CRS("+proj=longlat +datum=WGS84"))
#plot(bwwaBR,add=T)
gwwaBR<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters//gwwa_new_range_28JUL11_final.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/gwwa_new_range_28JUL11_final.shp", layer: "gwwa_new_range_28JUL11_final"
## with 4 features
## It has 1 fieldsgwwaBR<-spTransform(gwwaBR, CRS("+proj=longlat +datum=WGS84"))
#plot(gwwaBR,add=T)
bwwaNB<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BWNB_Dissolved.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BWNB_Dissolved.shp", layer: "BWNB_Dissolved"
## with 1 features
## It has 14 fields
## Integer64 fields read as strings: FID_bw_NB_ FID_Whemis POP2005bwwaNB<-spTransform(bwwaNB, CRS("+proj=longlat +datum=WGS84"))
#plot(bwwaNB,add=T)
gwwaNB<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/NBchrysRange100km_proj.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/NBchrysRange100km_proj.shp", layer: "NBchrysRange100km_proj"
## with 25 features
## It has 12 fieldsgwwaNB<-spTransform(gwwaNB, CRS("+proj=longlat +datum=WGS84"))
#plot(gwwaNB)
AM_BCR<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/AppMtn_BCR.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/AppMtn_BCR.shp", layer: "AppMtn_BCR"
## with 1 features
## It has 5 fieldsAM_BCR<-spTransform(AM_BCR, CRS("+proj=longlat +datum=WGS84"))
#plot(AM_BCR)
BHT_BCR<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BHTrans_BCR.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/BHTrans_BCR.shp", layer: "BHTrans_BCR"
## with 1 features
## It has 5 fieldsBHT_BCR<-spTransform(BHT_BCR, CRS("+proj=longlat +datum=WGS84"))
#plot(BHT_BCR)
CH_BCR<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/CH_BCR.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/CH_BCR.shp", layer: "CH_BCR"
## with 1 features
## It has 5 fieldsCH_BCR<-spTransform(CH_BCR, CRS("+proj=longlat +datum=WGS84"))
#plot(CH_BCR)
PH_BCR<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/PH_BCR.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/PH_BCR.shp", layer: "PH_BCR"
## with 1 features
## It has 5 fieldsPH_BCR<-spTransform(PH_BCR, CRS("+proj=longlat +datum=WGS84"))
#plot(PH_BCR)
winterpts<-readOGR("/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/winterPopCoords.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/Users/gunnarkramer/Desktop/FINAL DATA FOR DRUM LAND ECOL/Rasters/winterPopCoords.shp", layer: "winterPopCoords"
## with 90 features
## It has 7 fields
## Integer64 fields read as strings: spr fallwinterpts<-spTransform(winterpts, CRS("+proj=longlat +datum=WGS84"))
#plot(winterpts, add=TRUE, pch=5, col="black")
winterpts@data$color<-ifelse(winterpts$pop=="PHT", rgb(0/255, 146/255, 146/255, 1), ifelse(winterpts$pop=="CH", rgb(255/255, 9/255, 231/255, 1), ifelse(winterpts$pop=="bwAM"|winterpts$pop=="AM", rgb(255/255, 157/255, 9/255, 1), rgb(146/255,0/255,0/255, 1))))
winterpts@data$pch1<-ifelse(winterpts$seas=="wint", 4, 3)
rlist=list.files(getwd(), pattern="tif$", full.names=FALSE)
for(i in rlist) { assign(unlist(strsplit(i,"[.]"))[1], raster(i)) }
data(wrld_simpl)
#####
#colorblind friendly pallete
black<-rgb(0/255,0/255,0/255, 1)
darkgreen<-rgb(0/255,0/73,0/73, 1)
teal<-rgb(0/255,146/255,146/255, 1)
pink<-rgb(255/255,109/255,182/255, 1)
lightpink<-rgb(255/255,182/255,219/255, 1)
purple<-rgb(73/255,0/255,146/255, 1)
darkblue<-rgb(0/255,109/255,219/255, 1)
lightpurple<-rgb(182/255,109/255,255/255, 1)
lightblue<-rgb(109/255,182/255,255/255, 1)
paleblue<-rgb(182/255,219/255,255/255, 1)
maroon<-rgb(146/255,0/255,0/255, 1)
brown<-rgb(146/255,73/255,0/255, 1)
orange<-rgb(219/255,109/255,0/255, 1)
green<-rgb(36/255,255/255,36/255, 1)
yellow<-rgb(255/255,255/255,109/255, 1)
####Fig 3 - Range maps with population trends and study sites
#inset
#png("BWWA_poptrendmap.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
par(mar=c(5,4,4,5.5))
plot(1,1, pty="s", xlim=c(-116,-54), ylim=c(0,60), main="BWWA Sites", ylab="", xlab="", yaxs="i")
plot(bwwaBR, add=TRUE, border="gray65", col="gray65", lwd=.75, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", col="gray45", lwd=1, lty=1)
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(CH_BCR, add=TRUE, border=rgb(255/255, 9/255, 231/255, 1), col=rgb(255/255, 9/255, 231/255, .3), lwd=.75, lty=1)
plot(AM_BCR, add=TRUE, border=rgb(255/255, 157/255, 9/255, 1), col=rgb(255/255, 157/255, 9/255, .3), lwd=.75, lty=1)
plot(PH_BCR, add=TRUE, border=rgb(0/255, 146/255, 146/255, 1), col=rgb(0/255, 146/255, 146/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(nbwFall))>0), pch=21, bg=c(bwtrendcolList), lwd=.5, col="gray20", cex=as.numeric(as.character(bwcexlist@data$sitesizebwBINFall)))
#points(subset(winterpts, as.character(sp)=="BW"), pch=4, col=subset(winterpts$color, winterpts$sp=="BW"), lwd=1, cex=.5)
box()
#dev.off()
#png("GWWA_poptrendmapBCRgray_reduced_noBCR.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
par(mar=c(5,4,4,5.5))
plot(1,1, pty="s", xlim=c(-116,-54), ylim=c(0,60), main="GWWA Sites", ylab="", xlab="", yaxs="i")
plot(gwwaBR, add=TRUE, border="gray65", col="gray65", lwd=.75, lty=1)
plot(gwwaNB, add=TRUE, border="gray45", col="gray45", lwd=1, lty=1)
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(AM_BCR, add=TRUE, border=rgb(255/255, 157/255, 9/255, 1), col=rgb(255/255, 157/255, 9/255, .3), lwd=.75, lty=1)
plot(BHT_BCR, add=TRUE, border=rgb(146/255,0/255,0/255, 1), col=rgb(146/255,0/255,0/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(ngwFall))>0), pch=21, bg=c(gwtrendcolList), lwd=.5, col="gray20", cex=as.numeric(as.character(gwcexlist@data$sitesizegwBINFall)))
sites$sitetrendgwCOL## [1] "#3F85BB" NA NA NA NA "#7E1900" NA
## [8] "#8CDED9" NA NA "#E5E598" "#E5E598" "#7E1900" "#3F85BB"
## [15] "#8CDED9" "#3F85BB" "#8CDED9" NA NA "#7E1900" "#7E1900"
## [22] "#E5E598" "#7E1900" NA NA NA NA NA#points(subset(winterpts, as.character(sp)=="GW"), pch=subset(winterpts$pch1, winterpts$sp=="GW"), col=subset(winterpts$color, winterpts$sp=="GW"), lwd=1, cex=.5)
#dev.off()Plot population-specific space use maps (based on Bird Conservation Regions)
Blue-winged warblers, PHT, Fall
#####
#BWfallPHT
######
BWWA.Fall.PHT<-(WIB05_fall_x2_clip + MIB02_fall_x2_clip+ MIB05_fall_x2_clip+ MIB12_fall_x2_clip+ MIB14_fall_x2_clip+ WIB11_2016_fall_x2_clip+ WIB11_2017_fall_x2_clip + WIB19_fall_x2_clip)
BWWA.Fall.PHT.agg<-aggregate(BWWA.Fall.PHT, fact=1.5)
BWWA.Fall.PHT.agg<-BWWA.Fall.PHT.agg/cellStats(BWWA.Fall.PHT.agg, sum)
######
#png("BW_PHT_Fall_BCRtrend.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(BWWA.Fall.PHT.agg, col=cramp1(50), main="BWWA PHT Fall")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(bwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(PH_BCR, add=TRUE, border="black", lwd=.75, lty=1)
#plot(PH_BCR, add=TRUE, border=rgb(0/255, 146/255, 146/255, 1), col=rgb(0/255, 146/255, 146/255, .3), lwd=.75, lty=1)
#points(subset(sites, as.numeric(as.character(nbw))>0 & lat>40 & lon< -80), pch=21, bg=gray(.3, .2), col="black", cex=.75)
#points(subset(sites, as.numeric(as.character(nbwFall))>0 & lat>40 & lon< -80), pch=21, bg=c(bwtrendcolList[c(5,6,13)]), col="black", cex=as.numeric(as.character(bwcexlist@data[c(5,6,13),14])))
points(subset(sites, as.numeric(as.character(nbwFall))>0 & lat>40 & lon< -80), pch=21, bg=rgb(0/255, 146/255, 146/255, 1))
#points(subset(winterpts, as.character(pop)=="PHT" & fall==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="PHT"), lwd=1, cex=1)
#dev.off()Blue-winged warbler, CH, fall
#BWfallCH
#######
BWWA.Fall.CH<-(ILB04_fall_x2_clip + ILB06_fall_x2_clip+ ILB08_fall_x2_clip+ ILB18_fall_x2_clip + TNB09_fall_x2_clip)
BWWA.Fall.CH.agg<-aggregate(BWWA.Fall.CH, fact=1.5)
BWWA.Fall.CH.agg<-BWWA.Fall.CH.agg/cellStats(BWWA.Fall.CH.agg, sum)
#######
#png("BW_CH_Fall_BCRtrend.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(BWWA.Fall.CH.agg, col=cramp1(50), main="BWWA CH Fall")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(bwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(CH_BCR, add=TRUE, border=rgb(255/255, 9/255, 231/255, 1), col=rgb(255/255, 9/255, 231/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(nbwFall))>0 & lat<40 & lon< -85), pch=21, bg=rgb(255/255, 9/255, 231/255, 1), col="black")
#points(subset(winterpts, as.character(pop)=="CH" & fall==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="CH"), lwd=1, cex=1)
#points(subset(sites, as.numeric(as.character(nbw))>0 & lat<40 & lon< -85), pch=21, bg=gray(.3, .2), col="black", cex=.75)
#dev.off()Blue-winged warblers, AM, fall
#BWfallAM
######
BWWA.Fall.AM<-(TNB03_fall_x2_clip + KYB14_fall_x2_clip+ MAB05_fall_x2_clip+ MAB06_fall_x2_clip+ MAB09_fall_x2_clip+ MAB12_fall_x2_clip+ NCB01_spring_x2_clip+ ONB02_fall_x2_clip+ ONB05_fall_x2_clip+ PAB01_fall_x2_clip+ PAB03_fall_x2_clip+ PAB05_fall_x2_clip+ PAB07_fall_x2_clip)
BWWA.Fall.AM.agg<-aggregate(BWWA.Fall.AM, fact=1.5)
BWWA.Fall.AM.agg<-BWWA.Fall.AM.agg/cellStats(BWWA.Fall.AM.agg, sum)
#######
#png("BW_AM_Fall_BCRtrend.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(BWWA.Fall.AM.agg, col=cramp1(50), main="BWWA AM Fall")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(bwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(AM_BCR, add=TRUE, border=rgb(255/255, 157/255, 9/255, 1), col=rgb(255/255, 157/255, 9/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(nbwFall))>0 & lat<40 & lon> -85), pch=21, bg=rgb(255/255, 157/255, 9/255, 1), col="black")
points(subset(sites, as.numeric(as.character(nbwFall))>0 & lon> -80), pch=21, bg=rgb(255/255, 157/255, 9/255, 1), col="black")
#points(subset(winterpts, as.character(pop)=="bwAM" & fall==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="bwAM"), lwd=1, cex=1)
#points(subset(sites, as.numeric(as.character(nbw))>0 & lat<40 & lon< -85), pch=21, bg=gray(.3, .2), col="black", cex=.75)
#dev.off()Blue-winged warblers, PHT, spring
#BW_SPRING_PHT
######
BWWA.Spring.PHT<-(WIB05_spring_x2_clip + MIB02_spring_x2_clip+ MIB12_spring_x2_clip+ MIB14_spring_x2_clip + WIB11_2016_spring_x2_clip+ WIB11_2017_spring_x2_clip + WIB19_spring_x2_clip)
BWWA.Spring.PHT.agg<-aggregate(BWWA.Spring.PHT, fact=1.5)
BWWA.Spring.PHT.agg<-BWWA.Spring.PHT.agg/cellStats(BWWA.Spring.PHT.agg, sum)
######
#png("BW_PHT_Spring_BCRtrend.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(BWWA.Spring.PHT.agg, col=cramp1(50), main="BWWA PHT Spring")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(bwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(PH_BCR, add=TRUE, border=rgb(0/255, 146/255, 146/255, 1), col=rgb(0/255, 146/255, 146/255, .3), lwd=.75, lty=1)
#points(subset(sites, as.numeric(as.character(nbw))>0 & lat>40 & lon< -80), pch=21, bg=gray(.3, .2), col="black", cex=.75)
points(subset(sites, as.numeric(as.character(nbwSpr))>0 & lat>40 & lon< -80), pch=21, bg=rgb(0/255, 146/255, 146/255, 1), col="black")
#points(subset(winterpts, as.character(pop)=="PHT" & spr==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="PHT"), lwd=1, cex=1)
#dev.off()Blue-winged warblers, CH, spring
#BW CH SPRING
######
BWWA.Spring.CH<-(ILB04_spring_x2_clip + ILB06_spring_x2_clip+ ILB08_spring_x2_clip+ ILB18_spring_x2_clip+ TNB09_spring_x2_clip)
BWWA.Spring.CH.agg<-aggregate(BWWA.Spring.CH, fact=1.5)
BWWA.Spring.CH.agg<-BWWA.Spring.CH.agg/cellStats(BWWA.Spring.CH.agg, sum)
#####
#png("BW_CH_Spring_BCRtrend.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(BWWA.Spring.CH.agg, col=cramp1(50), main="BWWA CH Spring")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(bwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
#plot(CH_BCR, add=TRUE, border=rgb(255/255, 9/255, 231/255, 1), col=rgb(255/255, 9/255, 231/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(nbwSpr))>0 & lat<40 & lon< -85), pch=21, bg=rgb(255/255, 9/255, 231/255, 1), col="black")
#points(subset(winterpts, as.character(pop)=="CH" & spr==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="CH"), lwd=1, cex=1)
#points(subset(sites, as.numeric(as.character(nbw))>0 & lat<40 & lon< -85), pch=21, bg=gray(.3, .2), col="black", cex=.75)
#dev.off()Blue-winged warblers, AM, spring
#bw spring am
######
BWWA.Spring.AM<-(MAB05_spring_x2_clip+ MAB06_spring_x2_clip+ MAB09_spring_x2_clip+ MAB12_spring_x2_clip + NCB01_spring_x2_clip+ ONB02_spring_x2_clip+ ONB05_spring_x2_clip+ PAB01_spring_x2_clip+ PAB03_spring_x2_clip+ PAB05_spring_x2_clip+ PAB07_spring_x2_clip)
BWWA.Spring.AM.agg<-aggregate(BWWA.Spring.AM, fact=1.5)
BWWA.Spring.AM.agg<-BWWA.Spring.AM.agg/cellStats(BWWA.Spring.AM.agg, sum)
######
#png("BW_AM_Spring.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(BWWA.Spring.AM.agg, col=cramp1(50), main="BWWA AM Spring")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(bwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(AM_BCR, add=TRUE, border=rgb(255/255, 157/255, 9/255, 1), col=rgb(255/255, 157/255, 9/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(nbwSpr))>0 & lat<40 & lon>-82), pch=21, bg=rgb(255/255, 157/255, 9/255, 1), col="black")
points(subset(sites, as.numeric(as.character(nbwSpr))>0 & lon> -80), pch=21, bg=rgb(255/255, 157/255, 9/255, 1), col="black")
#points(subset(winterpts, as.character(pop)=="bwAM" & spr==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="bwAM"), lwd=1, cex=1)
#dev.off()Golden-winged warblers, BHT, fall
#gw gl fall
######
GWWA.Fall.GL<-(DMG03_fall_x2_clip + MIG01_fall_x2_clip + RLG12_fall_x2_clip + RLG19_fall_x2_clip + RLG18_fall_x2_clip + ONG03_fall_x2_clip + ONG05_fall_x2_clip + ONG10_fall_x2_clip + RLG16_fall_x2_clip + RLG23_fall_x2_clip + RLGF15_fall_x2_clip + RMG02_fall_x2_clip + SHG07_fall_x2_clip + SLG12_fall_x2_clip + SLG18_fall_x2_clip + TAG07_fall_x2_clip + TAG09_fall_x2_clip + TAG14_fall_x2_clip + WIG01_fall_x2_clip + WIG02_fall_x2_clip + MN05_fall_x2_clip + MN06_fall_x2_clip + MN11_fall_x2_clip + MN12_fall_x2_clip + MN14_fall_x2_clip + MN15_fall_x2_clip + MN16_fall_x2_clip + MN20_fall_x2_clip + MN25_fall_x2_clip + MN29_fall_x2_clip + MN36_fall_x2_clip + B1_fall_x2_clip + B2_fall_x2_clip + B3_fall_x2_clip + B4_fall_x2_clip + B5_fall_x2_clip + B6_fall_x2_clip + B7_fall_x2_clip + B8_fall_x2_clip + B9_fall_x2_clip + B10_fall_x2_clip + B11_fall_x2_clip + B12_fall_x2_clip + B13_fall_x2_clip + B14_fall_x2_clip + B15_fall_x2_clip + B16_fall_x2_clip + B17_fall_x2_clip + B18_fall_x2_clip + B20_fall_x2_clip)
GWWA.Fall.GL.agg<-aggregate(GWWA.Fall.GL, fact=1.5)
GWWA.Fall.GL.agg<-GWWA.Fall.GL.agg/cellStats(GWWA.Fall.GL.agg, sum)
######
#png("GW_GL_Fall.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(GWWA.Fall.GL.agg, col=cramp1(50), main="GWWA GL Fall")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(bwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(bwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(BHT_BCR, add=TRUE, border=rgb(146/255,0/255,0/255, 1), col=rgb(146/255,0/255,0/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(ngwFall))>0 & lat>43), pch=21, bg=rgb(146/255,0/255,0/255, 1), col="black")
points(subset(sites, as.numeric(as.character(ngwFall))>0 & lat<20), pch=21, bg=rgb(146/255,0/255,0/255, 1))
#points(subset(sites, as.numeric(as.character(ngwFall))>0 & lon> -80), pch=21, bg=c(gwtrendcolList[c(4,8,9,10)]), col="black", cex=as.numeric(as.character(gwcexlist@data[c(4,8,9,10),15])))
#points(subset(winterpts, as.character(pop)=="BHT" & fall==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="BHT"), lwd=1, cex=1)
#dev.off()Golden-winged warbler, AM, fall
#gw fall am
######
GWWA.Fall.AM<-(VAG01_fall_x2_clip + TNG13_fall_x2_clip + PAG12_fall_x2_clip + VAG02_fall_x2_clip + PA05_fall_x2_clip + PA11_fall_x2_clip + TN05_fall_x2_clip + TN06_fall_x2_clip + TN09_fall_x2_clip +TN10_fall_x2_clip + TN13_2015_fall_x2_clip + TN13_fall_x2_clip + TN16_fall_x2_clip)
GWWA.Fall.AM.agg<-aggregate(GWWA.Fall.AM, fact=1.5)
GWWA.Fall.AM.agg<-GWWA.Fall.AM.agg/cellStats(GWWA.Fall.AM.agg, sum)
######
#png("GW_AM_Fall.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(GWWA.Fall.AM.agg, col=cramp1(50), main="GWWA AM Fall")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(gwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(gwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(AM_BCR, add=TRUE, border=rgb(255/255, 157/255, 9/255, 1), col=rgb(255/255, 157/255, 9/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(ngwFall))>0 & lat<43 & lat>25), pch=21, bg=rgb(255/255, 157/255, 9/255, 1), col="black")
#points(subset(sites, as.numeric(as.character(ngwFall))>0 & lat<20), pch=21, col="black", cex=as.numeric(as.character(gwcexlist@data[c(14:19),15])))
#points(subset(sites, as.numeric(as.character(ngwFall))>0 & lon> -80), pch=21, bg=c(gwtrendcolList[c(4,8,9,10)]), col="black", cex=as.numeric(as.character(gwcexlist@data[c(4,8,9,10),15])))
#points(subset(winterpts, as.character(pop)=="AM" & fall==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="AM"), lwd=1, cex=1)
#dev.off()Golden-winged warbler, BHT, spring
#GW GL SPRING
#####
GWWA.Spring.GL<-(MIG01_spring_x2_clip + RLG12_spring_x2_clip + RLG18_spring_x2_clip + ONG03_spring_x2_clip + ONG05_spring_x2_clip + ONG10_spring_x2_clip + RLG23_spring_x2_clip + RLGF15_spring_x2_clip + RMG02_spring_x2_clip + SHG07_spring_x2_clip + SLG12_spring_x2_clip + SLG18_spring_x2_clip + TAG07_spring_x2_clip + TAG09_spring_x2_clip + TAG14_spring_x2_clip + WIG01_spring_x2_clip + WIG02_spring_x2_clip + MN05_spring_x2_clip + MN06_spring_x2_clip + MN11_spring_x2_clip + MN12_spring_x2_clip + MN16_spring_x2_clip + MN20_spring_x2_clip + MN25_spring_x2_clip + MN29_spring_x2_clip + MN36_spring_x2_clip + B1_spring_x2_clip + B2_spring_x2_clip + B3_spring_x2_clip + B4_spring_x2_clip + B5_spring_x2_clip + B6_spring_x2_clip + B7_spring_x2_clip + B8_spring_x2_clip + B9_spring_x2_clip + B10_spring_x2_clip + B11_spring_x2_clip + B12_spring_x2_clip + B13_spring_x2_clip + B14_spring_x2_clip + B15_spring_x2_clip + B16_spring_x2_clip + B17_spring_x2_clip + B18_spring_x2_clip + B19_spring_x2_clip + B20_spring_x2_clip)
GWWA.Spring.GL.agg<-aggregate(GWWA.Spring.GL, fact=1.5)
GWWA.Spring.GL.agg<-GWWA.Spring.GL.agg/cellStats(GWWA.Spring.GL.agg, sum)
#####
#png("GW_GL_Spring.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(GWWA.Spring.GL.agg, col=cramp1(50), main="GWWA GL Spring")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(gwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(gwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
#plot(BHT_BCR, add=TRUE, border=rgb(146/255,0/255,0/255, 1), col=rgb(146/255,0/255,0/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(ngwSpr))>0 & lat>43), pch=21, bg=rgb(146/255,0/255,0/255, 1), col="black")
points(subset(sites, as.numeric(as.character(ngwSpr))>0 & lat<20), pch=21, bg=rgb(146/255,0/255,0/255, 1), col="black")
#points(subset(sites, as.numeric(as.character(ngwFall))>0 & lon> -80), pch=21, bg=c(gwtrendcolList[c(4,8,9,10)]), col="black", cex=as.numeric(as.character(gwcexlist@data[c(4,8,9,10),15])))
#points(subset(winterpts, as.character(pop)=="BHT" & spr==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="BHT"), lwd=1, cex=1)
#dev.off()Golden-winged warbler, AM, spring
#GW AM SPRING
#####
GWWA.Spring.AM<-(TNG13_spring_x2_clip + VAG01_spring_x2_clip + PAG12_spring_x2_clip + VAG02_spring_x2_clip + PA05_spring_x2_clip + PA11_spring_x2_clip + TN05_spring_x2_clip + TN06_spring_x2_clip + TN09_spring_x2_clip +TN10_spring_x2_clip + TN13_2015_spring_x2_clip + TN13_spring_x2_clip + TN16_spring_x2_clip)
GWWA.Spring.AM.agg<-aggregate(GWWA.Spring.AM, fact=1.5)
GWWA.Spring.AM.agg<-GWWA.Spring.AM.agg/cellStats(GWWA.Spring.AM.agg, sum)
#####
#png("GW_AM_Spring.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(GWWA.Spring.AM.agg, col=cramp1(50), main="GWWA AM Spring")
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
plot(gwwaBR, add=TRUE, border="gray65", lwd=1, lty=1)
plot(gwwaNB, add=TRUE, border="gray45", lwd=1, lty=1)
maps::map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
#plot(AM_BCR, add=TRUE, border=rgb(255/255, 157/255, 9/255, 1), col=rgb(255/255, 157/255, 9/255, .3), lwd=.75, lty=1)
points(subset(sites, as.numeric(as.character(ngwSpr))>0 & lat<43 & lat>25), pch=21, bg=rgb(255/255, 157/255, 9/255, 1), col="black")
#points(subset(sites, as.numeric(as.character(ngwFall))>0 & lat<20), pch=21, col="black", cex=as.numeric(as.character(gwcexlist@data[c(14:19),15])))
#points(subset(sites, as.numeric(as.character(ngwFall))>0 & lon> -80), pch=21, bg=c(gwtrendcolList[c(4,8,9,10)]), col="black", cex=as.numeric(as.character(gwcexlist@data[c(4,8,9,10),15])))
#points(subset(winterpts, as.character(pop)=="AM" & spr==1), col="gray50", pch=subset(winterpts$pch1, winterpts$pop=="AM"), lwd=1, cex=1)
#dev.off()Extract top 25th percentile areas for each population/season
#################
#TOP 5 & 10 & 25 % of each population/season overlain
##########
#BW CH FALL
pvec <- sort(getValues(BWWA.Fall.CH.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
BW.CH.FA.25 <- BWWA.Fall.CH.agg # Initialize the output
BW.CH.FA.25[BWWA.Fall.CH.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
BW.CH.FA.25=BW.CH.FA.25/cellStats(BW.CH.FA.25, sum)
plot(BW.CH.FA.25, col=rgb(255/255, 9/255, 231/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#BW CH SPRING
pvec <- sort(getValues(BWWA.Spring.CH.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
BW.CH.SP.25 <- BWWA.Spring.CH.agg # Initialize the output
BW.CH.SP.25[BWWA.Spring.CH.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
BW.CH.SP.25=BW.CH.SP.25/cellStats(BW.CH.SP.25, sum)
plot(BW.CH.SP.25, col=rgb(255/255, 9/255, 231/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#BW PHT FALL
pvec <- sort(getValues(BWWA.Fall.PHT.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
BW.PH.FA.25 <- BWWA.Fall.PHT.agg # Initialize the output
BW.PH.FA.25[BWWA.Fall.PHT.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
BW.PH.FA.25=BW.PH.FA.25/cellStats(BW.PH.FA.25, sum)
plot(BW.PH.FA.25, col=rgb(0/255, 146/255, 146/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#BW PHT SPRING
pvec <- sort(getValues(BWWA.Spring.PHT.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
BW.PH.SP.25 <- BWWA.Spring.PHT.agg # Initialize the output
BW.PH.SP.25[BWWA.Spring.PHT.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
BW.PH.SP.25=BW.PH.SP.25/cellStats(BW.PH.SP.25, sum)
plot(BW.PH.SP.25, col=rgb(0/255, 146/255, 146/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#BW AM FALL
pvec <- sort(getValues(BWWA.Fall.AM.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
BW.AM.FA.25 <- BWWA.Fall.AM.agg # Initialize the output
BW.AM.FA.25[BWWA.Fall.AM.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
BW.AM.FA.25=BW.AM.FA.25/cellStats(BW.AM.FA.25, sum)
plot(BW.AM.FA.25, col=rgb(255/255, 157/255, 9/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#BW AM SPRING
pvec <- sort(getValues(BWWA.Spring.AM.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
BW.AM.SP.25 <- BWWA.Spring.AM.agg # Initialize the output
BW.AM.SP.25[BWWA.Spring.AM.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
BW.AM.SP.25=BW.AM.SP.25/cellStats(BW.AM.SP.25, sum)
plot(BW.AM.SP.25, col=rgb(255/255, 157/255, 9/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#GW AM FALL
pvec <- sort(getValues(GWWA.Fall.AM.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
GW.AM.FA.25 <- GWWA.Fall.AM.agg # Initialize the output
GW.AM.FA.25[GWWA.Fall.AM.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
GW.AM.FA.25=GW.AM.FA.25/cellStats(GW.AM.FA.25, sum)
plot(GW.AM.FA.25, col=rgb(255/255, 157/255, 9/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#GW AM SPRING
pvec <- sort(getValues(GWWA.Spring.AM.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
GW.AM.SP.25 <- GWWA.Spring.AM.agg # Initialize the output
GW.AM.SP.25[GWWA.Spring.AM.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
GW.AM.SP.25=GW.AM.SP.25/cellStats(GW.AM.SP.25, sum)
plot(GW.AM.SP.25, col=rgb(255/255, 157/255, 9/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#GW BHT FALL
pvec <- sort(getValues(GWWA.Fall.GL.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
GW.BH.FA.25 <- GWWA.Fall.GL.agg # Initialize the output
GW.BH.FA.25[GWWA.Fall.GL.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
GW.BH.FA.25=GW.BH.FA.25/cellStats(GW.BH.FA.25, sum)
plot(GW.BH.FA.25, col=rgb(146/255,0/255,0/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
#GW BHT SPRING
pvec <- sort(getValues(GWWA.Spring.GL.agg), decreasing=TRUE) # The probabilities, sorted
d <- cumsum(pvec) # Cumulative probabilities
dpos <- d[d <= 0.25] # Position to stop
GW.BH.SP.25 <- GWWA.Spring.GL.agg # Initialize the output
GW.BH.SP.25[GWWA.Spring.GL.agg < pvec[length(dpos)]] <- NA # Exclude the last 25% of the probability
GW.BH.SP.25=GW.BH.SP.25/cellStats(GW.BH.SP.25, sum)
plot(GW.BH.SP.25, col=rgb(146/255,0/255,0/255, .5)) # Display the result
plot(wrld_simpl,add=T, border = "grey70")
######
#individual plots of 5 10 and 20%
###########
#BWWA AGGREGATE (FALL)
#png("BW_FALL_25perc_Aggregate.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="", xlab="", main="Blue-winged warbler, fall")
raster::plot(BW.PH.FA.25, col=rgb(0/255, 146/255, 146/255, .5), add=TRUE, legend=FALSE)
raster::plot(BW.AM.FA.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
raster::plot(BW.CH.FA.25, col=rgb(255/255, 9/255, 231/255, .5), add=TRUE, legend=FALSE)
map(col = "grey70", lwd=.5, add=TRUE)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(nbwFall))>0), pch=24, col="gray30", bg=subset(sites@data$BCRcolBW, as.numeric(as.character(sites@data$nbwFall))>0), lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="BW" & fall==1), pch=4, col=subset(winterpts$color, winterpts$sp=="BW" & winterpts$fall==1), lwd=.5, cex=1)
box(cex=2)
#dev.off()
#BWWA AGGREGATE (Spring)
#png("BW_SPRING_25perc_Aggregate.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="", xlab="", main="Blue-winged warbler, spring")
plot(BW.PH.SP.25, col=rgb(0/255, 146/255, 146/255, .5), add=TRUE, legend=FALSE)
plot(BW.AM.SP.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
plot(BW.CH.SP.25, col=rgb(255/255, 9/255, 231/255, .5), add=TRUE, legend=FALSE)
plot(wrld_simpl, add=TRUE, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(nbwSpr)) > 0), pch=24, col="gray30", bg=subset(sites@data$BCRcolBW, as.numeric(as.character(sites@data$nbwSpr)) > 0), lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="BW" & spr==1), pch=4, col=subset(winterpts$color, winterpts$sp=="BW" & winterpts$spr==1), lwd=.5, cex=1)
box(cex=2)
#dev.off()
#GWWA AGGREGATE (FALL)
#png("GW_FALL_25perc_Aggregate.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="", xlab="", main="Golden-winged warbler, fall")
plot(GW.AM.FA.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
plot(GW.BH.FA.25, col=rgb(146/255,0/255,0/255, .5), add=TRUE, legend=FALSE)
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(ngwFall)) > 0), pch=24, col="gray30", bg=subset(sites@data$BCRcolGW, as.numeric(as.character(sites@data$ngwFall)) > 0), lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="GW" & fall==1), pch=4, col=subset(winterpts$color, winterpts$sp=="GW" & winterpts$fall==1), lwd=.5, cex=1)
box(cex=2)
#dev.off()
#GWWA AGGREGATE (Spring)
#png("GW_SPRING_25perc_Aggregate.png", width = 6.85039, height = 6.85039, units = 'in', res = 600)
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="", xlab="", main="Golden-winged warbler, spring")
plot(GW.AM.SP.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
plot(GW.BH.SP.25, col=rgb(146/255,0/255,0/255, .5), add=TRUE, legend=FALSE)
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(ngwSpr)) > 0 & lat>20), pch=24, col="gray30", bg=subset(sites@data$BCRcolGW, as.numeric(as.character(sites@data$ngwSpr)) > 0), lwd=.5, cex=1)
points(subset(sites, as.numeric(as.character(ngwSpr)) > 0 & lat<20), pch=25, col="gray30", bg=sites@data$BCRcolGW[1], lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="GW" & spr==1), pch=4, col=subset(winterpts$color, winterpts$sp=="GW" & winterpts$spr==1), lwd=.5, cex=1)
box(cex=2)
#dev.off()###25 %ile shapefile conversion
BW.PH.FA.25.x<-BW.PH.FA.25 > -Inf
BW.PH.FA.25.shp<-rasterToPolygons(BW.PH.FA.25.x, dissolve=TRUE)
BW.PH.SP.25.x<-BW.PH.SP.25 > -Inf
BW.PH.SP.25.shp<-rasterToPolygons(BW.PH.SP.25.x, dissolve=TRUE)
BW.CH.FA.25.x<-BW.CH.FA.25 > -Inf
BW.CH.FA.25.shp<-rasterToPolygons(BW.CH.FA.25.x, dissolve=TRUE)
BW.CH.SP.25.x<-BW.CH.SP.25 > -Inf
BW.CH.SP.25.shp<-rasterToPolygons(BW.CH.SP.25.x, dissolve=TRUE)
BW.AM.FA.25.x<-BW.AM.FA.25 > -Inf
BW.AM.FA.25.shp<-rasterToPolygons(BW.AM.FA.25.x, dissolve=TRUE)
BW.AM.SP.25.x<-BW.AM.SP.25 > -Inf
BW.AM.SP.25.shp<-rasterToPolygons(BW.AM.SP.25.x, dissolve=TRUE)
GW.BH.FA.25.x<-GW.BH.FA.25 > -Inf
GW.BH.FA.25.shp<-rasterToPolygons(GW.BH.FA.25.x, dissolve=TRUE)
GW.BH.SP.25.x<-GW.BH.SP.25 > -Inf
GW.BH.SP.25.shp<-rasterToPolygons(GW.BH.SP.25.x, dissolve=TRUE)
GW.AM.FA.25.x<-GW.AM.FA.25 > -Inf
GW.AM.FA.25.shp<-rasterToPolygons(GW.AM.FA.25.x, dissolve=TRUE)
GW.AM.SP.25.x<-GW.AM.SP.25 > -Inf
GW.AM.SP.25.shp<-rasterToPolygons(GW.AM.SP.25.x, dissolve=TRUE)
###plot
####Putting it all together
#png("25perc_Aggregate_COMBINED_with75perc.png", width = (6.85039), height = (6.85039), units = 'in', res = 600)
par(mfrow=c(2,2), pty="s")
#bw fall
par(mar=c(2,5,3,0))
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="Blue-wined warbler", xlab="", main="Autumn", cex.lab=1.5, font.lab=2, cex.main=1.5)
raster::plot(BW.PH.FA.25, col=rgb(0/255, 146/255, 146/255, .5), add=TRUE, legend=FALSE)
raster::plot(BW.AM.FA.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
raster::plot(BW.CH.FA.25, col=rgb(255/255, 9/255, 231/255, .5), add=TRUE, legend=FALSE)
# raster::plot(BW.PH.FA.75.shp, border=rgb(0/255, 146/255, 146/255, .2), add=TRUE, lwd=1, lty=1)
# raster::plot(BW.CH.FA.75.shp, border=rgb(255/255, 9/255, 231/255, .2), add=TRUE, lwd=1, lty=1)
# raster::plot(BW.AM.FA.75.shp, border=rgb(255/255, 157/255, 9/255, .2), add=TRUE, lwd=1, lty=1)
map(col = "grey70", lwd=.5, add=TRUE)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(nbwFall))>0), pch=24, col="gray30", bg=subset(sites@data$BCRcolBW, as.numeric(as.character(sites@data$nbwFall))>0), lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="BW" & fall==1), pch=4, col=subset(winterpts$color, winterpts$sp=="BW" & winterpts$fall==1), lwd=.5, cex=1)
legend(x=-120, y=10, bty='n', legend=c("Prairie-Hardwoods Transition BCR", "Central Hardwoods BCR", "Appalachian Mountain BCR"), pch=15, col=c(rgb(0/255, 146/255, 146/255, .5),rgb(255/255, 9/255, 231/255, .5),rgb(255/255, 157/255, 9/255, .5)), cex=.75, pt.cex=1.5)
box(cex=2)
#bw spring
par(mar=c(2,2,3,3))
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="", xlab="", yaxt='n', xaxt='n', main= "Spring", cex.lab=1.5, font.lab=2, cex.main=1.5)
plot(BW.PH.SP.25, col=rgb(0/255, 146/255, 146/255, .5), add=TRUE, legend=FALSE)
plot(BW.AM.SP.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
plot(BW.CH.SP.25, col=rgb(255/255, 9/255, 231/255, .5), add=TRUE, legend=FALSE)
# raster::plot(BW.PH.SP.75.shp, border=rgb(0/255, 146/255, 146/255, .2), add=TRUE, lwd=1, lty=1)
# raster::plot(BW.CH.SP.75.shp, border=rgb(255/255, 9/255, 231/255, .2), add=TRUE, lwd=1, lty=1)
# raster::plot(BW.AM.SP.75.shp, border=rgb(255/255, 157/255, 9/255, .2), add=TRUE, lwd=1, lty=1)
plot(wrld_simpl, add=TRUE, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(nbwSpr)) > 0), pch=24, col="gray30", bg=subset(sites@data$BCRcolBW, as.numeric(as.character(sites@data$nbwSpr)) > 0), lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="BW" & spr==1), pch=4, col=subset(winterpts$color, winterpts$sp=="BW" & winterpts$spr==1), lwd=.5, cex=1)
legend(x=-120, y=7, bty='n', legend=c("Geolocator deployment site", "Geolocator-derived nonbreeding site estimate"), pch=c(2,4), col="gray30", cex=.75, pt.cex=.75)
box(cex=2)
#gw fall
par(mar=c(4,5,1,0))
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="Golden-winged warbler", xlab="", yaxt='n', xaxt='n', cex.lab=1.5, font.lab=2, cex.main=1.5)
plot(GW.AM.FA.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
plot(GW.BH.FA.25, col=rgb(146/255,0/255,0/255, .5), add=TRUE, legend=FALSE)
# raster::plot(GW.BH.FA.75.shp, border=rgb(146/255,0/255,0/255, .2), add=TRUE, lwd=1, lty=1)
# raster::plot(GW.AM.FA.75.shp, border=rgb(255/255, 157/255, 9/255, .2), add=TRUE, lwd=1, lty=1)
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(ngwFall)) > 0 & lat>20), pch=24, col="gray30", bg=subset(sites@data$BCRcolGW, as.numeric(as.character(sites@data$ngwFall)) > 0), lwd=.5, cex=1)
points(subset(sites, as.numeric(as.character(ngwFall)) > 0 & lat<20), pch=24, col="gray30", bg=rgb(146/255,0/255,0/255, 1), lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="GW" & fall==1), pch=subset(winterpts@data$pch1, winterpts@data$sp=="GW" & winterpts@data$fall==1), col=subset(winterpts$color, winterpts$sp=="GW" & winterpts$fall==1), lwd=.5, cex=1)
legend(x=-120, y=7, bty='n', legend=c("Boreal-Hardwoods Transition BCR", "Appalachian Mountain BCR"), pch=15, col=c(rgb(146/255,0/255,0/255, .5), rgb(255/255, 157/255, 9/255, .5)), cex=.75, pt.cex=1.5)
box(cex=2)
#gw spring
par(mar=c(4,2,1,3))
plot(1,1, type='n', xlim=c(-120,-60), ylim=c(0,60), ylab="", xlab="", yaxt='n', xaxt='n')
plot(GW.AM.SP.25, col=rgb(255/255, 157/255, 9/255, .5), add=TRUE, legend=FALSE)
plot(GW.BH.SP.25, col=rgb(146/255,0/255,0/255, .5), add=TRUE, legend=FALSE)
# raster::plot(GW.BH.SP.75.shp, border=rgb(146/255,0/255,0/255, .2), add=TRUE, lwd=1, lty=1)
# raster::plot(GW.AM.SP.75.shp, border=rgb(255/255, 157/255, 9/255, .2), add=TRUE, lwd=1, lty=1)
plot(wrld_simpl, add=T, border = "grey70", lwd=.5)
map('state', add=TRUE, col = "grey70", lwd=.5)
plot(canada, add=TRUE, border="grey70", lwd=.5)
points(subset(sites, as.numeric(as.character(ngwSpr)) > 0 & lat>20), pch=24, col="gray30", bg=subset(sites@data$BCRcolGW, as.numeric(as.character(sites@data$ngwSpr)) > 0), lwd=.5, cex=1)
points(subset(sites, as.numeric(as.character(ngwSpr)) > 0 & lat<20), pch=24, col="gray30", bg=rgb(146/255,0/255,0/255, 1), lwd=.5, cex=1)
points(subset(winterpts, as.character(sp)=="GW" & spr==1), pch=subset(winterpts@data$pch1, winterpts@data$sp=="GW" & winterpts@data$spr==1), col=subset(winterpts$color, winterpts$sp=="GW" & winterpts$spr==1), lwd=.5, cex=1)
legend(x=-120, y=5, bty='n', legend=c("Geolocator-derived breeding site estimate"), pch=c(3), col="gray30", cex=.75, pt.cex=.75)
box(cex=2)
#dev.off()##### OVERLAP CALCULATION
BW.PH.FA.25.POLY<-rasterToPolygons(BW.PH.FA.25, dissolve=TRUE)
BW.PH.SP.25.POLY<-rasterToPolygons(BW.PH.SP.25, dissolve=TRUE)
BW.CH.FA.25.POLY<-rasterToPolygons(BW.CH.FA.25, dissolve=TRUE)
BW.CH.SP.25.POLY<-rasterToPolygons(BW.CH.SP.25, dissolve=TRUE)
BW.AM.FA.25.POLY<-rasterToPolygons(BW.AM.FA.25, dissolve=TRUE)
BW.AM.SP.25.POLY<-rasterToPolygons(BW.AM.SP.25, dissolve=TRUE)
GW.AM.FA.25.POLY<-rasterToPolygons(GW.AM.FA.25, dissolve=TRUE)
GW.AM.SP.25.POLY<-rasterToPolygons(GW.AM.SP.25, dissolve=TRUE)
GW.BH.FA.25.POLY<-rasterToPolygons(GW.BH.FA.25, dissolve=TRUE)
GW.BH.SP.25.POLY<-rasterToPolygons(GW.BH.SP.25, dissolve=TRUE)
areaPHfall<-sum(area(BW.PH.FA.25.POLY)/1000000)
areaPHspr<-sum(area(BW.PH.SP.25.POLY)/1000000)
areaCHfall<-sum(area(BW.CH.FA.25.POLY)/1000000)
areaCHspr<-sum(area(BW.CH.SP.25.POLY)/1000000)
areabwAMfall<-sum(area(BW.AM.FA.25.POLY)/1000000)
areabwAMspr<-sum(area(BW.AM.SP.25.POLY)/1000000)
areagwAMfall<-sum(area(GW.AM.FA.25.POLY)/1000000)
areagwAMspr<-sum(area(GW.AM.SP.25.POLY)/1000000)
areaBHfall<-sum(area(GW.BH.FA.25.POLY)/1000000)
areaBHspr<-sum(area(GW.BH.SP.25.POLY)/1000000)
#intersect between populaitons
#bw
int.ph.ch.fall<-intersect(BW.PH.FA.25.POLY,BW.CH.FA.25.POLY)
int.ph.bwam.fall<-intersect(BW.PH.FA.25.POLY,BW.AM.FA.25.POLY)
int.bwam.ch.fall<-intersect(BW.AM.FA.25.POLY,BW.CH.FA.25.POLY)
areaPH.CHfall<-sum(area(int.ph.ch.fall)/1000000)
areaPH.bwAMfall<-sum(area(int.ph.bwam.fall)/1000000)
areabwAM.CHfall<-sum(area(int.bwam.ch.fall)/1000000)
int.ph.ch.spr<-intersect(BW.PH.SP.25.POLY,BW.CH.SP.25.POLY)
int.ph.bwam.spr<-intersect(BW.PH.SP.25.POLY,BW.AM.SP.25.POLY)
int.bwam.ch.spr<-intersect(BW.AM.SP.25.POLY,BW.CH.SP.25.POLY)
areaPH.CHspr<-sum(area(int.ph.ch.spr)/1000000)
areaPH.bwAMspr<-sum(area(int.ph.bwam.spr)/1000000)
areabwAM.CHspr<-sum(area(int.bwam.ch.spr)/1000000)
#gw
int.bh.gwam.fall<-intersect(GW.BH.FA.25.POLY,GW.AM.FA.25.POLY)
#int.bh.gwam.spr<-intersect(GW.BH.SP.25.POLY,GW.AM.SP.25.POLY) #No overlap
areaBH.gwAMfall<-sum(area(int.bh.gwam.fall)/1000000)
#areaBH.gwAMspr<-sum(area(int.bh.gwam.spr)/1000000) #No overlap
#all verm combos
int.ph.gwam.fall<-intersect(BW.PH.FA.25.POLY,GW.AM.FA.25.POLY)
int.ph.bh.fall<-intersect(BW.PH.FA.25.POLY,GW.BH.FA.25.POLY)
int.ch.gwam.fall<-intersect(BW.CH.FA.25.POLY,GW.AM.FA.25.POLY)
int.ch.bh.fall<-intersect(BW.CH.FA.25.POLY,GW.BH.FA.25.POLY)
int.bwam.gwam.fall<-intersect(BW.AM.FA.25.POLY,GW.AM.FA.25.POLY)
int.bwam.bh.fall<-intersect(BW.AM.FA.25.POLY,GW.BH.FA.25.POLY)
areaPH.gwAMfall<-sum(area(int.ph.gwam.fall)/1000000)
areaPH.BHfall<-sum(area(int.ph.bh.fall)/1000000)
areaCH.gwAMfall<-sum(area(int.ch.gwam.fall)/1000000)
areaCH.BHfall<-sum(area(int.ch.bh.fall)/1000000)
areabwAM.gwAMfall<-sum(area(int.bwam.gwam.fall)/1000000)
areabwAM.BHfall<-sum(area(int.bwam.bh.fall)/1000000)
int.ph.gwam.spr<-intersect(BW.PH.SP.25.POLY,GW.AM.SP.25.POLY)
int.ph.bh.spr<-intersect(BW.PH.SP.25.POLY,GW.BH.SP.25.POLY)
int.ch.gwam.spr<-intersect(BW.CH.SP.25.POLY,GW.AM.SP.25.POLY)
int.ch.bh.spr<-intersect(BW.CH.SP.25.POLY,GW.BH.SP.25.POLY)
#int.bwam.gwam.spr<-intersect(BW.AM.SP.25.POLY,GW.AM.SP.25.POLY) #No overlap
int.bwam.bh.spr<-intersect(BW.AM.SP.25.POLY,GW.BH.SP.25.POLY)
areaPH.gwAMspr<-sum(area(int.ph.gwam.spr)/1000000)
areaPH.BHspr<-sum(area(int.ph.bh.spr)/1000000)
areaCH.gwAMspr<-sum(area(int.ch.gwam.spr)/1000000)
areaCH.BHspr<-sum(area(int.ch.bh.spr)/1000000)
#areabwAM.gwAMspr<-sum(area(int.bwam.gwam.spr)/1000000) #NO OVERLAP
areabwAM.BHspr<-sum(area(int.bwam.bh.spr)/1000000)
#intersect same population -- different seasons
int.ph.fall.spr<-intersect(BW.PH.FA.25.POLY,BW.PH.SP.25.POLY)
int.ch.fall.spr<-intersect(BW.CH.FA.25.POLY,BW.CH.SP.25.POLY)
int.bwam.fall.spr<-intersect(BW.AM.FA.25.POLY,BW.AM.SP.25.POLY)
int.bh.fall.spr<-intersect(GW.BH.FA.25.POLY,GW.BH.SP.25.POLY)
int.gwam.fall.spr<-intersect(GW.AM.FA.25.POLY,GW.AM.SP.25.POLY)
areaPH.FS<-sum(area(int.ph.fall.spr)/1000000)
areaCH.FS<-sum(area(int.ch.fall.spr)/1000000)
areabwAM.FS<-sum(area(int.bwam.fall.spr)/1000000)
areaBH.FS<-sum(area(int.bh.fall.spr)/1000000)
areagwAM.FS<-sum(area(int.gwam.fall.spr)/1000000)