Occupancy and Abundance of Eastern Whip-poor-will in USFWS Northeast Region

Overview

Estimation of probability of site occupancy and abundance for eastern whip-poor-will (Antrostomus vociferus) in the U.S. Fish and Wildlife (USFWS) Northeast Region to support the management planning and guidance developed by the At-Risk Pine Barrens Working group.

Purpose

The purpose of this modeling effort is to identify key areas of high conservation importance for eastern whip-poor-will and where these areas are in relation to the imperriled frosted elfin (Callophrys irus)

Objectives

1. Estimate species occupancy
   
2. Estimate species abundance
   
3. Build predictive models using remotely sensed information (e.g., Landfire vegetation data) that influence the occupancy and abundance for the focal species
   
4. Identify key areas of imporance for species

Methods

Using citizen-science based eBird data (https://ebird.org/home) we estimated the probability of site occupancy and relative abundance for the eastern whip-poor-will throughout the USFWS Northeast Region. These data were downloaded from eBird checklists where observers detected whip-poor-wills during the breeding season in June, July, and August annualy from 2016-2021. These data were then sampled to generate a spatially-balanced data set, and analyzed using hierarchical N-mixture models using the ‘unmarked’ package in program R (Core Team R 2022). Combined with Landfire data consisiting of 12 unique ecosystem types and elevation, we then modeled probability of site occupancy and abundance while accounting for imperfect detection. Finally, we predicted (mean and SE) the spatially-explicit probability of whip-poor-will occupancy and abundance across a 5-km^2 resolution surface of the Northeast region to aid in identifying key areas where whip-poor-will occur (and do not occur), to help inform next steps for the Pine Barrens Habitat At-risk working group’s action plan.

Load required packages

#clear environment
rm(list=ls())

#load packages
suppressMessages({
library(knitr)
library(readxl)
library(unmarked)
library(ggplot2)
library(ggrepel)
library(ggmap)
library(lubridate)
library(data.table)
library(reshape2)
library(plyr)
library(kableExtra)
library(pander)
library(raster)
library(pals)
library(rgdal)
library(gdalUtils)
library(sp)
library(lattice)
library(dplyr)
library(rgeos)
library(exactextractr)
library(purrr)
library(tidyr)
library(dismo)
library(ggsn)
library(MuMIn)
library(rebird)
library(sf)
library(dggridR)

# resolve namespace conflicts
select <- dplyr::select
map <- purrr::map
projection <- raster::projection
  
#set working directory
opts_knit$set(root.dir = "/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/")
#setwd("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/")

#register Google Maps API key
register_google(key="AIzaSyBLb0jdTnIwfWzduauW5sjfxE91hnqRvfE")

#load functions
source("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Source/makePOccu_umf.R")
source("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Source/makePOccu_umf_with_covs.R")
source("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Source/makePCount_umf.R")
source("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Source/makePCount_umf_with_covs.R")
source("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Source/AICcmodavg.gof.test_source.R")
source("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Source/xyToLatLong.R")

})

Download eBird detection data using ‘rebird’ package from Jun, Jul, and August 2017-2021

#define API key from https://ebird.org/api/keygen
myKey<-"u9drtus3ammg"

#example query for historic data
test.data<-ebirdhistorical(loc ="US-ME",date='2019-06-15', fieldSet='full', key=myKey)

#iterate through dates in June to get detection data

#first create list of dates (last 5 years)
dates2017<-seq.Date(from=as.Date("2017-06-01"), to=as.Date("2017-09-1"), by=1)
dates2018<-seq.Date(from=as.Date("2018-06-01"), to=as.Date("2018-09-1"), by=1)
dates2019<-seq.Date(from=as.Date("2019-06-01"), to=as.Date("2019-09-1"), by=1)
dates2020<-seq.Date(from=as.Date("2020-06-01"), to=as.Date("2020-09-1"), by=1)
dates2021<-seq.Date(from=as.Date("2021-06-01"), to=as.Date("2021-09-1"), by=1)

#combine into a single list
datesList<-as.character(c(dates2017, dates2018, dates2019, dates2020, dates2021))

#get list of locations
locList<-c("US-CT","US-DE","US-ME","US-MD","US-MA","US-NH","US-NJ","US-NY","US-PA","US-RI","US-VT","US-VA","US-WV")

#empty list to capture data
save.data<-list()

#for loop
for(i in 1:length(locList)){
  
  #new location
  new.loc<-locList[i]
  print(new.loc)
  
  #empty list to save date data
  date.data<-list()
  
  for(j in 1:length(datesList)){
    
    #new date
    new.date<-datesList[j]
    print(new.date)
    
    #get new.data
    new.data<-ebirdhistorical(loc = new.loc ,date=new.date, fieldSet='full', key=myKey)
    
    date.data<-rbind(date.data, new.data)
  }
  
  #now row bind location data
  save.data<-rbind(save.data,date.data)
}

#save data
write.csv(save.data, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/eBird_USFWS_LR5_data.csv",row.names = FALSE)

Get state polygons for Northeast Region

#get state polygons
stateMap <- raster::getData('GADM', country='USA', level=2)

#list of states
statesList<-c("Connecticut","Delaware","Maine","Maryland","Massachusetts","New Hampshire","New Jersey","New York","Pennsylvania","Rhode Island","Vermont","Virginia","West Virginia")

#subset to only include desired states
states.sub<-stateMap[stateMap@data$NAME_1 %in% statesList,]

states.poly<-fortify(states.sub, region="NAME_1")
states.poly <- states.poly[order(states.poly$order),]
states.poly$id<-as.factor(states.poly$id)

Subset data to include only those points within state polygons.

#read in save.data
save.data<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/eBird_USFWS_LR5_data.csv",header=TRUE)

#remove unverified checklists
save.data.sub<-subset(save.data, obsReviewed==TRUE)

#keep all checklists
save.data.sub<-save.data

# replace any non-whip-poor-will data with 0s
save.data.sub$howMany<-ifelse(save.data.sub$speciesCode=="easwpw1",save.data.sub$howMany,0)

#make column for WPWI_Present a factor
save.data.sub$WPWI_Present<-ifelse(save.data.sub$howMany>0,"WPWI","NO_WPWI")

#levels(save.data.simp.2$WPWI_Present)

#filter checklists where detections occur within CONUS (remove any lists from marine/offshore)

#get only points that intersect state polygons
myPoints<-st_as_sf(save.data.sub, coords = c("lng", "lat"), crs = 4326)

myPolys<-st_as_sf(states.sub)
#plot(myPolys)

#get points intersecting polygons
pts_in<-st_intersects(myPoints, myPolys)

#length(pts_in)

#grab all rows for points that intersect the state polygons
pts_in.df<-list()
for(i in 1:length(pts_in)){
  #print(i)
  new.value<-ifelse(length(unlist(pts_in[i]))==0,0,unlist(pts_in[i]))
  new.data<-data.frame("row.id"= i ,"num_pts"=new.value)
  pts_in.df<-rbind(pts_in.df, new.data)
}

#subset to remove any rows with num_pts==0
pts_in.df.sub<-subset(pts_in.df, num_pts!=0)

#make row.id a character
pts_in.df.sub$row.id<-as.character(pts_in.df.sub$row.id)

#merge with save.data.simp

#add row names as row.id field in save.data.simp
save.data.sub$row.id<-seq(1, nrow(save.data.sub), by=1)

save.data.merge<-merge(save.data.sub, pts_in.df.sub, by="row.id",all.y=TRUE)

#remove NAs
save.data.merge<-na.omit(save.data.merge)

#export save.data.merge
write.csv(save.data.merge, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/eBird_USFWS_LR5_data_filter.csv",row.names=FALSE)

Genearate basemap for figures

save.data.merge<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/eBird_USFWS_LR5_data_filter.csv")

#define map projection
map_proj <- "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

#get bounding box coords from eBird data
x.min=min(save.data.merge$lng, na.rm=TRUE)-0.4
x.max=max(save.data.merge$lng, na.rm=TRUE)+0.4
y.min=min(save.data.merge$lat, na.rm=TRUE)-0.3
y.max=max(save.data.merge$lat, na.rm=TRUE)+0.3

#compile bounding box corner coordinates
ebird.bbox = c(left = x.min, bottom = y.min, right = x.max, top = y.max)

#define map area and type
basemap <- get_map(location =ebird.bbox, maptype="toner-lite", zoom=8)

## maptype = "toner-lite" is only available with source = "stamen".

## resetting to source = "stamen"...

## Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

## 156 tiles needed, this may take a while (try a smaller zoom).

#ggmap(basemap)

# #get map extent
mapExtent<-attr(basemap, "bb")

Figure 1. Spatial sampling of data

save.data.merge<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/eBird_USFWS_LR5_data_filter.csv")

#add Year and Month and Week columns
save.data.merge$Year<-year(as.Date(save.data.merge$obsDt, format="%Y-%m-%d"))
save.data.merge$Month<-month(as.Date(save.data.merge$obsDt, format="%Y-%m-%d"))
save.data.merge$Week<-week(as.Date(save.data.merge$obsDt, format="%Y-%m-%d"))

#iterate over years and months to get repeated 3 repeated visits per year
yearList<-sort(unique(save.data.merge$Year))
monthList<-sort(unique(save.data.merge$Month))

#subset out WPWI detections
wpwi.data<-subset(save.data.merge, howMany>0)

non.wpwi.data<-subset(save.data.merge, howMany<1)


#empty list to save spatially-sampled data (non-WPWI)
save.spatial.data<-list()

for(i in 1:length(yearList)){
  
  print(yearList[i])
  
  count.data.year<-subset(non.wpwi.data, Year==yearList[i])

  save.spatial.month<-list()
  for(j in 1:length(monthList)){
    
    print(monthList[j])
    
    count.data.month<-subset(count.data.year, Month==monthList[j])
    
    #get weekList
    weekList<-unique(count.data.month$Week)
    
    save.spatial.week<-list()
    for(z in 1:length(weekList)){
      
      count.data.week<-subset(count.data.month, Week==weekList[z])
  
      # generate hexagonal grid with ~ 35 km betweeen cells
      dggs <- dgconstruct(spacing = 50)
      
      # get hexagonal cell id for each site
      occ_wide_cell <- count.data.week %>% 
        mutate(cell = dgGEO_to_SEQNUM(dggs, lng, lat)$seqnum)
      
      # sample one site per grid cell
      occ_ss <- occ_wide_cell %>%
        group_by(cell) %>%
        sample_n(size = 1) %>%
        ungroup() %>%
        select(-cell)
      # calculate the percent decrease in the number of sites
      #1 - nrow(occ_ss) / nrow(occ_wide)
      
      #add Hexagon back in 
      occ_hex<-merge(occ_ss, occ_wide_cell[,c("row.id","cell")],by="row.id",all.x=TRUE)
    
      #compile
      save.spatial.week<-rbind(save.spatial.week, occ_hex)
    }
  #compile
    save.spatial.month<-rbind(save.spatial.month, save.spatial.week)
  }   
  #compile
save.spatial.data<-rbind(save.spatial.data, save.spatial.month)
}


# get hexagonal cell id for wpwi
wpwi.cell <- wpwi.data %>% 
        mutate(cell = dgGEO_to_SEQNUM(dggs, lng, lat)$seqnum)

#add wpwi.cell to sampled data
save.spatial.data.all<-rbind(save.spatial.data,wpwi.cell)

#count how many detections and non-detections in full data and sampled data
sampledDetections<-subset(save.spatial.data.all, howMany>0)
nrow(sampledDetections)

sampledDetectionsUnique<-unique(sampledDetections)

fullDetections<-subset(save.data.merge, howMany>0)
nrow(fullDetections)

    
#define Points column
save.data.merge$Points<-paste("All (n = ",length(unique(save.data.merge$obsId)), ")",sep="")

save.spatial.data.all$Points<-paste("Sampled (n = ", length(unique(save.spatial.data.all$obsId)), ")", sep="")

#remove cell column to produce comparison figure
save.spatial.data.1<-save.spatial.data.all
save.spatial.data.1$cell<-NULL
occ_comb<-rbind(save.data.merge, save.spatial.data.1)
occ_comb$Points<-as.factor(occ_comb$Points)


plotColors<-c("gray20","blue")

subSamplePlot<-ggmap(basemap, darken=c(0.7,"white"))+
  geom_point(data=occ_comb, aes(x=lng,y=lat,color=Points),size=0.2,alpha=0.7, shape=16, inherit.aes = FALSE)+
  geom_sf(data = hexagons,inherit.aes = FALSE,color="black",alpha=0.5,size=0.3)+
  theme(panel.border=element_rect(color="black",fill="transparent"), panel.background = element_rect(fill="white"))+
  scale_color_manual(values=plotColors)+
  labs(x="Longitude",y="Latitude")+
  guides(alpha="none")
  
subSampleFacet<-subSamplePlot+facet_wrap(.~Points)
subSampleFacet

ggsave(subSampleFacet, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_1_WPWI_spatial_sampling_eBird.png",
       width=9, height=6, dpi=600)

#Load the saved image of raster veg layer
include_graphics("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_1_WPWI_spatial_sampling_eBird.png")

Figure 2. Map of eBird checklists with detections and non-detections.

###############################################################################
#make a map of presence absence data


#plot sampling locations over study area basemap
studyAreaMap<-ggmap(basemap)+
  geom_map(data=states.poly, map=states.poly, aes(y=lat, x=long, map_id=id),color="black",fill="gray90",alpha=0.4, inherit.aes = FALSE)+
  geom_point(data=save.spatial.data.all, aes(x=lng, y=lat, color=WPWI_Present), shape=16,size=1)+
  scale_color_manual(values=c(alpha("gray40",0.3),"red"), name="Detected")+
  labs(x="Longitude",y="Latitude")+
  theme(axis.text=element_text(size=12),
    axis.title=element_text(size=14),
    panel.border=element_rect(fill="transparent", color="black"),
    legend.position = "right")+
  ggtitle("Map of detections and non-detections")+
  #geom_text_repel(check_overlap = FALSE, size=3, nudge.y=-10, nudge.x=-10,aes(x=Longitude, y=Latitude,color=PatchID, label=PatchID),data=PatchID_labels)+
  #guides(fill = guide_legend(title = "National Wildlife Refuge", label = FALSE))+
  ggsn::scalebar(x.min=mapExtent$ll.lon, x.max=mapExtent$ur.lon, y.min=mapExtent$ll.lat, y.max=mapExtent$ur.lat,transform=TRUE, dist_unit="km",dist=200, st.dist=0.04, height = 0.015, anchor=c(x=-68,y=38.6), border.size=0.5,st.size=3)

studyAreaMap

ggsave(studyAreaMap, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_2_WPWI_presence_absence_eBird.png",width=9, height=9, dpi=600)

#Load the saved image of raster veg layer
include_graphics("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_2_WPWI_presence_absence_eBird.png")

Format data

#add hexagonal cell id for each site
save.spatial.data_hex <- save.spatial.data.all
names(save.spatial.data_hex)[names(save.spatial.data_hex)=="cell"]<-"Hexagon"

#Set up dataframe with repeated visits

#create Date column
save.spatial.data_hex$Date<-as.Date(save.spatial.data_hex$obsDt, format = "%Y-%m-%d")

#add Time column
save.spatial.data_hex$obsDt<-as.POSIXct(save.spatial.data_hex$obsDt, format = "%Y-%m-%d %H:%M")
save.spatial.data_hex$Time<-format(save.spatial.data_hex$obsDt, format = "%H:%M")

#now aggregate over weeks and get max counts at each hexagon
data.agg<-aggregate(howMany~Year+Month+Time+Hexagon+lat+lng, FUN=max, data=save.spatial.data_hex)

#make Hexagon a factor
data.agg$Hexagon<-as.factor(as.character(data.agg$Hexagon))

#cast data to put Weeks in columns
data.agg.cast<-reshape2::dcast(data.agg, Hexagon~Month,value.var="howMany",fun.aggregate = max )

#subset to get checklists with at least 3 repeated visits
data.agg.cast.sub<-subset(data.agg.cast, c(`6`!="-Inf" & `7` !="-Inf" & `8`!="-Inf"  ))

#add in Hexagon Latitude and Longitudes
#cast data to put Weeks in columns
hex.lat.cast<-stats::aggregate(lat~Hexagon, FUN=mean, data=data.agg)
hex.lng.cast<-stats::aggregate(lng~Hexagon, FUN=mean, data=data.agg)
#combine
hex.coords<-merge(hex.lat.cast, hex.lng.cast, by="Hexagon",all.x=TRUE)

#merge with count data
data.agg.merge<-merge(data.agg.cast.sub, hex.coords, by="Hexagon",all.x=TRUE)

#save data
write.csv(data.agg.merge, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/WPWI_data_5years_wide.csv",row.names=FALSE)

Results

Table 1. Estimated probability of species occupancy (with simple model).

count.data<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/WPWI_data_5years_wide.csv")
#rename columns
names(count.data)<-c("Hexagon","Jun","Jul","Aug","Latitude","Longitude")

#add visit columns
count.data$Visit_1<-"Visit_1"
count.data$Visit_2<-"Visit_2"
count.data$Visit_3<-"Visit_3"

#convert data to an unmarked frame

#define counts
y=data.frame(y.1=count.data$Jun, y.2=count.data$Jul, y.3=count.data$Aug)
  
#step to convert all y values > 1 to 1
y$y.1<-ifelse(y$y.1 > 1, 1, y$y.1)
y$y.2<-ifelse(y$y.2 > 1, 1, y$y.2)
y$y.3<-ifelse(y$y.3 > 1, 1, y$y.3)

#define site covariates - FOR MULTIPLE YEARS
site.covs<- data.frame(PointID=as.factor(count.data$Hexagon),
                       Latitude=as.numeric(count.data$Latitude),
                       Longitude=as.numeric(count.data$Longitude))

#define observation covariates
obs.covs<-list(Visit=count.data[,c("Visit_1","Visit_2","Visit_3")])

#set up unmarked FRAME for PCount
occu.umf<-unmarkedFrameOccu(y=y,siteCovs = site.covs,obsCovs= obs.covs)

#fit example model to test
occu.mod<-NULL
occu.mod<-occu(~Visit ~1, data=occu.umf)
  
suppressMessages({

#run GOF test
mod.gof<- try(mb.gof.test(occu.mod, nsim = 10, plot.hist = FALSE, report = NULL, parallel = FALSE))
  
#get c-hat (goodness of fit test)
mod.c.hat<-mod.gof$c.hat.est

#get predicted estimates
occu.mod.predict<-unique(predict(occu.mod, type="state",appendData=TRUE))
})

#remove columns
occu.mod.predict<-unique(occu.mod.predict[,c("Predicted","SE")])

#remove row.names
row.names(occu.mod.predict)<-NULL

#add GOF statisitc
occu.mod.predict$GOF_c_hat<- mod.c.hat

colnames(occu.mod.predict)<-c("Mean p(Occupancy)","SE p(Occupancy)","GOF (c-hat)")
  

panderOptions('table.alignment.default', function(df) ifelse(sapply(df, is.numeric), 'right', 'left'))
panderOptions("table.split.table", Inf)

pander(occu.mod.predict)

Mean p(Occupancy)	SE p(Occupancy)	GOF (c-hat)
0.4757	0.03312	1.735

Table 2. Estimated species abundance (with simple model).

#define counts
y=data.frame(y.1=count.data$Jun, y.2=count.data$Jul, y.3=count.data$Aug)

abun.umf<-unmarkedFramePCount(y=y, 
                               siteCovs = site.covs,
                               obsCovs= obs.covs)


#fit example model to test
abun.mod<-NULL
abun.mod<-pcount(~Visit ~1, data=abun.umf)
  
  suppressMessages({
  
  #run GOF test
  mod.gof<- Nmix.gof.test(abun.mod, nsim = 20, plot.hist = FALSE,report = NULL, parallel = FALSE)
  
  #get c-hat
  mod.c.hat<-mod.gof$c.hat.est
  #set to NA for now, as it prints message otherwise (need to suppress somehow)
  #mod.c.hat<-NA
  
  #get predicted estimates

  abun.mod.predict<-unique(predict(abun.mod, type="state",appendData=TRUE))
  
  })
  
#remove columns
abun.mod.predict<-unique(abun.mod.predict[,c("Predicted","SE")])

#remove row.names
row.names(abun.mod.predict)<-NULL

#add c-hat goodness of fit
abun.mod.predict$GOF_c_hat<-mod.c.hat

colnames(abun.mod.predict)<-c("Mean abundance\n(birds per sampling location)","SE abundance","GOF (c-hat)")
  
panderOptions('table.alignment.default', function(df) ifelse(sapply(df, is.numeric), 'right', 'left'))
panderOptions("table.split.table", Inf)

pander(abun.mod.predict)

Mean abundance (birds per sampling location)	SE abundance	GOF (c-hat)
2.174	0.2008	3.286

Format and save vegetation data layer from 2019 Landfire dataset (30-m resolution)

#load vegetation raster layer
veg.raster<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LF2016_EVT_200_CONUS/Tif/LC16_EVT_200.tif")

#get extent of basemap
basemapExt<-attr(basemap, "bb")

#convert map extent from lat/long to UTM
basemapExt.df<-data.frame("long"=c(basemapExt$ll.lon, basemapExt$ur.lon),"lat"=c(basemapExt$ll.lat, basemapExt$ur.lat))

basemap_extent_sp<-SpatialPoints(cbind(basemapExt.df$long, basemapExt.df$lat), proj4string=CRS("+proj=longlat"))

#project to AEA
basemap_extent_AEA <- spTransform(basemap_extent_sp, projection(veg.raster))

#crop veg raster
veg.raster_crop<-crop(veg.raster, extent(basemap_extent_AEA))

#bring in States
#plot(states.sub)
#plot(veg.raster_crop)
#plot(sharp_patches_sub, col="blue")

#projection(sharp_patches_sub)

#reproject sharp_patches_sub
states_proj<-spTransform(states.sub, CRS(projection(veg.raster_crop)))

#dissolve states into one polygon, but first fix self-intersection regions with gBuffer, width=0
states_proj_valid<-gBuffer(states_proj, byid=TRUE, width=0)

#dissolve using sp::aggregate
states_proj_dissolve<-aggregate(states_proj_valid, dissolve = TRUE)

#plot(states_proj_dissolve)

#add buffer around polygons to give a little extra wiggle room when masking veg raster
states_proj_dissolve_buffer<-gBuffer(states_proj_dissolve, byid=TRUE, width=50)

#plot(states_proj_dissolve_buffer)

#now mask raster
veg.raster.mask<-mask(veg.raster_crop, states_proj_dissolve_buffer)

plot(veg.raster.mask)

#save cropped raster
writeRaster(veg.raster_crop, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_AEA_crop.tif",overwrite=TRUE)

#use gdalUtils::gdalwarp; GDAL is better for reprojecting rasters, and much faster than projectRaster(), when files are large
gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_AEA_crop.tif",
                    dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_latlong_crop.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",overwrite=TRUE,verbose=TRUE)

#save masked raster
writeRaster(veg.raster.mask, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_AEA_mask.tif",overwrite=TRUE)

#use gdalUtils::gdalwarp; GDAL is better for reprojecting rasters, and much faster than projectRaster(), when files are large
gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_AEA_mask.tif",
                    dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_latlong_mask.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",overwrite=TRUE,verbose=TRUE)

Figure 3. Landfire Ecosystem Class raster layer.

#read in cropped raster in lat/long
veg.raster.mask.latlong<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_latlong_mask.tif")

#ratify raster (function for defining raster values as categorical, i.e., factor levels)
r <- veg.raster.mask.latlong

r <- ratify(r)

rat <- levels(r)[[1]]

#read in attribute table
LFattributes<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LF2016_EVT_200_CONUS/CSV_Data/LF16_EVT_200.csv")

#subset LFattributes
LFattributes.sub<-subset(LFattributes, VALUE %in% rat$ID)

#create numeric column for EVT_CLASS
LFattributes.sub$EVT_CLASS_value<-as.integer(as.factor(LFattributes.sub$EVT_CLASS))
unique(LFattributes.sub$EVT_CLASS_value)

#get data.frame of Values and EVT_CLASS_value
EVT_CLASS.df<-unique(LFattributes.sub[,c("EVT_CLASS","VALUE","EVT_CLASS_value")])

#reclassify values
reclass.df <- data.frame(id=EVT_CLASS.df$VALUE, v=EVT_CLASS.df$EVT_CLASS_value)

#fix NA
#reclass.df$v<-ifelse(is.na(reclass.df$v),-9999,reclass.df$v)
reclass.r <- subs(r, reclass.df,subsWithNA=FALSE)

reclass.rat<- ratify(reclass.r)

rat <- levels(reclass.rat)[[1]]

#fix NA
LFattributes.sub$EVT_CLASS_value<-ifelse(is.na(LFattributes.sub$EVT_CLASS_value),-9999,LFattributes.sub$EVT_CLASS_value)

#make small df to merge
LFattributes.sub.small<-unique(LFattributes.sub[,c("EVT_CLASS","EVT_CLASS_value")])
names(LFattributes.sub.small)<-c("Ecosystem_type","ID")

rat_merge<-merge(rat, LFattributes.sub.small, by="ID",all.x=TRUE)
#rat$ID<-as.integer(unique(LFattributes.sub$EVT_CLASS_value))
#rat$Ecosystem_type<-unique(as.character(EVT_CLASS.df$EVT_CLASS))
rat_merge$band<-as.integer(rat_merge$ID)

#set factor levels
levels(reclass.r) <- rat_merge

#write raster
writeRaster(reclass.r, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask.tif",overwrite=TRUE)

#define veg class colors
myColors<-c("forestgreen", "lightgreen","seagreen", "goldenrod2", "purple","blue","orange","skyblue","salmon","cadetblue","gray40")

#save to not have to evaluate later (takes long time!)
png(filename="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_3_Land_Fire_2016_veg_layer.png",width=7, height=8, units="in", res=300)

  plot(reclass.r,col=myColors,legend=FALSE)
  legend("topleft", legend=unique(rat_merge$Ecosystem_type),fill=myColors)

dev.off()

#Load the saved image of raster veg layer
include_graphics("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_3_Land_Fire_2016_veg_layer.png")

Figure 4. Format and save Landfire elevation dataset (30-m resolution)

#load elevation raster layer
elev.raster<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LF2020_Elev_220_CONUS/Tif/LC20_Elev_220.tif")

#get extent of basemap
basemapExt<-attr(basemap, "bb")

#convert map extent from lat/long to UTM
basemapExt.df<-data.frame("long"=c(basemapExt$ll.lon, basemapExt$ur.lon),"lat"=c(basemapExt$ll.lat, basemapExt$ur.lat))

basemap_extent_sp<-SpatialPoints(cbind(basemapExt.df$long, basemapExt.df$lat), proj4string=CRS("+proj=longlat"))

#project to AEA
basemap_extent_AEA <- spTransform(basemap_extent_sp, projection(elev.raster))

#crop veg raster
elev.raster_crop<-crop(elev.raster, extent(basemap_extent_AEA))

#reproject sharp_patches_sub
states_proj<-spTransform(states.sub, CRS(projection(elev.raster_crop)))

#dissolve states into one polygon, but first fix self-intersection regions with gBuffer, width=0
states_proj_valid<-gBuffer(states_proj, byid=TRUE, width=0)

#dissolve using sp::aggregate
states_proj_dissolve<-aggregate(states_proj_valid, dissolve = TRUE)

#plot(states_proj_dissolve)

#add buffer around polygons to give a little extra wiggle room when masking veg raster
states_proj_dissolve_buffer<-gBuffer(states_proj_dissolve, byid=TRUE, width=50)

#plot(sharp_patches_sub_proj_dissolve_buffer)

#now mask raster
elev.raster.mask<-mask(elev.raster_crop, states_proj_dissolve_buffer)

#save cropped raster
writeRaster(elev.raster_crop, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_AEA_crop.tif",overwrite=TRUE)

#use gdalUtils::gdalwarp; GDAL is better for reprojecting rasters, and much faster than projectRaster(), when files are large
gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_AEA_crop.tif",
                    dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_crop.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",overwrite=TRUE,verbose=TRUE)

#save masked raster
writeRaster(elev.raster.mask, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_AEA_mask.tif",overwrite=TRUE)

#use gdalUtils::gdalwarp; GDAL is better for reprojecting rasters, and much faster than projectRaster(), when files are large
gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_AEA_mask.tif",
                    dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_mask.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",overwrite=TRUE,verbose=TRUE)

#plot elevation
elev.raster.mask.latlong<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_mask.tif")

#save to not have to evaluate later (takes long time!)
png(filename="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_4_Land_Fire_2016_elev_layer.png",width=7, height=8, units="in", res=300)

  plot(elev.raster.mask.latlong, col=parula(100))

dev.off()

#Load the saved image of raster veg layer
include_graphics("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_4_Land_Fire_2016_elev_layer.png")

Extract proportions of each veg class within 5000-m radius buffer around each sampling point.

#add PointID column
count.data$PointID<-paste("Point",seq(1,nrow(count.data)),sep="_")

#get all unique sampling locations
sampling.locs<-unique(count.data[,c("PointID","Longitude","Latitude")])
sampling.locs$ID<-seq(1,length(row.names(sampling.locs)),by=1)

#add buffers of 750-m radii to sampling locs
# circular buffers with a radius of 1000 m
buffers <- dismo::circles(sampling.locs[c("Longitude","Latitude")], d=5000, lonlat=TRUE, dissolve=FALSE)

#coerce to polygons
buffers.polygon <- polygons(buffers)

reclass.r<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask.tif")

#extract values (IDs are circular buffer polygons)
extracted.veg.cell.counts<-raster::extract(reclass.r, buffers.polygon, factors=TRUE, df=TRUE)
names(extracted.veg.cell.counts)[names(extracted.veg.cell.counts)=="category"]<-"Ecosystem_type"

#make ID and Vegetation_class factors
extracted.veg.cell.counts$ID<-as.factor(extracted.veg.cell.counts$ID)
extracted.veg.cell.counts$Ecosystem_type<-as.factor(extracted.veg.cell.counts$Ecosystem_type)

#add CellCount column
extracted.veg.cell.counts$CellCount<-1

#get total cell count per ID
total.cell.count<-aggregate(CellCount~ID, data=extracted.veg.cell.counts, FUN=sum)
names(total.cell.count)[names(total.cell.count)=="CellCount"]<-"TotalCellCount"

#aggregate to get proportions
extracted.veg.aggregate<-stats::aggregate(CellCount~ID+Ecosystem_type, data=extracted.veg.cell.counts, FUN=sum)

#merge
extracted.merge<-merge(extracted.veg.aggregate, total.cell.count, by=c("ID"),all.x=TRUE)

#Compute proportions (CellCount/TotalCellCount)
extracted.merge$Ecosystem_type_proportion<-extracted.merge$CellCount/extracted.merge$TotalCellCount

#Now merge back with sampling locs
veg.proportions.df<-merge(sampling.locs, extracted.merge, by="ID",all.x=TRUE)

#melt and cast
veg.data.melt<-melt(veg.proportions.df, id.vars=c("PointID","ID","Longitude","Latitude", "Ecosystem_type"), measure.vars=c("Ecosystem_type_proportion"))

#cast and fill any NAs with 0s
pred.veg.proportions.df.cast<-dcast(PointID+ID+Longitude+Latitude~Ecosystem_type,data=veg.data.melt, fill=0)

#remove NA column
pred.veg.proportions.df.cast$`NA`<-NULL

#merge with count data
count.data.merge<-merge(count.data, veg.proportions.df[,c("PointID","Ecosystem_type","Ecosystem_type_proportion")], by=c("PointID"),all.x=TRUE)

#save .csv
write.csv(count.data.merge, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Detection_Data_with_EcoSys_covs.csv",row.names=FALSE)

Extract elevation (km) at each sampling point.

#load elevation data
elev.raster.mask.latlong<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_mask.tif")

# get elevation at sampling locations
spatialPts<- SpatialPoints(cbind(count.data$Longitude, count.data$Latitude), proj4string=CRS(projection(elev.raster.mask.latlong)))

#extract values (IDs are circular buffer polygons)
extracted.elev.cell.values<-raster::extract(elev.raster.mask.latlong, spatialPts, fun=mean, factors=FALSE, df=TRUE)

#make ID and Vegetation_class factors
extracted.elev.cell.values$ID<-as.factor(extracted.elev.cell.values$ID)

#add PointID
extracted.elev.cell.values$PointID<-paste("Point",seq(1,nrow(count.data)),sep="_")

#merge with count data
elev.data.merge<-merge(count.data, extracted.elev.cell.values[,c("PointID","LandFire_Elevation_latlong_mask")], by=c("PointID"),all.x=TRUE)

#save .csv
write.csv(extracted.elev.cell.values, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Sampling_Locs_elevation.csv",row.names=FALSE)

Add covariates to count data (e.g., elevation).

#read in new data with Veg layer covariates
detectionData_with_veg_covs<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Detection_Data_with_EcoSys_covs.csv")

#melt and cast
landcover.melt<-melt(detectionData_with_veg_covs, id.vars=c("PointID","Hexagon","Jun", "Jul","Aug","Longitude","Latitude","Visit_1","Visit_2","Visit_3", "Ecosystem_type"), measure.vars=c("Ecosystem_type_proportion"))

#cast and fill any NAs with 0s
lancover.proportions.cast<-dcast(PointID+Hexagon+Longitude+Latitude+Jun+Jul+Aug+Visit_1+Visit_2+Visit_3~Ecosystem_type,data=landcover.melt, fill=0)

#load elevation data
elev.df<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Sampling_Locs_elevation.csv")

#add other covs
covs.merge.1<-merge(lancover.proportions.cast, elev.df, by="PointID",all.x=TRUE)

#fix names
names(covs.merge.1)[names(covs.merge.1)=="LandFire_Elevation_latlong_mask"]<-"Elevation_m"

#remove Covariate == NA
covs.merge.1$`NA`<-NULL

#save covs.merge.1
write.csv(covs.merge.1, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Detection_data_ALL_Covs.csv",row.names=FALSE)

Build prediction surface from landcover and elevation rasters

#build prediction surface (and save .csv)

#create a new blank raster with all values = 1 with dimensions for prediction
veg.raster.mask.latlong.class<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask.tif")

#create a raster that is not quite so hi-res as original raster (to use for extraction)
#compute 100-m resolution
res100<-res(veg.raster.mask.latlong.class)[1]*100/30

#downsample raster with gdalwarp (use r="mode" for categorical data)
gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask.tif",
                    dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask_100m.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",tr=c(res100,res100),r="mode",overwrite=TRUE,verbose=TRUE)

veg.raster.mask.latlong_100<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask_100m.tif")

#now compute 5000-m resolution to get polygons for extraction
new.res<-res(veg.raster.mask.latlong.class)[1]*5000/30

#downsample raster with gdalwarp (use r="mode" for categorical data)
gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask.tif",dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask_5000m.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",tr=c(new.res,new.res),r="mode",overwrite=TRUE,verbose=TRUE)

veg.raster.mask.latlong_5000<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_EcoSys_EVT_Class_latlong_mask_5000m.tif")

plot(veg.raster.mask.latlong_5000, col=parula(13))
# 
# plot(veg.raster.mask.latlong_300)
# dim(veg.raster.mask.latlong_300)
# res(veg.raster.mask.latlong_300)

#ratify raster (function for defining raster values as categorical, i.e., factor levels)
r_5000 <- veg.raster.mask.latlong_5000
r_5000_rat <- ratify(r_5000)

rat_5000 <- levels(r_5000_rat)[[1]]

#read in attribute table
rat_merge_5000<-merge(rat_5000, LFattributes.sub.small, by="ID",all.x=TRUE)

rat_merge_5000$band<-as.integer(rat_merge_5000$ID)

#set factor levels
levels(veg.raster.mask.latlong_5000) <- rat_merge_5000

new.raster.df<-as.data.frame(rasterToPoints(veg.raster.mask.latlong_5000))
colnames(new.raster.df)<-c("Longitude","Latitude","Value")

###############################################################################
#get all unique sampling locations
pred.sampling.locs<-unique(new.raster.df[,c("Longitude","Latitude")])
pred.sampling.locs$ID<-seq(1,length(row.names(pred.sampling.locs)),by=1)
pred.sampling.locs$ID<-as.factor(as.character(pred.sampling.locs$ID))

# circular buffers with a radius of 100m
pred.buffers <- dismo::circles(pred.sampling.locs[,c("Longitude","Latitude")], d=5000, lonlat=TRUE, dissolve=FALSE)

#coerce to polygons
pred.buffers.polygon <- polygons(pred.buffers)

#save as Rdata
# saveRDS(pred.buffers.polygon, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Predicted_Polygons.rds")
# 
# pred.buffers.polygon_test<-readRDS("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Predicted_Polygons.rds")
##############################################################################
#Vegetation Class proportions
#create empty list
save.proportions<-list()

#loop through and extract proportions
for(i in 1:length(pred.buffers.polygon)){
  print(i)
  new.poly<-pred.buffers.polygon[i]
  pred.veg.cell.counts.new<-raster::extract(veg.raster.mask.latlong.class, new.poly, factors=TRUE, df=TRUE,na.rm=TRUE)

  #make ID (cell location ID) and Vegetation_class factors
  pred.veg.cell.counts.new$ID<-as.character(i)
  names(pred.veg.cell.counts.new)[names(pred.veg.cell.counts.new)=="category"]<-"Ecosystem_type"

  #add CellCount column
  pred.veg.cell.counts.new$CellCount<-1

  #get total cell count per Vegetation_class
  total.cell.count.pred.new<-stats::aggregate(CellCount~ID, data=na.omit(pred.veg.cell.counts.new), FUN=sum)
  names(total.cell.count.pred.new)[names(total.cell.count.pred.new)=="CellCount"]<-"TotalCellCount"

  #aggregate to get proportions
  pred.veg.aggregate.new<-stats::aggregate(CellCount~ID+Ecosystem_type, data=pred.veg.cell.counts.new, FUN=sum)

#merge
pred.veg.merge<-merge(pred.veg.aggregate.new, total.cell.count.pred.new, by=c("ID"),all.x=TRUE)

#Compute proportions (CellCount/TotalCellCount)
pred.veg.merge$Ecosystem_type_proportion<-pred.veg.merge$CellCount/pred.veg.merge$TotalCellCount

save.proportions<-rbind(save.proportions, pred.veg.merge)

}

#save as .csv
write.csv(save.proportions, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Predicted_Surface_proportions.csv", row.names=FALSE)

pred.surface.proportions<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Predicted_Surface_proportions.csv")

#extract values (IDs are circular buffer polygons)
#pred.veg.cell.counts<-raster::extract(veg.raster.mask.latlong.class, pred.buffers.polygon, factors=TRUE, df=TRUE)

#Now merge back with sampling locs
pred.veg.proportions.df<-merge(pred.surface.proportions, pred.sampling.locs, by="ID",all.x=TRUE)

#melt and cast
pred.veg.proportions.df.melt<-melt(pred.veg.proportions.df, id.vars=c("ID","Longitude","Latitude", "Ecosystem_type"), measure.vars=c("Ecosystem_type_proportion"))

#cast and fill any NAs with 0s
pred.veg.proportions.df.cast<-dcast(ID+Longitude+Latitude~Ecosystem_type,data=pred.veg.proportions.df.melt, fill=0)

#remove NA column
pred.veg.proportions.df.cast$`NA`<-NULL
#dim(pred.veg.proportions.df.cast)
#unique(pred.veg.proportions.df.cast$ID)

#save predictive surface as .csv
write.csv(pred.veg.proportions.df.cast, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Predicted_Surface__Landfire_proportions_5000m.tif", row.names=FALSE)

#elevation raster
elev.raster.mask.latlong<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_mask.tif")

#downsample from 30 to 1000m resolution
gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_mask.tif",
                    dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_mask_5000.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",tr=c(new.res,new.res),r="average",overwrite=TRUE,verbose=TRUE)

elev.raster.mask.latlong_5000<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/LandFire_Elevation_latlong_mask_5000.tif")

elev.raster.df<-as.data.frame(rasterToPoints(elev.raster.mask.latlong_5000))
colnames(elev.raster.df)<-c("Longitude","Latitude","Elevation_m")

#merge covariates 
all.covs.merge<-Reduce(function(x, y) merge(x, y, by=c("Latitude","Longitude"), all.x=TRUE), list(pred.veg.proportions.df.cast, elev.raster.df))

write.csv(all.covs.merge, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/ALL_covs_pred_surface.csv", row.names=FALSE)

Figure 5. Explore covariates for colinearity (Pearson’s correlation) using the unmarkedFrame

# #read in new data with Veg layer covariates
all.data<-na.omit(read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Detection_data_ALL_Covs.csv"))


#reduce to include only variables we model (USE LFattributes for 12 landcover types)
#Open Tree Canopy = Evergreen Open Tree Canopy
covariate.df<-all.data[,c("Elevation_m",
                          "Closed.tree.canopy",
                           #"Dwarf.shrubland",
                           "Herbaceous.grassland",
                           "Herbaceous.shrub.steppe",
                           "No.Dominant.Lifeform",
                           "Non.vegetated",
                           "Open.tree.canopy",
                           "Open.Tree.Canopy",
                           "Shrubland",
                           "Sparse.tree.canopy",
                           "Sparsely.vegetated"
                             )]
#scale covs
covariate.df<-sapply(covariate.df, function(x) scale(x) )  

#create a matrix of these data
cov_matrix <-as.matrix(covariate.df)

#get correlation matrix
cov_cor_matrix <-Hmisc::rcorr(cov_matrix)

cor_mat.df<-as.data.frame(cov_cor_matrix[1])

#plot the correlation matrix with GGally package
corrPlot<-GGally::ggcorr(cor_mat.df, midpoint=0,label=TRUE)

#save plot
jpeg(file = "/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_5_Correlation_matrix_covs.jpeg", width=12, height=12, units="in",res=600)

print(corrPlot)

dev.off()

#show plot by loading saved image
include_graphics("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_5_Correlation_matrix_covs.jpeg")

Figure 6. Explore covariate relationships to detection probabilty

#read in detection data with covariates
all.data<-na.omit(read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Detection_data_ALL_Covs.csv"))

#make list of detection process covariates (non-factors only)
detformulaList<-c(~Elevation_m, ~Closed.tree.canopy, ~Herbaceous.grassland,~Herbaceous.shrub.steppe, ~No.Dominant.Lifeform, ~Non.vegetated, ~Open.tree.canopy, ~Open.Tree.Canopy,~Shrubland,~Sparse.tree.canopy,
                  ~Sparsely.vegetated)


#build new unmarked frame
new.umf<-makePOccu_umf_with_covs(DataIn=all.data)

#get obsCovs as data.frame
obsCovs.df<-as.data.frame(new.umf@obsCovs)

#empty list to fill with results
occu.mod.save<-list()

#create empty list to save predicted detection Probs
DetCovariatePred<-list()

#for loop to iterate over detection covariates and produce predicted estimates for detection prob
for(j in 1:length(detformulaList)){
    #print(detformulaList[[j]])
    
    #fit model
    new.occ_model <- occu(detformulaList[[j]]~ 1,  data=new.umf)
    
    #get det covariate name
    varName<-gsub("~","",as.character(detformulaList[[j]]))[2]
  
    #get range of values from data to simulate new data
    varData<-obsCovs.df[,names(obsCovs.df)==varName]
  
    #simulate new data
    sim.data<-data.frame(seq(min(varData,na.rm=TRUE), max(varData,na.rm=TRUE),length.out=100))
    colnames(sim.data)<-c(varName)
    #get predicted values using sim.data
    pred.df<-predict(new.occ_model, type="det", newdata=sim.data, appendData=TRUE)
    colnames(pred.df)[5]<-"det_cov"
    
    #get original covariate data
    cov.data<-all.data[,c(print(varName))]
    meanCov <- mean(cov.data)
    sdCov <- sd(cov.data)
    
    #backtransform scaled variable
    pred.df$det_cov_bt<- t((t(pred.df$det_cov) * sdCov) + meanCov)
    
    #add column with covariate name
    pred.df$CovariateName<-varName
    
    DetCovariatePred<-rbind(DetCovariatePred, pred.df)
    
  #compile covariates for each species
occu.mod.save<-rbind(occu.mod.save, DetCovariatePred)
}

## [1] "Elevation_m"
## [1] "Closed.tree.canopy"
## [1] "Herbaceous.grassland"
## [1] "Herbaceous.shrub.steppe"
## [1] "No.Dominant.Lifeform"
## [1] "Non.vegetated"
## [1] "Open.tree.canopy"
## [1] "Open.Tree.Canopy"
## [1] "Shrubland"
## [1] "Sparse.tree.canopy"
## [1] "Sparsely.vegetated"

#make CovariateName a factor
occu.mod.save$CovariateName<-factor(occu.mod.save$CovariateName, levels=c("Elevation_m","Closed.tree.canopy",  "Herbaceous.grassland","Herbaceous.shrub.steppe","Open.tree.canopy","Open.Tree.Canopy","Shrubland","Sparse.tree.canopy","Sparsely.vegetated","Non.vegetated","No.Dominant.Lifeform"))

#now plot 
  #plot p(Occupancy) vs. Vegetation class proportions
  det.prob.plot<-ggplot(data=occu.mod.save, aes(x=det_cov_bt, y=Predicted))+
    geom_point(aes(color=CovariateName),alpha=0.9,size=1)+
    geom_smooth(method = "glm", formula = y~x, method.args = list(family = binomial), aes(color=CovariateName, fill=CovariateName),alpha=0.2,fullrange=FALSE)+
    # scale_color_manual(values=myColors)+
    # scale_fill_manual(values=myColors)+
    #ylim(c(0,1))+
    #coord_cartesian(expand = FALSE, xlim=c(0,1),ylim=c(0.99,1))+
theme(
      panel.background = element_rect(fill="white",color="black"),
      panel.border = element_rect(fill="transparent",color="black"),
      axis.text.x = element_text(size=8),
      axis.text.y = element_text(size=8),
      legend.position="bottom"
    )+
    labs(x="Observation Covariate", y="Detection Probability")
  
  #facet over Vegetation_class
  det.prob.plot.facet<-det.prob.plot+facet_wrap(CovariateName~., scales="free",ncol=4)
  print(det.prob.plot.facet)

  #create filename with timestamp
  new.filename<-paste("Fig_6_Detection_Probability_covariates",".png",sep="")
  
  ggsave(det.prob.plot.facet, file=paste("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/",new.filename, sep=""),width=9, height=9, dpi=600)

Fit occupancy models with site-level covariates.

#build new unmarked frame
new.umf<-makePOccu_umf_with_covs(DataIn=all.data)

  #fit model with Veg covariates
  mod.global<-NULL
  mod.global<-occu(~Visit ~Elevation_m+
                           Closed.tree.canopy+
                           #Dwarf.shrubland+
                           Herbaceous.grassland+
                           #Herbaceous.shrub.steppe+
                           Open.tree.canopy+
                           Open.Tree.Canopy+
                           #Shrubland+
                           #Sparse.tree.canopy+
                           Sparsely.vegetated+
                           #Non.vegetated+
                           No.Dominant.Lifeform
                           , data=new.umf,se=TRUE)
  
# #model selection step
# dredge.mod <- dredge(mod.global)
# dredge.mod
# dredge.mod.top<-subset(dredge.mod, delta < 4)

# Visualize the model selection table:
# par(mar = c(3,5,6,4))
# plot(dredge.mod.top, labAsExpr = TRUE)


# Model average models with delta AICc < 4
# mod.average<-model.avg(dredge.mod.top)
# summary(mod.average)

# Model-averaged coefficients:  
# (full average) 
#                             Estimate Std. Error z value Pr(>|z|)    
# psi(Int)                  -5.455e-01  3.168e-01   1.722  0.08503 .  
# psi(Open.tree.canopy)      2.971e+00  1.243e+00   2.390  0.01687 *  
# p(Int)                     8.212e-01  2.200e-01   3.732  0.00019 ***
# p(VisitVisit_2)           -3.394e-01  2.623e-01   1.294  0.19565    
# p(VisitVisit_3)           -1.754e+00  2.788e-01   6.291  < 2e-16 ***
# psi(Open.Tree.Canopy)      2.348e+00  4.627e+00   0.507  0.61182    
# psi(No.Dominant.Lifeform) -1.958e-01  5.778e-01   0.339  0.73475    
# psi(Sparsely.vegetated)    1.143e+01  3.553e+01   0.322  0.74773    
# psi(Closed.tree.canopy)    3.487e-02  2.479e-01   0.141  0.88814    
# psi(Herbaceous.grassland) -1.499e-02  2.625e-01   0.057  0.95445    
# psi(Non.vegetated)        -1.087e-02  4.705e-01   0.023  0.98157    
# psi(Elevation_m)          -6.114e-06  2.349e-04   0.026  0.97923    
#  
# (conditional average) 
#                             Estimate Std. Error z value Pr(>|z|)    
# psi(Int)                  -5.455e-01  3.168e-01   1.722  0.08503 .  
# psi(Open.tree.canopy)      2.971e+00  1.243e+00   2.390  0.01687 *  
# p(Int)                     8.212e-01  2.200e-01   3.732  0.00019 ***
# p(VisitVisit_2)           -3.394e-01  2.623e-01   1.294  0.19565    
# p(VisitVisit_3)           -1.754e+00  2.788e-01   6.291  < 2e-16 ***
# psi(Open.Tree.Canopy)      6.622e+00  5.663e+00   1.169  0.24226    
# psi(No.Dominant.Lifeform) -7.887e-01  9.367e-01   0.842  0.39979    
# psi(Sparsely.vegetated)    4.352e+01  5.841e+01   0.745  0.45617    
# psi(Closed.tree.canopy)    2.125e-01  5.803e-01   0.366  0.71425    
# psi(Herbaceous.grassland) -1.251e-01  7.490e-01   0.167  0.86738    
# psi(Non.vegetated)        -9.173e-02  1.364e+00   0.067  0.94639    
# psi(Elevation_m)          -4.508e-05  6.364e-04   0.071  0.94352    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#refit using significant covariates
  #fit model with Veg covariates
  mod.global<-NULL
  mod.global<-occu(~Visit ~Elevation_m+
                           Closed.tree.canopy+
                           #Dwarf.shrubland+
                           #Herbaceous.grassland+
                           #Herbaceous.shrub.steppe+
                           Open.tree.canopy+
                           Open.Tree.Canopy
                           #Shrubland+
                           #Sparse.tree.canopy+
                           #Sparsely.vegetated+
                           #Non.vegetated
                           #No.Dominant.Lifeform
                           , data=new.umf,se=TRUE)


  mod.gof<-NULL
  mod.gof<- try(mb.gof.test(mod.global, nsim = 10, plot.hist = FALSE, report = NULL, parallel = FALSE))
  
  #get c-hat
  if(!is.null(mod.gof)){
    mod.c.hat<-mod.gof$c.hat.est
  } else {
    mod.c.hat<-NA
  }
                   
  #get predicted estimates
  mod.predict<-predict(mod.global, type="state",appendData=TRUE)
  #add original unscaled covariates
  #mod.predict<-cbind(mod.predict, all.data)
  
  #add GOF_cHat
  mod.predict$GOF_c_hat<-mod.c.hat
  
  #simplify
  #mod.predict.sub<-unique(mod.predict[,c("PointID","Longitude","Latitude","Predicted","SE","lower","upper","GOF_c_hat","Elevation_m","Closed.tree.canopy","Herbaceous.grassland","Herbaceous.shrub.steppe","Open.tree.canopy","Open.Tree.Canopy","Shrubland","Sparse.tree.canopy","Sparsely.vegetated","Non.vegetated","No.Dominant.Lifeform")])
  
 mod.predict.sub<-unique(mod.predict[,c("PointID","Longitude","Latitude","Predicted","SE","lower","upper","GOF_c_hat","Elevation_m","Closed.tree.canopy","Open.tree.canopy","Open.Tree.Canopy")])

Figure 7. Plot of occupancy model-estimated beta coefficients for site-level covariates. Points and bars of mean and upper and lower 95% CIs are shown in blue for non-zero crossing and gray for zero crossing covariates.

#########################################################################
#Plot the beta coefficients 

  #get model summary
  mod.summary<-summary(mod.global)

## 
## Call:
## occu(formula = ~Visit ~ Elevation_m + Closed.tree.canopy + Open.tree.canopy + 
##     Open.Tree.Canopy, data = new.umf, se = TRUE)
## 
## Occupancy (logit-scale):
##                     Estimate       SE      z P(>|z|)
## (Intercept)        -1.045626 0.497773 -2.101 0.03568
## Elevation_m        -0.000223 0.000711 -0.313 0.75403
## Closed.tree.canopy  0.490350 0.677592  0.724 0.46927
## Open.tree.canopy    4.696161 1.643803  2.857 0.00428
## Open.Tree.Canopy    7.815849 6.676846  1.171 0.24176
## 
## Detection (logit-scale):
##              Estimate    SE     z  P(>|z|)
## (Intercept)     0.821 0.220  3.73 1.90e-04
## VisitVisit_2   -0.339 0.262 -1.29 1.96e-01
## VisitVisit_3   -1.754 0.279 -6.29 3.17e-10
## 
## AIC: 786.7018 
## Number of sites: 266
## optim convergence code: 0
## optim iterations: 58 
## Bootstrap iterations: 0

  new.coefs<-as.data.frame(mod.summary$state)
  new.coefs$Variable<-row.names(new.coefs)
  row.names(new.coefs)<-NULL
  
  #coompute 95% CIs
  new.coefs$`Lower_95%_CI` <- new.coefs$Estimate - 1.96*new.coefs$SE 
  new.coefs$`Upper_95%_CI` <- new.coefs$Estimate + 1.96*new.coefs$SE

#first remove intercepts
coefs.occu.save<-subset(new.coefs, Variable != "(Intercept)")

#add ZeroCross column
coefs.occu.save$ZeroCross<-ifelse(coefs.occu.save$`Lower_95%_CI`>0 | coefs.occu.save$`Upper_95%_CI` < 0,"No_Zero_Cross","Zero_Cross")

#make a factor
coefs.occu.save$ZeroCross<-as.factor(coefs.occu.save$ZeroCross)

#sort Variables
coefs.occu.save.order<-coefs.occu.save[order(coefs.occu.save$Estimate, decreasing = FALSE),]

#set factor levels for Variable
coefs.occu.save.order$Variable<-factor(coefs.occu.save.order$Variable, levels=unique(coefs.occu.save.order$Variable))

#set colors
myBetaColors<-c("royalblue2","gray60")

  #plot beta coefficients
  betaCoefPlot_occu<-ggplot(data=coefs.occu.save.order, aes(x=Variable, y=Estimate))+
        geom_hline(yintercept=0,linetype=2)+
        geom_errorbar(aes(ymin=`Lower_95%_CI`, ymax=`Upper_95%_CI`, color = ZeroCross),width=0)+
        geom_point(aes(color=ZeroCross),alpha=0.8)+
        scale_y_continuous(limits=c(-5000,5000))+
        #coord_cartesian(ylim=c(-180,180))+
        scale_color_manual(values=myBetaColors)+
        coord_flip(ylim=c(-10,10))+
#scale_color_gradient2(low="red",mid="gray40",high="royalblue2",midpoint = 0)+
theme(legend.position="none",
              panel.background = element_rect(fill="white",color="black"),
              panel.border = element_rect(fill="transparent",color="black"))+
    labs(y="Beta Coefficient Estimate", x="Site-level Covariate")
  
betaCoefPlot_occu

ggsave(betaCoefPlot_occu, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_7_Beta_coefficient_plot.png",width=6, height=6,dpi=600)

Figure 8. Probability of occupancy in relation to site-level covariates.

#melt
mod.predict.melt<-melt(mod.predict.sub, id.vars=c("PointID","Predicted","SE","lower","upper"), measure.vars=c("Elevation_m","Closed.tree.canopy","Open.tree.canopy","Open.Tree.Canopy"), variable.name = "Variable", value.name="value")

#define veg class colors
#myColors<-c("#538233","#a9d08e", "#f9bf8e", "#feff00", "#01b0f0","#0070c0","#000000","#a6a6a6","magenta2","royalblue3","seagreen4","orangered3")

#levels(mod.predict.melt$Variable)

  #plot p(Occupancy) vs. Vegetation class proportions
  mod.occu.plot<-ggplot(data=mod.predict.melt, aes(x=value, y=Predicted,group=Variable))+
    geom_point(aes(color=Variable),alpha=0.3,size=1)+
     stat_smooth(method = "glm", formula = y~x, method.args = list(family = binomial), aes(color=Variable, fill=Variable),alpha=0.2,fullrange=FALSE,se=TRUE)+
    # scale_color_manual(values=myColors)+
    # scale_fill_manual(values=myColors)+
    #ylim(c(0,1))+
theme(
      panel.background = element_rect(fill="white",color="black"),
      panel.border = element_rect(fill="transparent",color="black"),
      axis.text.x = element_text(size=8),
      axis.text.y = element_text(size=8),
      legend.position="bottom"
    )+
    labs(x="Variable", y="Predicted p(Occupancy)")
  
  #facet over Vegetation_class
  mod.occu.plot.facet<-mod.occu.plot+facet_wrap(Variable~.,scales="free",ncol=2)
   
  mod.occu.plot.facet

  #create filename with timestamp
  new.filename<-paste("Fig_8_Predicted_Occupancy_models_with_veg_class",".png",sep="")
  
  ggsave(mod.occu.plot.facet, file=paste("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/",new.filename, sep=""),width=12, height=12, dpi=600)

Use occupancy model with covs to get predicted estimates

#load predictive surface
pred_surface<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/ALL_covs_pred_surface.csv")

#remove unwanted columns
# pred_surface$Dwarf.shrubland<-NULL
# pred_surface$Herbaceous.shrub.steppe<-NULL
# pred_surface$No.Dominant.Lifeform<-NULL
# pred_surface$Non.vegetated<-NULL
# pred_surface$Shrubland<-NULL
# pred_surface$Sparse.tree.canopy<-NULL
# pred_surface$Sparsely.vegetated<-NULL


#get predicted estimates of avian occupancy as a function of vegetation cover on map
occ_pred <- predict(mod.global,
                    newdata = as.data.frame(pred_surface),
                    type = "state", appendData=TRUE)

suppressMessages({

# add to prediction surface
pred_occ <- bind_cols(pred_surface, 
                      occ_prob = occ_pred$Predicted, 
                      occ_se = occ_pred$SE) %>% 
  select(Latitude, Longitude, occ_prob, occ_se)

  })

 #melt to plot in ggplot
pred_occ.melt<-melt(pred_occ, id.vars=c("Latitude","Longitude"),
                    variable.name = "Measure",value.name="p(Occupancy)")

Figure 9. Predicted probability of whip-poor-will occupancy (mean ± SE) in USFWS Northeast region.

  map.occu.plot<-ggmap(basemap,darken = c(0.4, "black"))+
    geom_tile(data=pred_occ.melt, aes(x=Longitude, y=Latitude, fill=`p(Occupancy)`,color=`p(Occupancy)`))+
  geom_map(data=states.poly, map=states.poly, aes(y=lat, x=long, map_id=id),color=alpha("white",0.3),fill="transparent", inherit.aes = FALSE,size=0.1)+
    scale_fill_gradientn(colors=ocean.thermal(100),na.value = "transparent",limits=c(0,1))+
    scale_color_gradientn(colors=ocean.thermal(100),na.value = "transparent",limits=c(0,1))+
  labs(x="Longitude",y="Latitude")+
    theme(axis.text=element_text(size=12),
      axis.title=element_text(size=14),
      axis.text.x=element_text(size=8),
      axis.text.y=element_text(size=8),
      panel.border=element_rect(fill="transparent", color="black"),
      legend.position = "right")+
    ggtitle(paste("WPWI"," - Predicted probability of occupancy",sep=""))+
    #geom_text_repel(check_overlap = FALSE, size=3, nudge.y=-10, nudge.x=-10,aes(x=Longitude, y=Latitude,color=PatchID, label=PatchID),data=PatchID_labels)+
    #guides(fill = guide_legend(title = "National Wildlife Refuge", label = FALSE))+
    ggsn::scalebar(x.min=mapExtent$ll.lon, x.max=mapExtent$ur.lon, y.min=mapExtent$ll.lat, y.max=mapExtent$ur.lat,transform=TRUE, dist_unit="km",dist=100, st.dist=0.03, height = 0.015, anchor=c(x=-70.5,y=38.55), border.size=0.5,st.size=3)

  map.occu.plot.facet.save<-map.occu.plot+facet_wrap(Measure~.,ncol=1)
  
  #create filename with timestamp
  new.filename<-paste(paste("Fig_9_WPWI_","Map_Predicted_Occupancy",sep=""),".png",sep="")
    
  ggsave(map.occu.plot.facet.save, file=paste("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/",new.filename, sep=""),width=14, height=8, dpi=600)
  
print(map.occu.plot.facet.save)

Fit abundance models and explore covariate relationships.

#read in detection data with covariates
all.data<-na.omit(read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/Detection_data_ALL_Covs.csv"))

#build new unmarked frame
new.umf<-makePCount_umf_with_covs(DataIn=all.data)

  #fit model with Veg covariates
  mod.global<-NULL
  mod.global<-pcount(~Visit ~Elevation_m+
                           Closed.tree.canopy+
                           #Dwarf.shrubland+
                           Herbaceous.grassland+
                           #Herbaceous.shrub.steppe+
                           Open.tree.canopy+
                           Open.Tree.Canopy+
                           #Shrubland+
                           #Sparse.tree.canopy+
                           Sparsely.vegetated+
                           #Non.vegetated+
                           No.Dominant.Lifeform
                           , data=new.umf,se=TRUE)
  
# #model selection step
# dredge.mod <- dredge(mod.global)
# dredge.mod
# dredge.mod.top<-subset(dredge.mod, delta < 4)

# Visualize the model selection table:
# par(mar = c(3,5,6,4))
# plot(dredge.mod.top, labAsExpr = TRUE)

# Model average models with delta AICc < 4
# mod.average<-model.avg(dredge.mod.top)
# summary(mod.average)

# Model-averaged coefficients:  
# (full average) 
#                             Estimate Std. Error z value Pr(>|z|)    
# lam(Int)                   2.096e-02  3.008e-01   0.070  0.94444    
# lam(Closed.tree.canopy)    3.274e-01  3.491e-01   0.938  0.34829    
# lam(Open.tree.canopy)      3.552e+00  5.409e-01   6.567  < 2e-16 ***
# lam(Open.Tree.Canopy)      5.794e+00  2.589e+00   2.238  0.02522 *  
# lam(Sparsely.vegetated)    3.367e+01  1.542e+01   2.183  0.02902 *  
# p(Int)                    -4.936e-01  1.751e-01   2.818  0.00483 ** 
# p(VisitVisit_2)           -3.846e-01  1.205e-01   3.191  0.00142 ** 
# p(VisitVisit_3)           -1.705e+00  1.621e-01  10.517  < 2e-16 ***
# lam(Herbaceous.grassland) -1.406e-01  3.312e-01   0.425  0.67111    
# lam(No.Dominant.Lifeform) -5.116e-02  2.792e-01   0.183  0.85460    
# lam(Elevation_m)           1.305e-05  7.074e-05   0.184  0.85368    
#  
# (conditional average) 
#                             Estimate Std. Error z value Pr(>|z|)    
# lam(Int)                   0.0209584  0.3007513   0.070  0.94444    
# lam(Closed.tree.canopy)    0.5377469  0.2950078   1.823  0.06833 .  
# lam(Open.tree.canopy)      3.5520350  0.5408537   6.567  < 2e-16 ***
# lam(Open.Tree.Canopy)      6.0753976  2.3065417   2.634  0.00844 ** 
# lam(Sparsely.vegetated)   35.1468185 14.0129940   2.508  0.01214 *  
# p(Int)                    -0.4935698  0.1751444   2.818  0.00483 ** 
# p(VisitVisit_2)           -0.3845809  0.1205179   3.191  0.00142 ** 
# p(VisitVisit_3)           -1.7051791  0.1621320  10.517  < 2e-16 ***
# lam(Herbaceous.grassland) -0.4167398  0.4582780   0.909  0.36316    
# lam(No.Dominant.Lifeform) -0.2162472  0.5419927   0.399  0.68990    
# lam(Elevation_m)           0.0002180  0.0001973   1.105  0.26918    

#refit using significant covariates
  #fit model with Veg covariates
  mod.global<-NULL
  mod.global<-pcount(~Visit ~
                           #Elevation_m+
                           #Closed.tree.canopy+
                           #Dwarf.shrubland+
                           Herbaceous.grassland+
                           #Herbaceous.shrub.steppe+
                           Open.tree.canopy
                           #Open.Tree.Canopy+
                           #Shrubland+
                           #Sparse.tree.canopy+
                           #Sparsely.vegetated+
                           #Non.vegetated+
                           #No.Dominant.Lifeform
                           , data=new.umf,se=TRUE)


  mod.gof<-NULL
  mod.gof<- Nmix.gof.test(mod.global, nsim = 20, plot.hist = FALSE,report = NULL, parallel = FALSE)
  
  #get c-hat
  if(!is.null(mod.gof)){
    mod.c.hat<-mod.gof$c.hat.est
  } else {
    mod.c.hat<-NA
  }
                   
  #get predicted estimates
  mod.predict<-predict(mod.global, type="state",appendData=TRUE)
  #add original unscaled covariates
  #mod.predict<-cbind(mod.predict, all.data)
  
  #add GOF_cHat
  mod.predict$GOF_c_hat<-mod.c.hat
  
  #simplify
  mod.predict.sub<-unique(mod.predict[,c("PointID","Longitude","Latitude","Predicted","SE","lower","upper","GOF_c_hat","Elevation_m","Closed.tree.canopy","Herbaceous.grassland","Herbaceous.shrub.steppe","Open.tree.canopy","Open.Tree.Canopy","Shrubland","Sparse.tree.canopy","Sparsely.vegetated","Non.vegetated","No.Dominant.Lifeform")])
  
 # mod.predict.sub<-unique(mod.predict[,c("PointID","Longitude","Latitude","Predicted","SE","lower","upper","GOF_c_hat","Open.tree.canopy","Herbaceous.grassland")])

Figure 10. Plot of abundance model-estimated beta coefficients for site-level covariates. Points and bars of mean and upper and lower 95% CIs are shown in purple for non-zero crossing and gray for zero crossing covariates.

#########################################################################
#Plot the beta coefficients 

  #get model summary
  mod.summary<-summary(mod.global)

## 
## Call:
## pcount(formula = ~Visit ~ Herbaceous.grassland + Open.tree.canopy, 
##     data = new.umf, se = TRUE)
## 
## Abundance (log-scale):
##                      Estimate    SE     z  P(>|z|)
## (Intercept)             0.421 0.137  3.08 2.08e-03
## Herbaceous.grassland   -0.515 0.378 -1.36 1.74e-01
## Open.tree.canopy        2.926 0.384  7.62 2.63e-14
## 
## Detection (logit-scale):
##              Estimate    SE      z  P(>|z|)
## (Intercept)    -0.473 0.173  -2.73 6.29e-03
## VisitVisit_2   -0.389 0.121  -3.22 1.30e-03
## VisitVisit_3   -1.712 0.162 -10.54 5.49e-26
## 
## AIC: 1841.402 
## Number of sites: 266
## optim convergence code: 0
## optim iterations: 31 
## Bootstrap iterations: 0

  new.coefs<-as.data.frame(mod.summary$state)
  new.coefs$Variable<-row.names(new.coefs)
  row.names(new.coefs)<-NULL
  
  #coompute 95% CIs
  new.coefs$`Lower_95%_CI` <- new.coefs$Estimate - 1.96*new.coefs$SE 
  new.coefs$`Upper_95%_CI` <- new.coefs$Estimate + 1.96*new.coefs$SE

#first remove intercepts
coefs.abun.save.sub<-subset(new.coefs, Variable != "(Intercept)")

#add ZeroCross column
coefs.abun.save.sub$ZeroCross<-ifelse(coefs.abun.save.sub$`Lower_95%_CI`>0 | coefs.abun.save.sub$`Upper_95%_CI` < 0,"No_Zero_Cross","Zero_Cross")

#make a factor
coefs.abun.save.sub$ZeroCross<-as.factor(coefs.abun.save.sub$ZeroCross)

#set factor levels for Variable
coefs.abun.save.sub$Variable<-factor(coefs.abun.save.sub$Variable, levels=unique(coefs.abun.save.sub$Variable))

#set colors
myBetaColors<-c("magenta","gray60")

  #plot beta coefficients
  betaCoefPlot_abun<-ggplot(data=coefs.abun.save.sub, aes(x=Variable, y=Estimate))+
        geom_hline(yintercept=0,linetype=2)+
        geom_errorbar(aes(ymin=`Lower_95%_CI`, ymax=`Upper_95%_CI`, color = ZeroCross),width=0)+
        geom_point(aes(color=ZeroCross),alpha=0.8)+
        scale_y_continuous(limits=c(-5000,5000))+
        #coord_cartesian(ylim=c(-180,180))+
        scale_color_manual(values=myBetaColors)+
        coord_flip(ylim=c(-100,100))+
#scale_color_gradient2(low="red",mid="gray40",high="royalblue2",midpoint = 0)+
theme(legend.position="none",
              panel.background = element_rect(fill="white",color="black"),
              panel.border = element_rect(fill="transparent",color="black"))+
    labs(y="Beta Coefficient Estimate", x="Site-level Covariate")
  
betaCoefPlot_abun

ggsave(betaCoefPlot_abun, file="/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/Fig_10_Beta_coefficient_plot_abundance.png",width=7, height=6,dpi=600)

Figure 11. Estimated relative abundance in relation to covariates.

#melt
# mod.predict.melt<-melt(mod.predict.sub, id.vars=c("PointID","Predicted","SE","lower","upper"), measure.vars=c("Elevation_m","Closed.tree.canopy","Herbaceous.grassland","Herbaceous.shrub.steppe","Open.tree.canopy","Open.Tree.Canopy","Shrubland","Sparse.tree.canopy","Sparsely.vegetated","Non.vegetated","No.Dominant.Lifeform"), variable.name = "Variable", value.name="value")

mod.predict.melt<-melt(mod.predict.sub, id.vars=c("PointID","Predicted","SE","lower","upper"), measure.vars=c("Herbaceous.grassland","Open.tree.canopy"), variable.name = "Variable", value.name="value")


  #plot Abundance vs. Vegetation class proportions
  abun.plot<-ggplot(data=mod.predict.melt, aes(x=value, y=Predicted))+
    geom_point(aes(color=Variable),alpha=0.1,size=1)+
    geom_smooth(method = "glm", formula = y~x, method.args = list(family = poisson(link = 'log')), aes(color=Variable, fill=Variable),alpha=0.2,fullrange=FALSE)+
    # scale_color_manual(values=myColors)+
    # scale_fill_manual(values=myColors)+
theme(
      panel.background = element_rect(fill="white",color="black"),
      panel.border = element_rect(fill="transparent",color="black"),
      axis.text.x = element_text(size=8),
      axis.text.y = element_text(size=8),
      legend.position="bottom"
    )+
    labs(x="Abundance Covariate", y="Predicted relative abundance")
  
  #facet over Vegetation_class
  abun.plot.facet<-abun.plot+facet_wrap(Variable~.,scales="free",ncol=4)
  abun.plot.facet

  #create filename with timestamp
  new.filename<-paste("Fig_11_Predicted_Abundance_models_with_veg_class",".png",sep="")
  
  ggsave(abun.plot.facet, file=paste("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/",new.filename, sep=""),width=8, height=8, dpi=600)

Use abundance model with landcover covs to get predicted estimates

#load predictive surface
pred_surface<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Data/ALL_covs_pred_surface.csv")

# #remove unwanted columns
#pred_surface$Dwarf.shrubland<-NULL
# pred_surface$Herbaceous.shrub.steppe<-NULL
# pred_surface$No.Dominant.Lifeform<-NULL
# pred_surface$Non.vegetated<-NULL
# pred_surface$Shrubland<-NULL
# pred_surface$Sparse.tree.canopy<-NULL
# pred_surface$Elevation_m<-NULL

#pred_surface$Sparsely.vegetated<-NULL


#get predicted estimates of avian occupancy as a function of vegetation cover on map
abun_pred <- predict(mod.global,
                    newdata = as.data.frame(pred_surface),
                    type = "state", appendData=TRUE)

#replace any estimates > max predicted value from model with max
# range(abun_pred$Predicted)
# range(abun_pred$SE)
# max.mod.value<-max(mod.predict.sub$Predicted)
# 
# abun_pred$Predicted<-ifelse(abun_pred$Predicted>max.mod.value, max.mod.value, abun_pred$Predicted)

suppressMessages({

# add to prediction surface
pred_abun <- bind_cols(pred_surface, 
                      abun_estimate = abun_pred$Predicted, 
                      abun_se = abun_pred$SE) %>% 
  select(Latitude, Longitude, abun_estimate, abun_se)

  })

 #melt to plot in ggplot
pred_abun.melt<-melt(pred_abun, id.vars=c("Latitude","Longitude"),
                    variable.name = "Measure",value.name="Relative_abundance")

Figures 12. Maps of predicted species abundance (mean ± SE) around DE Bay.

  map.abun.plot<-ggmap(basemap,darken = c(0.4, "black"))+
    geom_tile(data=pred_abun.melt, aes(x=Longitude, y=Latitude, fill=`Relative_abundance`,color=`Relative_abundance`))+
  geom_map(data=states.poly, map=states.poly, aes(y=lat, x=long, map_id=id),color=alpha("white",0.3),fill="transparent", inherit.aes = FALSE,size=0.1)+
    scale_fill_gradientn(colors=viridis(100),na.value = "transparent")+
    scale_color_gradientn(colors=viridis(100),na.value = "transparent")+
  labs(x="Longitude",y="Latitude")+
    theme(axis.text=element_text(size=12),
      axis.title=element_text(size=14),
      axis.text.x=element_text(size=8),
      axis.text.y=element_text(size=8),
      panel.border=element_rect(fill="transparent", color="black"),
      legend.position = "right")+
    ggtitle(paste("WPWI"," - Estimated relative abundance",sep=""))+
    #geom_text_repel(check_overlap = FALSE, size=3, nudge.y=-10, nudge.x=-10,aes(x=Longitude, y=Latitude,color=PatchID, label=PatchID),data=PatchID_labels)+
    #guides(fill = guide_legend(title = "National Wildlife Refuge", label = FALSE))+
    ggsn::scalebar(x.min=mapExtent$ll.lon, x.max=mapExtent$ur.lon, y.min=mapExtent$ll.lat, y.max=mapExtent$ur.lat,transform=TRUE, dist_unit="km",dist=100, st.dist=0.03, height = 0.015, anchor=c(x=-70.5,y=38.55), border.size=0.5,st.size=3)

 map.abun.plot.facet.save<-map.abun.plot+facet_wrap(Measure~., ncol=1)
  
  #create filename with timestamp
  new.filename<-paste(paste("Fig_12_WPWI_","Map_Predicted_relative_abundance",sep=""),".png",sep="")
    
  ggsave(map.abun.plot.facet.save, file=paste("/Users/zach/Dropbox (ZachTeam)/Projects/USFWS_WPWI/Results/Figures/",new.filename, sep=""),width=14, height=8, dpi=600)
  
print(map.abun.plot.facet.save)

References