Predicting the Probability of Occupancy of Shorebirds in Delaware Bay, U.S.A.

Overview

The Delaware Bay in Delaware and New Jersey, U.S.A. provides substantial tidal marsh and mudflat habitat for many shorebird species of conservation concern. These include the American Oystercatcher (Haematopus palliatus), Greater Yellowlegs (Tringa melanoleuca), Lesser Yellowlegs (Tringa flavipes), Piping Plover (Charadrius melodus), Ruddy Turnstone (Arenaria interpres), Sanderling (Calidris alba), and Semipalmated Sandpiper (Calidris pusilla).

Throughout, species are indicated by corresponding 4-letter American Ornithological Society codes:

• American Oystercatcher = AMOY
• Greater Yellowlegs = GRYE
• Lesser Yellowlegs = LEYE
• Ruddy Turnstone = RUTU
• Semipalmated Sandpiper = SESA

Sea-level rise is altering tidal marsh hydrology and reducing habitat quality for these species. We analyzed extensive bird monitoring data to estimate the probability of occupancy of shorebird species. We can use the occupancy models to identify high quality areas, high preference areas, or areas of high importance. and see which habitat types are most important to each species and to shorebirds and the tidal bird community as a whole. Understanding the factors that influence the estimated occupancy for these restricted and declining species in Delaware Bay will provide land managers and restoration practitioners with more detailed information to identify sites that are currently most important and sites that are potentially the most important in the future for these species.

Purpose

The goal of this exploratory analysis is to use the existing Saltmarsh and Habitat Avian Research Program (SHARP; https://www.tidalmarshbirds.org/) tidal marsh bird survey data in Delaware Bay to estimate the probability of site occupancy for focal shorebird species of conservation concern with sufficient detection data, and explore the factors that may influence these population parameters.

Objectives

1. Estimate focal species probability of site occupancy
   
2. Model the factors using remotely sensed information (e.g., vegetation, sea-level rise, land cover) that influence the probability of site occupancy for the focal species
   
3. Identify priority tidal marsh habitat and sites for focal species in Delaware Bay

Methods

We used existing tidal marsh bird survey data from 2011-2020 collected at > 300 locations in Delaware Bay tidal marshes to estimate occupancy for six focal shorebird species (SHARP 2017). Trained observers visited each location at least twice during each breeding season (May – Aug) and counted the number of individuals for detected during specified time intervals and distance bands. The field sampling protocol (Wiest et al. 2016) we used allows for the use of hierarchical models that can simultaneously estimate detection probability to adjust the occupancy estimates (Fiske and Chandler 2011, Ladin et al. 2020). Our occupancy estimates are spatially explicit so that we can model the factors that influence population parameters using remotely sensed data. We used a high resolution tidal marsh vegetation layer (Correll et al. 2019), a tidal marsh ditching index, a tidal marsh restriction index, and sea-level rise estimates (NOAA 2020) as independent variables in our models. We used a model selection approach to determine the most important covariates for the occupancy of each of the focal species and predict these parameters throughout the tidal marsh habitat in Delaware Bay. We then generated maps for the predicted occupancy for each focal species to help identify priority sites for conservation and management.

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(sf)

# resolve namespace conflicts
select <- dplyr::select
map <- purrr::map
projection <- raster::projection
  
#set working directory
opts_knit$set(root.dir = "./")

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

#load functions
source("./Source/makePOccu_umf.R")
source("./Source/makePOccu_umf_with_covs.R")
source("./Source/makePCount_umf.R")
source("./Source/makePCount_umf_with_covs.R")
source("./Source/AICcmodavg.gof.test_source.R")
source("./Source/xyToLatLong.R")

})

Load SHARP point count data, subset, and export as .csv

#load SHARP data (be carefull. . .this takes a while to load!)
sharp_master<-read_xlsx("./Data/SHARP Survey Data Master File 20200917_callback_ok.xlsx", sheet="vw_SurveyData")

# #subset data
callback_DE_NJ<-subset(sharp_master, c(State=="DE" | State=="NJ"))

# #rename POINT_X and POINT_Y as Longitude and Latitude
names(callback_DE_NJ)[names(callback_DE_NJ)=="Point_X"]<-"Longitude"
names(callback_DE_NJ)[names(callback_DE_NJ)=="Point_Y"]<-"Latitude"

# #get list locations with missing coords
missing_coords_locations<-subset(unique(callback_DE_NJ[,c("PointID","Longitude","Latitude")]), is.na(Longitude))
#0

# #remove all points with Longitude > -74.5
callback_DE_NJ_sub<-subset(callback_DE_NJ, Longitude < -74.5)

#make PatchID a factor
callback2020_DE_NJ$PatchID<-as.factor(callback2020_DE_NJ$PatchID)

# #make PatchID_labels to see which patches to keep on map
PatchID_labels_long<-aggregate(Longitude~PatchID, FUN=mean, data=callback2020_DE_NJ)
PatchID_labels_lat<-aggregate(Latitude~PatchID, FUN=mean, data=callback2020_DE_NJ)
PatchID_labels<-merge(PatchID_labels_lat, PatchID_labels_long, by="PatchID")

# #get mean Lat/Long to center basemap
meanLong<-mean(callback_DE_NJ_sub$Longitude, na.rm=TRUE)
meanLat<-mean(callback_DE_NJ_sub$Latitude, na.rm=TRUE)
#centroid of DE Bay
meanLong<--75.250147
meanLat<-39.133269

# #get basemap initially to determine ballpark extent
basemap<-get_map(location=c(meanLong, meanLat),zoom=09, maptype="terrain")
#ggmap(basemap)

# #get map extent
mapExtent<-attr(basemap, "bb")
#38.80475   -75.62104   39.48637    -74.74214

#refine map extent to include all sampling locations around DE Bay
basemap<-get_map(location=c(-75.7,38.7,-74.5,39.7),maptype="terrain")

#plot sampling locations over study area basemap
studyAreaMap<-ggmap(basemap)+
  geom_point(data=callback_DE_NJ_sub, aes(x=Longitude, y=Latitude,color=PatchID), shape=16,size=1, alpha=0.4)+
  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("Delaware Bay SHARP point count locations")+
  #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=15, st.dist=0.04, height = 0.015, anchor=c(x=-74.5,y=38.6), border.size=0.5,st.size=4)

studyAreaMap

# #make PatchID a factor
callback2020_DE_NJ$PatchID<-as.factor(as.character(callback2020_DE_NJ$PatchID))

PatchRemoveList<-c("9121","9105", "8938","8878", "8572","8511","8423","8392","8259", "8162","7585","7257","7239","7221","7220","7186","7118","6181","6153")

# #subset points in NJ that aren't adjacent to DE bay
callback2020_DE_NJ_sub<-subset(callback2020_DE_NJ, c(Longitude < -74.75 & ! PatchID %in% PatchRemoveList))

#save example data
write.csv(callback_DE_NJ_sub, file="./Data/Callback_DE_NJ_sub_wide.csv",row.names=FALSE)

Results

Figure 1. Map of study area showing tidal marshes and sampling locations (colored by PatchID) around Delaware Bay.

#read in subset data
callback2020_DE_NJ_sub<-read.csv("./Data/Callback_DE_NJ_sub_wide.csv")

#get mean Lat/Long to center basemap
meanLong<-mean(callback2020_DE_NJ_sub$Longitude, na.rm=TRUE)
meanLat<-mean(callback2020_DE_NJ_sub$Latitude, na.rm=TRUE)
#centroid of DE Bay
# meanLong<--75.250147
# meanLat<-39.133269

#get basemap initially to determine ballpark extent
basemap<-get_map(location=c(meanLong, meanLat),zoom=9, maptype="terrain")
#ggmap(basemap)

#get map extent
mapExtent<-attr(basemap, "bb")
#38.80475   -75.62104   39.48637    -74.74214

#refine map extent to include all sampling locations around DE Bay
#basemap<-get_map(location=c(-75.7,38.7,-74.5,39.7),maptype="terrain")

#Make PatchID a factor
callback2020_DE_NJ_sub$PatchID<-as.factor(callback2020_DE_NJ_sub$PatchID)

#plot sampling locations over study area basemap
studyAreaMap<-ggmap(basemap, darken = c(0.5, "black") )+
  geom_point(data=callback2020_DE_NJ_sub, aes(x=Longitude, y=Latitude,color=PatchID), shape=16,size=0.3, alpha=0.3)+
  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 = "none")+
  ggtitle("Delaware Bay SHARP point count locations")+
  #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=15, st.dist=0.04, height = 0.015, anchor=c(x=-74.5,y=38.55), border.size=0.5,st.size=4)

studyAreaMap

#save map
ggsave(studyAreaMap, filename="./Results/Figures/Fig_1_DE_Bay_study_area_map.png",width=8, height=6, dpi=600)

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

#List oof AFSI focal species: American Oystercatcher, American Golden-Plover, Snowy Plover, Wilson's Plover, Marbled Godwit, Ruddy Turnstone, Red Knot, Sanderling, Piping Plover, Greater Yellow-legs, Lesser Yellowlegs, Whimbrel, Purple Sandpiper, Semipalmated Sandpiper, Red-necked Phalarope, Eskimo Curlew

#get speciesList
speciesList<-c("AMOY","GRYE", "LEYE", "RUTU","SESA")

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

#for loop to iterate over species and estimate probability of occupancy
for(i in 1:length(speciesList)){
  
  new.AlphaCode<-speciesList[i]
  
  #use functions to build unmarkedFrame objects
  occu.umf<-makePOccu_umf(DataIn=callback2020_DE_NJ_sub, SpeciesIn=new.AlphaCode)

  #check count data
  #occu.umf.test<-occu.umf@y
  
  #fit example model to test
  mod<-NULL
  mod<-occu(~Visit ~1, data=occu.umf)
  
  suppressMessages({
  #run GOF test
  #mod.gof<- try(mb.gof.test(mod, nsim = 10, 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
  mod.predict<-unique(predict(mod, type="state",appendData=FALSE))
  
  })
  #add species
  mod.predict$AlphaCode<-new.AlphaCode
  
  #add GOF statisitc
  mod.predict$GOF_c.hat<-mod.c.hat
  
  #organize columns
  mod.predict<-mod.predict[,c("AlphaCode","Predicted","SE","lower","upper","GOF_c.hat")]
  colnames(mod.predict)<-c("Alpha Code","Mean p(Occupancy)","SE p(Occupancy)","Lower 95%CI","Upper 95%CI","GOF: c-hat")
  
  occu.save<-rbind(occu.save, mod.predict)

}

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

pander(occu.save)

Alpha Code	Mean p(Occupancy)	SE p(Occupancy)	Lower 95%CI	Upper 95%CI	GOF: c-hat
AMOY	0.04376	0.005574	0.03404	0.05608	NA
GRYE	0.2848	0.02377	0.2405	0.3335	NA
LEYE	0.3082	0.0628	0.2	0.4424	NA
RUTU	0.02624	0.006478	0.01613	0.04241	NA
SESA	0.3146	0.03266	0.2543	0.3818	NA

Gather and format land cover and spatial covariates for modeling

List of covariates at each sampling location:

1. Patch area (km^2)
   
2. Distance to water (m)
   
3. Distance to ocean (m)
   
4. Distance to upland vegetation (m)
   
5. Proportion of vegetation class (within 100-m radius buffer)
   
6. Density of marsh ditching
   
7. Tidal restriction effect (0/1)

Format vegetation data layer (3-m resolution)

#load vegetation raster layer
# veg.raster<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/High Marsh Layer/TidalMarshVegClass/TidalMarshVegClass/Regional_DEM.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 UTM
basemap_extent_UTM <- spTransform(basemap_extent_sp, projection(veg.raster))

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

#bring in Patch polygons
# sharp_patches<-readOGR(dsn="/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/Patches/", layer="SHARP_Patches_SR1-9_w_bird_dens")

#get PatchList
PatchList<-unique(as.character(callback2020_DE_NJ_sub$PatchID))

#subset sharp_patches .shp
sharp_patches_sub<-subset(sharp_patches, PatchID %in% PatchList)

#plot(veg.raster_crop)
#plot(sharp_patches_sub, col="blue")

#projection(sharp_patches_sub)

#reproject sharp_patches_sub
sharp_patches_sub_proj<-spTransform(sharp_patches_sub, CRS(projection(veg.raster_crop)))

#dissolve sharp_patches_sub_proj into one polygon, but first fix self-intersection regions with gBuffer, width=0
sharp_patches_sub_proj_valid<-gBuffer(sharp_patches_sub_proj, byid=TRUE, width=0)
#dissolve using sp::aggregate
sharp_patches_sub_proj_dissolve<-aggregate(sharp_patches_sub_proj_valid, dissolve = TRUE)

#plot(sharp_patches_sub_proj_dissolve)

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

#plot(sharp_patches_sub_proj_dissolve_buffer)

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

#plot(veg.raster.mask)

#save cropped raster
writeRaster(veg.raster_crop, file="./Data/raster_data/TidalMarshVegClass_UTM_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="./Data/raster_data/TidalMarshVegClass_UTM_crop.tif",
                    dstfile=file.path("./Data/raster_data/TidalMarshVegClass_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="./Data/raster_data/TidalMarshVegClass_UTM_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="./Data/raster_data/TidalMarshVegClass_UTM_mask.tif",
                    dstfile=file.path("./Data/raster_data/TidalMarshVegClass_latlong_mask.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",overwrite=TRUE,verbose=TRUE)

Figure 2. Tidal marsh vegetation class raster layer.

#read in cropped raster in lat/long
veg.raster.crop.latlong<-raster("./Data/raster_data/TidalMarshVegClass_latlong_crop.tif")

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

rat <- levels(r)[[1]]
rat$Vegetation_class <- c("High_Marsh","Low_Marsh", "Mudflat","Phragmites",
                          "Pool_Pan","Stream","Terrestrial_Border","Upland")
rat$band <- c(1,2,4,5,6,7,8,9)
levels(r) <- rat

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

# plot(r,col=myColors,legend=FALSE)
# legend("bottomright", legend=c("High Marsh","Low Marsh", "Mudflat","Phragmites",
#                         "Pool/Pan","Stream","Terrestrial Border","Upland"),fill=myColors)

#save to not have to evaluate later (takes long time!)
png(filename="./Results/Figures/Fig_2_TidalMarsh_veg_layer.png",width=8, height=6, units="in", res=600)

  plot(r,col=myColors,legend=FALSE)
  legend("bottomright", legend=c("High Marsh","Low Marsh", "Mudflat","Phragmites",
                          "Pool/Pan","Stream","Terrestrial Border","Upland"),fill=myColors)
dev.off()

#Load the saved image of raster veg layer
include_graphics("./Results/Figures/Fig_2_TidalMarsh_veg_layer.png")

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

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

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

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

#plot(buffers.polygon)

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

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

#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+Vegetation_class, 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$Veg_class_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)

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

#save .csv
write.csv(count.data.merge, file="./Data/Callback_2020_DE_NJ_withVegCovs.csv",row.names=FALSE)

Figure 3. Map showing line used to compute distance from each sampling point to ocean.

#add spatial line to delineate "ocean" from DE bay (Add points to extend line south along DE inland bays)
#oceanLine<-data.frame(x=c(-75.094463, -74.933379), y=c(38.805033,38.930453))

oceanLine<-data.frame(x=c(-75.0, -75.0, -74.4), 
                      y=c(38.47, 38.76, 39.3))


#plot line on map to make sure it's lined up
#plot sampling locations over study area basemap
studyAreaMap_oceanLine<-ggmap(basemap, darken=c(0.5, "black"))+
  geom_point(data=callback2020_DE_NJ_sub, aes(x=Longitude, y=Latitude,color=PatchID), shape=16,size=0.3, alpha=0.3)+
  geom_path(data=oceanLine, aes(x=x, y=y),color="red",linetype=2)+
  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 = "none")+
  ggtitle("Delaware Bay SHARP point count locations")+
  #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=15, st.dist=0.04, height = 0.015, anchor=c(x=-74.5,y=38.55), border.size=0.5,st.size=4)

studyAreaMap_oceanLine

#convert to spatialpoints, and change CRS to utm
oceanLine.pts <- sp::SpatialPoints(cbind(oceanLine$x,oceanLine$y), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
#convert to utm
oceanLine.pts.utm<-spTransform(oceanLine.pts, CRS("+init=epsg:26918"))

# use as to convert to line
oceanLine.utm <- as(oceanLine.pts.utm,"SpatialLines")
#crs(oceanLine.sp_line)

#make points into spatialPoints (put into UTM coords)
sampling_pts.sp <- sp::SpatialPoints(unique(cbind(callback2020_DE_NJ_sub$Longitude,callback2020_DE_NJ_sub$Latitude)), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

#convet to utm
sampling_pts.utm <-spTransform(sampling_pts.sp, CRS("+init=epsg:26918"))

#now compute distance from each point to oceanLine
dist_to_ocean_m<-data.frame("PointID" = unique(callback2020_DE_NJ_sub$PointID),
                                "Distance"= gDistance(oceanLine.utm, sampling_pts.utm, byid=TRUE))
names(dist_to_ocean_m)[names(dist_to_ocean_m)=="ID"]<-"Distance"

range(dist_to_ocean_m$Distance)

## [1]  4734.865 98060.006

#1363.037 91708.477

#save distance to ocean as .csv file
write.csv(dist_to_ocean_m, file="./Data/Dist_to_ocean_m.csv",row.names=FALSE)

Compute distance from each sampling to nearest water.

#convert to points
r.stream.points<-rasterToPoints(r, function(x) x == 7)
r.stream.points.df<-unique(as.data.frame(r.stream.points))

#convert points to utm
#make points into spatialPoints (put into UTM coords)
r.stream.points.sp <- sp::SpatialPoints(cbind(r.stream.points.df$x,r.stream.points.df$y), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

plot(r.stream.points.sp)

#convet to utm
r.stream.points.utm <-spTransform(r.stream.points.sp, CRS("+init=epsg:26918"))

#convert to spatialPolygon
r.stream.points.utm.df<-as(r.stream.points.utm,"SpatialPointsDataFrame")
r.stream.poly<-Polygon(r.stream.points.utm.df)

r.stream.poly.2 <- sp::Polygons(list(r.stream.poly), ID = "ID")
r.stream.poly.sp <- sp::SpatialPolygons(list(r.stream.poly.2))

#now compute distance from each point to oceanLine
dist_to_water_m<-unique(data.frame("PointID" = unique(callback2020_DE_NJ_sub$PointID),
                                "Distance"= gDistance(r.stream.poly.sp, sampling_pts.utm, byid=TRUE)))
names(dist_to_water_m)[names(dist_to_water_m)=="ID"]<-"Distance"

#range(dist_to_water_m$Distance)
#0.0000 807.6468

#save distance to ocean as .csv file
write.csv(dist_to_water_m, file="./Data/Dist_to_water_m.csv",row.names=FALSE)

Compute distance from each sampling to nearest upland land cover.

#convert to points
r.upland.points<-rasterToPoints(r, function(x) x == 9)
r.upland.points.df<-unique(as.data.frame(r.upland.points))

#convert points to utm
#make points into spatialPoints (put into UTM coords)
r.upland.points.sp <- sp::SpatialPoints(cbind(r.upland.points.df$x,r.upland.points.df$y), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

#convet to utm
r.upland.points.utm <-spTransform(r.upland.points.sp, CRS("+init=epsg:26918"))

#convert to spatialPolygon
r.upland.points.utm.df<-as(r.upland.points.utm,"SpatialPointsDataFrame")
r.upland.poly<-Polygon(r.upland.points.utm.df)

r.upland.poly.2 <- sp::Polygons(list(r.upland.poly), ID = "ID")
r.upland.poly.sp <- sp::SpatialPolygons(list(r.upland.poly.2))

#now compute distance from each point to oceanLine
dist_to_upland_m<-unique(data.frame("PointID" = unique(callback2020_DE_NJ_sub$PointID),
                                "Distance"= gDistance(r.upland.poly.sp, sampling_pts.utm, byid=TRUE)))
names(dist_to_upland_m)[names(dist_to_upland_m)=="ID"]<-"Distance"

#range(dist_to_upland_m$Distance)
#0.000 992.268

#save distance to ocean as .csv file
write.csv(dist_to_upland_m, file="./Data/Dist_to_upland_m.csv",row.names=FALSE)

Extract marsh ditch density

#read in marsh ditch density raster
# marsh.ditch.density<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_ditching_Metric_2017/DSL_ditches_2010_v3.1/ditchgrid_2010_v3.1.tif")


#get projection of marsh.ditch.grid
crs(marsh.ditch.density)
# +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 

#reproject marsh.ditch.density to latlong (use gdalwarp for faster progress)
#marsh.ditch.density.latlong<-projectRaster(from=marsh.ditch.density, crs=crs(r))

# gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_ditching_Metric_2017/DSL_ditches_2010_v3.1/ditchgrid_2010_v3.1.tif",
#                     dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_ditching_Metric_2017/DSL_ditches_2010_v3.1/ditchgrid_2010_v3.1_latlong.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",overwrite=TRUE,verbose=TRUE)
# 
# marsh.ditch.density.latlong<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_ditching_Metric_2017/DSL_ditches_2010_v3.1/ditchgrid_2010_v3.1_latlong.tif")

#plot
marsh.ditch.density.latlong.crop<-crop(marsh.ditch.density.latlong, extent(r))
plot(marsh.ditch.density.latlong.crop, col=rev(ocean.solar(100)), main="Marsh Ditch Density")

#extract values (IDs are circular buffer polygons) (ASK Greg and Liz. . .what are units of these ditches per km^2?)
extracted.ditch.density<-raster::extract(marsh.ditch.density.latlong.crop, buffers.polygon, factors=FALSE, df=TRUE, fun=mean, na.rm=TRUE)
names(extracted.ditch.density)[names(extracted.ditch.density)=="ditchgrid_2010_v3.1_latlong"]<-"Ditch_density"

#simplify to get only unique rows
extracted.ditch.density.2<-unique(extracted.ditch.density)

#fill in 0s where NaN occur
extracted.ditch.density.2$Ditch_density<-ifelse(is.na(extracted.ditch.density.2$Ditch_density), 0, extracted.ditch.density.2$Ditch_density)

#make ID and Vegetation_class factors
extracted.ditch.density.2$ID<-as.factor(extracted.ditch.density.2$ID)

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

#save .csv
write.csv(ditch.density.df, file="./Data/Marsh_ditch_density.csv",row.names=FALSE)

Get tidal restriction data.

#read in tidal restriction raster (I used the 2020 data)
# tidal.restriction<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_tiderestrict_2010_v3.0/DSL_data_tidal_restrictions/tideres_2020_v5.0.tif")

#get projection of marsh.ditch.grid
crs(tidal.restriction)
# +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 

# gdalUtils::gdalwarp(srcfile="/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_tiderestrict_2010_v3.0/DSL_data_tidal_restrictions/tideres_2020_v5.0.tif",
#                     dstfile=file.path("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_tiderestrict_2010_v3.0/DSL_data_tidal_restrictions/tideres_2020_v5.0_latlong.tif"),t_srs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",overwrite=TRUE,verbose=TRUE)
# 
# tidal.restriction.latlong<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_tiderestrict_2010_v3.0/DSL_data_tidal_restrictions/tideres_2020_v5.0_latlong.tif")

#plot
tidal.restriction.latlong.crop<-crop(tidal.restriction.latlong, extent(r))
plot(tidal.restriction.latlong.crop, col=ocean.oxy(100),main="Tide Restriction Density")

#extract values (IDs are circular buffer polygons) (ASK Greg and Liz. . .what are units of these tidal restrictions per km^2?)
extracted.tidal.restriction.density<-raster::extract(tidal.restriction.latlong.crop, buffers.polygon, factors=FALSE, df=TRUE, fun=mean, na.rm=TRUE)
names(extracted.tidal.restriction.density)[names(extracted.tidal.restriction.density)=="tideres_2020_v5.0_latlong"]<-"Tidal_Restriction"

#make ID and Vegetation_class factors
extracted.tidal.restriction.density$ID<-as.factor(extracted.tidal.restriction.density$ID)

#fill in 0s where NaN occur
extracted.tidal.restriction.density$Tidal_Restriction<-ifelse(is.na(extracted.tidal.restriction.density$Tidal_Restriction), 0, extracted.tidal.restriction.density$Tidal_Restriction)


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

#save .csv
write.csv(extracted.tidal.restriction.density.df, file="./Data/Tidal_restriction_density.csv",row.names=FALSE)

Get patch area info.

#read in patch shapefile
# sharp_patches<-readOGR(dsn="/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/Patches/", layer="SHARP_Patches_SR1-9_w_bird_dens")

#get PatchList
PatchList<-unique(as.character(callback2020_DE_NJ_sub$PatchID))

#subset sharp_patches .shp
sharp_patches_sub<-subset(sharp_patches, PatchID %in% PatchList)

#plot(sharp_patches_sub, col="blue")

#reproject sharp_patches_sub
sharp_patches_sub_proj<-spTransform(sharp_patches_sub, CRS(projection(new.raster)))

#create raster as a function of patch area
patch.area.raster<- rasterize(sharp_patches_sub_proj, r,field="area_ha")

plot(patch.area.raster, col=parula(100))

#save raster
writeRaster(patch.area.raster, filename="./Data/raster_data/Patch_area_raster_latlong.tif",overwrite=TRUE)

#now get Patch attributes to include as site-level covariates in models
# patch.attributes<-read.csv("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/Patches/Patch_attributes.csv")

#simplify
patch.attributes.sub<-unique(patch.attributes[,c("PatchID","area_ha")])

#save
write.csv(patch.attributes.sub, file="./Data/Patch_area_info.csv",row.names=FALSE)

Add additional covariates to count data.

#read in new data with Veg layer covariates
callback2020_DE_NJ_veg<-unique(read.csv("./Data/Callback_2020_DE_NJ_withVegCovs.csv"))

#load patch area
patch.df<-read.csv("./Data/Patch_area_info.csv")

#add other covs
covs.merge.1<-merge(callback2020_DE_NJ_veg, patch.df, by="PatchID",all.x=TRUE)

#check for missing values
#area_check<-subset(covs.merge.1, is.na(area_ha))

#load dist to ocean
ocean.df<-read.csv("./Data/Dist_to_ocean_m.csv")
names(ocean.df)[names(ocean.df)=="Distance"]<-"Dist_to_ocean_m"

#add other covs
covs.merge.2<-merge(covs.merge.1, ocean.df, by="PointID",all.x=TRUE)

#check for missing values
dist_to_ocean_check<-subset(covs.merge.2, is.na(Dist_to_ocean_m))

#load dist to water
water.df<-read.csv("./Data/Dist_to_water_m.csv")
names(water.df)[names(water.df)=="Distance"]<-"Dist_to_water_m"

#add other covs
covs.merge.3<-merge(covs.merge.2, water.df, by="PointID",all.x=TRUE)

#check for missing values
dist_to_water_check<-subset(covs.merge.3, is.na(Dist_to_water_m))

#load dist to upland
upland.df<-read.csv("./Data/Dist_to_upland_m.csv")
names(upland.df)[names(upland.df)=="Distance"]<-"Dist_to_upland_m"

#add other covs
covs.merge.4<-merge(covs.merge.3, upland.df, by="PointID",all.x=TRUE)

#check for missing values
dist_to_upland_check<-subset(covs.merge.4, is.na(Dist_to_upland_m))

#load marsh ditch density
marsh.ditch.df<-read.csv("./Data/Marsh_ditch_density.csv")
marsh.ditch.df$ID<-as.factor(marsh.ditch.df$ID)

#add other covs
covs.merge.5<-merge(covs.merge.4, marsh.ditch.df, by=c("PointID","Longitude","Latitude"),all.x=TRUE)
#names(covs.merge.5)

#check for missing values
marsh_ditch_check<-subset(covs.merge.5, is.na(Ditch_density))


#load tidal restriction density
tidal.res.df<-read.csv("./Data/Tidal_restriction_density.csv")

#add other covs
covs.merge.6<-merge(covs.merge.5, tidal.res.df, by=c("PointID","ID","Longitude","Latitude"),all.x=TRUE)

#check for missing values
tidal_restrict_check<-subset(covs.merge.6, is.na(Tidal_Restriction))

#save covs.merge.6
write.csv(covs.merge.6, file="./Data/Callback_2020_DE_NJ_with_ALL_Covs.csv",row.names=FALSE)

Build prediction surface and save as .csv

#build prediction surface (and save .csv)

#create a new blank raster with all values = 1 with dimensions for prediction

veg.raster.crop.latlong<-raster("./Data/raster_data/TidalMarshVegClass_latlong_crop.tif")

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

rat <- levels(r)[[1]]
rat$Vegetation_class <- c("High_Marsh","Low_Marsh", "Mudflat","Phragmites",
                          "Pool_Pan","Stream","Terrestrial_Border","Upland")
rat$band <- c(1,2,4,5,6,7,8,9)
levels(r) <- rat

# plot(r)
# dim(r)

#define neighborhood radius (average 300x300 cell regions)
#neighborhood_radius <-  0.00463 * ceiling(max(res(r))) / 2
#define neighborhood radius (average 420x420 cell regions)
neighborhood_radius <-  0.00463 * ceiling(max(res(r))) / 2 * 0.5102041

#aggregate among veg class factor levels
agg_factor <- round(2 * neighborhood_radius / res(r))

#downsample veg.raster to 10-m resolution
new.raster<- veg.raster.crop.latlong %>%
  raster::aggregate(agg_factor)

#round new.raster values
values(new.raster)<-round(values(new.raster),0)

#check raster output
#plot(new.raster)
#dim(new.raster)

new.raster.df<-as.data.frame(rasterToPoints(new.raster))
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))

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

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

##############################################################################
#Vegetation Class proportions

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

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

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

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

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

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

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

#Now merge back with sampling locs
pred.veg.proportions.df<-merge(pred.veg.merge, 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", "Vegetation_class"), measure.vars=c("Veg_class_proportion"))

#cast and fill any NAs with 0s
pred.veg.proportions.df.cast<-dcast(ID+Longitude+Latitude~Vegetation_class,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="./Data/Predictive_surface/Vegetation_pred_surface.csv", row.names=FALSE)

write.csv(pred.veg.proportions.df.cast, file="./Data/Predictive_surface/Vegetation_pred_surface_hi_res.csv", row.names=FALSE)

#######################################################################################################
#Marsh ditch density

# ditch.raster<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_ditching_Metric_2017/DSL_ditches_2010_v3.1/ditchgrid_2010_v3.1_latlong.tif")

#crop to landcover extent (bew.raster)
ditch.raster.crop<-crop(ditch.raster, extent(new.raster))
plot(ditch.raster.crop, col=parula(100))

#extract values (IDs are circular buffer polygons)
ditch.raster.crop.df<-raster::extract(ditch.raster.crop, pred.buffers.polygon, factors=FALSE, df=TRUE, fun=mean,na.rm=TRUE)

#add Latitude and Longitude to data.frame
#dist.raster.crop.df.2<-cbind(dist.raster.crop.df, pred.sampling.locs[,c("Longitude","Latitude")])
ditch.raster.crop.df.2<-merge(ditch.raster.crop.df, pred.sampling.locs, by=c("ID"),all.x=TRUE)

#make ID (cell location ID) factor
ditch.raster.crop.df.2$ID<-as.factor(ditch.raster.crop.df.2$ID)

#modify column names
names(ditch.raster.crop.df.2)[names(ditch.raster.crop.df.2)=="ditchgrid_2010_v3.1_latlong"]<-"Ditch_density"

#fill NAs with 0s
ditch.raster.crop.df.2$Ditch_density<-ifelse(is.na(ditch.raster.crop.df.2$Ditch_density),0,ditch.raster.crop.df.2$Ditch_density)

# write.csv(ditch.raster.crop.df.2, file="./Data/Predictive_surface/Marsh_ditch_density_pred_surface.csv", row.names=FALSE)

write.csv(ditch.raster.crop.df.2, file="./Data/Predictive_surface/Marsh_ditch_density_pred_surface_hi_res.csv", row.names=FALSE)

#######################################################################################################
#Tidal restrictions
#read in tidal restriction density raster
# tidal.raster<-raster("/Users/zach/Dropbox (ZachTeam)/Projects/SHARP_UD/gis/DSL_tiderestrict_2010_v3.0/DSL_data_tidal_restrictions/tideres_2020_v5.0_latlong.tif")

#crop to landcover extent (bew.raster)
tidal.raster.crop<-crop(tidal.raster, extent(new.raster))
plot(tidal.raster.crop, col=parula(100))

#extract values (IDs are circular buffer polygons)
tidal.raster.crop.df<-raster::extract(tidal.raster.crop, pred.buffers.polygon, factors=FALSE, df=TRUE, fun=mean,na.rm=TRUE)

#add Latitude and Longitude to data.frame
tidal.raster.crop.df.2<-merge(tidal.raster.crop.df, pred.sampling.locs, by=c("ID"),all.x=TRUE)

#make ID (cell location ID) factor
tidal.raster.crop.df.2$ID<-as.factor(tidal.raster.crop.df.2$ID)

#modify column names
names(tidal.raster.crop.df.2)[names(tidal.raster.crop.df.2)=="tideres_2020_v5.0_latlong"]<-"Tidal_restriction_density"

#fill NAs with 0s
tidal.raster.crop.df.2$Tidal_restriction_density<-ifelse(is.na(tidal.raster.crop.df.2$Tidal_restriction_density),0,tidal.raster.crop.df.2$Tidal_restriction_density)

# write.csv(tidal.raster.crop.df.2, file="./Data/Predictive_surface/Tidal_restriction_density_pred_surface.csv", row.names=FALSE)

write.csv(tidal.raster.crop.df.2, file="./Data/Predictive_surface/Tidal_restriction_density_pred_surface_hi_res.csv", row.names=FALSE)

##################################################################################
#patch area

#read in tidal restriction density raster
area.raster<-raster("./Data/raster_data/Patch_area_raster_latlong.tif")

#crop to landcover extent (bew.raster)
area.raster.crop<-crop(area.raster, extent(new.raster))
plot(area.raster.crop, col=parula(100), main="Patch Area (ha)")

#extract values (IDs are circular buffer polygons)
area.raster.crop.df<-raster::extract(area.raster.crop, pred.buffers.polygon, factors=FALSE, df=TRUE, fun=max)

#add Latitude and Longitude to data.frame
area.raster.crop.df.2<-merge(area.raster.crop.df, pred.sampling.locs, by=c("ID"),all.x=TRUE)

#make ID (cell location ID) factor
area.raster.crop.df.2$ID<-as.factor(area.raster.crop.df.2$ID)

#modify column names
names(area.raster.crop.df.2)[names(area.raster.crop.df.2)=="Patch_area_raster_latlong"]<-"Patch_area_ha"

#fill NAs with 0s
area.raster.crop.df.2$Patch_area_ha<-ifelse(is.na(area.raster.crop.df.2$Patch_area_ha),0,area.raster.crop.df.2$Patch_area_ha)


# write.csv(area.raster.crop.df.2, file="./Data/Predictive_surface/Patch_area_pred_surface.csv", row.names=FALSE)

write.csv(area.raster.crop.df.2, file="./Data/Predictive_surface/Patch_area_pred_surface_hi_res.csv", row.names=FALSE)


#######################################################################################################
#Distance to water

#create new blank raster with same dimensions as new.raster
new.raster.blank<-new.raster
new.raster.blank[!is.na(new.raster.blank)]<-1

#water_distances <- distanceFromPoints(object = new.raster.blank, xy = r.stream.points.sp)
water.raster<-new.raster
water.raster[water.raster!=7]<-NA

water.raster.distance<-distance(water.raster)
#plot(water.raster.distance, col=cubicyf(100),main="Distance to Water (m)")

#convert raster to data.frame and rename x and y
water.raster.distance.df<-as.data.frame(water.raster.distance,xy=TRUE)
water.raster.distance.df$ID<-row.names(water.raster.distance.df)
names(water.raster.distance.df)[names(water.raster.distance.df)=="x"]<-"Longitude"
names(water.raster.distance.df)[names(water.raster.distance.df)=="y"]<-"Latitude"
names(water.raster.distance.df)[names(water.raster.distance.df)=="layer"]<-"Dist_to_water_m"

# write.csv(water.raster.distance.df, file="./Data/Predictive_surface/Dist_to_water_pred_surface.csv", row.names=FALSE)

write.csv(water.raster.distance.df, file="./Data/Predictive_surface/Dist_to_water_pred_surface_hi_res.csv", row.names=FALSE)

#######################################################################################################
#Dist to upland
upland.raster<-new.raster
upland.raster[upland.raster!=9]<-NA

upland.raster.distance<-distance(upland.raster)
plot(upland.raster.distance, col=cubicyf(100),main="Distance to Upland (m)")

#convert raster to data.frame and rename x and y
upland.raster.distance.df<-as.data.frame(upland.raster.distance,xy=TRUE)
upland.raster.distance.df$ID<-row.names(upland.raster.distance.df)
names(upland.raster.distance.df)[names(upland.raster.distance.df)=="x"]<-"Longitude"
names(upland.raster.distance.df)[names(upland.raster.distance.df)=="y"]<-"Latitude"
names(upland.raster.distance.df)[names(upland.raster.distance.df)=="layer"]<-"Dist_to_upland_m"

# write.csv(upland.raster.distance.df, file="./Data/Predictive_surface/Dist_to_upland_pred_surface.csv", row.names=FALSE)

write.csv(upland.raster.distance.df, file="./Data/Predictive_surface/Dist_to_upland_pred_surface_hi_res.csv", row.names=FALSE)

#######################################################################################################
#Dist to ocean
oceanLine.sp<-spLines(cbind(oceanLine$x, oceanLine$y))
#class(oceanLine.sp)
#plot(oceanLine.sp)

#ocean.raster.distance = raster(extent(new.raster.blank),217,262)
#hi-res version, need to get new coords here: 217*262/90000*0.5102041
#ocean.raster.distance = raster(extent(new.raster.blank),427,516)

#dim(new.raster)
#modify dimension to match new raster
ocean.raster.distance = raster(extent(new.raster.blank),577,752)


#dim(ocean.raster.distance)
ocean.raster.distance[]<-1
p = as(ocean.raster.distance,"SpatialPoints")
d = gDistance(p, oceanLine.sp, byid=TRUE)
dim(d)
dmin = apply(d,2,min)
ocean.raster.distance[]=dmin
plot(ocean.raster.distance, col=cubicyf(100),main="Distance to Ocean (m)")
plot(oceanLine.sp,lwd=2, col='red',add=TRUE)

#convert raster to data.frame and rename x and y
ocean.raster.distance.df<-as.data.frame(ocean.raster.distance,xy=TRUE)
ocean.raster.distance.df$ID<-row.names(ocean.raster.distance.df)
names(ocean.raster.distance.df)[names(ocean.raster.distance.df)=="x"]<-"Longitude"
names(ocean.raster.distance.df)[names(ocean.raster.distance.df)=="y"]<-"Latitude"
names(ocean.raster.distance.df)[names(ocean.raster.distance.df)=="layer"]<-"Dist_to_ocean_m"

# write.csv(ocean.raster.distance.df, file="./Data/Predictive_surface/Dist_to_ocean_pred_surface.csv", row.names=FALSE)

##################################################################################
# move ocean line to capture coastal points that were outside of DE bay

#PatchRemoveList<-c("9121","9105", "8938","8878", "8572","8511","8423","8392","8259", "8162","7585","7257","7239","7221","7220","7186","7118","6181","6153")


###################################################################################

write.csv(ocean.raster.distance.df, file="./Data/Predictive_surface/Dist_to_ocean_pred_surface_hi_res.csv", row.names=FALSE)

Combine all covariates together

#read in pred.veg.proportions and add other covariates
# pred.veg.proportions.df.cast<-read.csv("./Data/Predictive_surface/Vegetation_pred_surface.csv")
pred.veg.proportions.df.cast<-read.csv("./Data/Predictive_surface/Vegetation_pred_surface_hi_res.csv")

pred.veg.proportions.df.cast$ID<-as.factor(pred.veg.proportions.df.cast$ID)

#Other covariates
# patch.area<-read.csv("./Data/Predictive_surface/Patch_area_pred_surface.csv")
patch.area<-read.csv("./Data/Predictive_surface/Patch_area_pred_surface_hi_res.csv")

names(patch.area)[names(patch.area)=="Patch_area_ha"]<-"Patch_Area"
patch.area$ID<-as.factor(patch.area$ID)
patch.area$Longitude<-NULL
patch.area$Latitude<-NULL

# tidal.res<-read.csv("./Data/Predictive_surface/Tidal_restriction_density_pred_surface.csv")
tidal.res<-read.csv("./Data/Predictive_surface/Tidal_restriction_density_pred_surface_hi_res.csv")

names(tidal.res)[names(tidal.res)=="Tidal_restriction_density"]<-"Tidal_Restriction"
tidal.res$ID<-as.factor(tidal.res$ID)
tidal.res$Longitude<-NULL
tidal.res$Latitude<-NULL

# marsh.ditches<-read.csv("./Data/Predictive_surface/Marsh_ditch_density_pred_surface.csv")
marsh.ditches<-read.csv("./Data/Predictive_surface/Marsh_ditch_density_pred_surface_hi_res.csv")

names(marsh.ditches)[names(marsh.ditches)=="Ditch_density"]<-"Marsh_Ditch"
marsh.ditches$ID<-as.factor(marsh.ditches$ID)
marsh.ditches$Longitude<-NULL
marsh.ditches$Latitude<-NULL

# dist_water<-read.csv("./Data/Predictive_surface/Dist_to_water_pred_surface.csv")
dist_water<-read.csv("./Data/Predictive_surface/Dist_to_water_pred_surface_hi_res.csv")

names(dist_water)[names(dist_water)=="Dist_to_water_m"]<-"Dist_to_water"
dist_water$ID<-as.factor(dist_water$ID)
dist_water$Longitude<-NULL
dist_water$Latitude<-NULL

# dist_upland<-read.csv("./Data/Predictive_surface/Dist_to_upland_pred_surface.csv")
dist_upland<-read.csv("./Data/Predictive_surface/Dist_to_upland_pred_surface_hi_res.csv")

names(dist_upland)[names(dist_upland)=="Dist_to_upland_m"]<-"Dist_to_upland"
dist_upland$ID<-as.factor(dist_upland$ID)
dist_upland$Longitude<-NULL
dist_upland$Latitude<-NULL

# dist_ocean<-read.csv("./Data/Predictive_surface/Dist_to_ocean_pred_surface.csv")
dist_ocean<-read.csv("./Data/Predictive_surface/Dist_to_ocean_pred_surface_hi_res.csv")

names(dist_ocean)[names(dist_ocean)=="Dist_to_ocean_m"]<-"Dist_to_ocean"
dist_ocean$ID<-as.factor(dist_ocean$ID)
dist_ocean$Longitude<-NULL
dist_ocean$Latitude<-NULL

#merge multiple data.frames
all.covs.merge<-Reduce(function(x, y) merge(x, y, by=c("ID"), all.x=TRUE), 
                       list(pred.veg.proportions.df.cast, patch.area, tidal.res, marsh.ditches, dist_water,
                            dist_ocean, dist_upland))

#save predictive surface as .csv
# write.csv(all.covs.merge, file="./Data/Predictive_surface/ALL_covs_pred_surface.csv", row.names=FALSE)

write.csv(all.covs.merge, file="./Data/Predictive_surface/ALL_covs_pred_surface_hi_res.csv", row.names=FALSE)

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

#read in new data with Veg layer covariates
callback2020_DE_NJ_covs<-read.csv("./Data/Callback_2020_DE_NJ_with_ALL_Covs.csv")

#create unmarkedFrame
new.umf<-makePOccu_umf_with_covs(DataIn=callback2020_DE_NJ_covs, SpeciesIn="SALS")

#get data.frame of just the covariates
covariate.df<-as.data.frame(new.umf@siteCovs)[,-1:-3]

#reduce to include only variables we model
covariate.df<-covariate.df[,c("High_Marsh","Low_Marsh","Mudflat","Phragmites",
                              "Pool_Pan","Stream","Terrestrial_Border","Upland","Marsh_Ditch",
                              "Tidal_Restriction")]

#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 = "./Results/Figures/Correlation_matrix_covs.jpeg", width=10, height=10, units="in",res=600)

print(corrPlot)

dev.off()

## quartz_off_screen 
##                 2

#show plot by loading saved image
include_graphics("./Results/Figures/Correlation_matrix_covs.jpeg")

Figure 5. Explore covariate relationships to detection probabilty

#read in new data with Veg layer covariates
callback2020_DE_NJ_covs<-read.csv("./Data/Callback_2020_DE_NJ_with_ALL_Covs.csv")

#get speciesList
speciesList<-c("AMOY","GRYE", "LEYE", "RUTU","SESA")

#make list of detection process covariates (non-factors only)
detformulaList<-c(~Visit, ~Day,~Month, ~Time, ~Temp, ~Wind, ~Sky, ~Noise,
  ~High_Marsh, ~Low_Marsh, ~Mudflat, ~Phragmites, ~Pool_Pan, ~Stream, 
  ~Terrestrial_Border, ~Upland)

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

#for loop to iterate over species and estimate probability of occupancy
for(i in 1:length(speciesList)){
  
  new.AlphaCode<-speciesList[i]
  #print(new.AlphaCode)

  #build new unmarked frame
  new.umf<-makePOccu_umf_with_covs(DataIn=callback2020_DE_NJ_covs, SpeciesIn=new.AlphaCode)
  
  new.umf.test<-new.umf@y

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

  #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"
    
    #add column with covariate name
    pred.df$CovariateName<-varName
    
    #add Species AlphaCode
    pred.df$AlphaCode<-new.AlphaCode
    
    DetCovariatePred<-rbind(DetCovariatePred, pred.df)
    
    }
  
  #compile covariates for each species
occu.mod.save<-rbind(occu.mod.save, DetCovariatePred)
}

#make CovariateName a factor
occu.mod.save$CovariateName<-factor(occu.mod.save$CovariateName, levels=c("Month","Day","Visit","Time","Temp","Sky","Wind","Noise","High_Marsh","Low_Marsh", "Mudflat","Phragmites","Pool_Pan","Stream","Terrestrial_Border","Upland"))

#now plot 

  #plot p(Occupancy) vs. Vegetation class proportions
  det.prob.plot<-ggplot(data=occu.mod.save, aes(x=det_cov, y=Predicted, group=AlphaCode))+
    geom_point(aes(color=CovariateName),alpha=0.1,size=1)+
    geom_smooth(method = "glm", formula = y~x, method.args = list(family = binomial(link = 'log')), 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,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~AlphaCode,scales="free",ncol=5)
  #det.prob.plot.facet<-det.prob.plot+facet_grid(CovariateName~AlphaCode, scales="free")
  print(det.prob.plot.facet)

  #create filename with timestamp
  new.filename<-paste("Detection_Probability_covariates_",Sys.time(),".png",sep="")
  
  ggsave(det.prob.plot.facet, file=paste("./Results/Figures/",new.filename, sep=""),width=10, height=25, dpi=600)

Fit occupancy models with site-level covariates.

#read in new data with Veg layer covariates
callback2020_DE_NJ_covs<-read.csv("./Data/Callback_2020_DE_NJ_with_ALL_Covs.csv")

#iterate over species and fit occupancy models with veg layers

#get speciesList
speciesList<-c("AMOY","GRYE", "LEYE", "RUTU","SESA")

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

#for loop to iterate over species and estimate probability of occupancy
for(i in 1:length(speciesList)){
  
  new.AlphaCode<-speciesList[i]
  #print(new.AlphaCode)

  #build new unmarked frame
  new.umf<-makePOccu_umf_with_covs(DataIn=callback2020_DE_NJ_covs, SpeciesIn=new.AlphaCode)
  
  #new.umf.test
  new.umf.test<-new.umf@y
  
  #######TESTING
  # #Check out Patch_Area variable with plot
  # patch.area.plot<-ggplot(data=new.umf@siteCovs, aes(x=Longitude, y=Latitude))+
  #   geom_point(aes(color=new.umf@siteCovs$Dist_to_upland))
  # 
  # patch.area.plot
  
  ##############

  #fit model with Veg covariates
  mod.global<-NULL
  mod.global<-occu(~Visit ~High_Marsh+Low_Marsh+Mudflat+Phragmites+Pool_Pan+Stream+Terrestrial_Border+Upland+Marsh_Ditch+Tidal_Restriction, data=new.umf)
    # mod.global<-occu(~Visit ~Dist_to_water, data=new.umf)
    # mod.global<-occu(~Visit ~Dist_to_upland, data=new.umf)
    # mod.global<-occu(~Visit ~Marsh_Ditch, data=new.umf)
    # mod.global<-occu(~Visit ~Tidal_Restriction, data=new.umf)
    # mod.global<-occu(~Visit ~Terrestrial_Border, data=new.umf)
    #Like Dist_to_ocean, also leads to 0.5 across the board. . .need to look into this
    # mod.global<-occu(~Visit ~Patch_Area, data=new.umf) 

  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 GOF_cHat
  mod.predict$GOF_cHat<-mod.c.hat
  
  
  #add Alphacode
  mod.predict$AlphaCode<-new.AlphaCode
  
  #simplify
  mod.predict.sub<-unique(mod.predict[,c("AlphaCode","GOF_cHat","PointID","Predicted","SE","lower","upper","High_Marsh","Low_Marsh","Mudflat","Phragmites","Pool_Pan","Stream","Terrestrial_Border","Upland","Patch_Area","Dist_to_ocean","Marsh_Ditch","Tidal_Restriction")])
  
  #save results
  occu.mod.save<-rbind(occu.mod.save, mod.predict.sub) 
}

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

  #melt
  mod.predict.melt<-melt(occu.mod.save, id.vars=c("AlphaCode","PointID","Predicted","SE","lower","upper"), measure.vars=c("High_Marsh","Low_Marsh","Mudflat","Phragmites","Pool_Pan","Stream","Terrestrial_Border","Upland", "Marsh_Ditch","Tidal_Restriction"), variable.name = "Variable", value.name="value")
  
    #make AlphaCode a factor
  mod.predict.melt$AlphaCode<-factor(mod.predict.melt$AlphaCode, levels=speciesList)

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

  #plot p(Occupancy) vs. Vegetation class proportions
  mod.occu.plot<-ggplot(data=mod.predict.melt, aes(x=value, y=Predicted, group=AlphaCode))+
    geom_point(aes(color=Variable),alpha=0.1,size=1)+
    geom_smooth(method = "glm", formula = y~x, method.args = list(family = binomial(link = 'log')), aes(color=Variable, fill=Variable),alpha=0.2,fullrange=FALSE)+
    #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~AlphaCode,scales="free",ncol=5)
  #mod.occu.plot.facet<-mod.occu.plot+facet_grid(Variable~AlphaCode,scales="free")
  mod.occu.plot.facet

  #create filename with timestamp
  new.filename<-paste("Fig_3_Predicted_Occupancy_models_with_veg_class_",Sys.time(),".png",sep="")
  
  ggsave(mod.occu.plot.facet, file=paste("./Results/Figures/",new.filename, sep=""),width=10, height=16, dpi=600)

Fit global model and explore variable importance

#read in new data with Veg layer covariates
callback2020_DE_NJ_covs<-read.csv("./Data/Callback_2020_DE_NJ_with_ALL_Covs.csv")

#iterate over species and fit occupancy models with veg layers

#get speciesList
speciesList<-c("AMOY","GRYE", "LEYE", "RUTU","SESA")

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

#empty list to gather coefs results
coefs.occu.save<-list()

#for loop to iterate over species and estimate probability of occupancy
for(i in 1:length(speciesList)){
  
  suppressMessages({

  new.AlphaCode<-speciesList[i]
  #print(new.AlphaCode)

  #build new unmarked frame
  new.umf<-makePOccu_umf_with_covs(DataIn=callback2020_DE_NJ_covs, SpeciesIn=new.AlphaCode)

  #fit model with Veg covariates
  mod.global<-NULL
  mod.global<-occu(~Visit ~High_Marsh+Low_Marsh+Mudflat+Phragmites+Pool_Pan+Stream+Upland+Marsh_Ditch+Tidal_Restriction, data=new.umf)
  
  #get model summary
  mod.summary<-summary(mod.global)
  
  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
  
  #add AlphaCode
  new.coefs$AlphaCode<-new.AlphaCode
  
  coefs.occu.save<-rbind(coefs.occu.save, new.coefs)
  
  # # get list of all possible terms, then subset to those we want to keep
  # det_terms <- getAllTerms(mod.global) %>% 
  # # retain the detection submodel covariates
  # discard(stringr::str_detect, pattern = "psi")
  # 
  # 
  # model_terms <- getAllTerms(mod.global)

  # # add to prediction surface
  # pred_occ <- bind_cols(pred_surface, occ_prob = occ_pred$fit,occ_se = occ_pred$se.fit) %>% 
  # select(Latitude, Longitude, occ_prob, occ_se)
  # 
  #  r_pred <- pred_occ %>% 
  #   # convert to spatial features
  #   sf::st_as_sf(coords = c("Longitude", "Latitude"), crs = crs(new.raster)) %>% 
  #   sf::st_transform(crs = raster::projection(r)) %>% 
  #   # rasterize
  #   rasterize(r)
  # r_pred <- r_pred[[c("occ_prob", "occ_se")]]
  # 
  # #save raster
  # writeRaster(r_pred[["occ_prob"]], 
  #           filename = file.path("./Results/Figures/", "occupancy-model_prob.tif"),
  #           overwrite = 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 GOF_cHat
  mod.predict$GOF_cHat<-mod.c.hat
  
  #add Alphacode
  mod.predict$AlphaCode<-new.AlphaCode
  
  #simplify
  mod.predict.sub<-unique(mod.predict[,c("AlphaCode","GOF_cHat","PointID","Predicted","SE","lower","upper","High_Marsh","Low_Marsh","Mudflat","Phragmites","Pool_Pan","Stream","Terrestrial_Border","Upland","Patch_Area","Dist_to_ocean","Marsh_Ditch","Tidal_Restriction")])
  
  #save results
  occu.mod.save<-rbind(occu.mod.save, mod.predict.sub) 
  
  })

}

Figure 7. Plot of occupancy model-estimated beta coefficients for site-level covariates among species. 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. Note differing scales on x-axis.

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

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

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

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

#set up factor levels for species and variables
coefs.occu.save.sub$AlphaCode<-factor(coefs.occu.save.sub$AlphaCode, levels=c("AMOY","GRYE", "LEYE", "RUTU","SESA"))

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

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

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

  #plot beta coefficients
  betaCoefPlot_occu<-ggplot(data=coefs.occu.save.sub.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)+
        coord_flip()+
        scale_color_manual(values=myBetaColors)+
        #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+facet_wrap(.~AlphaCode, scales="free_x",nrow = 1)

Use occupancy model with covs to get predicted estimates

#load predictive surface
pred_surface<-read.csv("./Data/Predictive_surface/ALL_covs_pred_surface_hi_res.csv")

#get speciesList
speciesList<-c("AMOY","GRYE", "LEYE","RUTU","SESA")

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

library(MuMIn)

# dredge.mod <- dredge(global.mod)
# dredge.mod
# dredge.mod.top<-subset(dredge.mod, delta < 4)
# 
# # Visualize the model selection table:
# 
# par(mar = c(3,5,6,4))
# plot(dd, labAsExpr = TRUE)
# 
# 
# # Model average models with delta AICc < 4
# model.avg(dd, subset = delta < 4)
# 
# #or as a 95% confidence set:
# model.avg(dd, subset = cumsum(weight) <= .95) # get averaged coefficients
# 
# #'Best' model
# summary(get.models(dd, 1)[[1]])



#for loop to iterate over species and estimate probability of occupancy
for(i in 1:length(speciesList)){
  
    suppressMessages({

  new.AlphaCode<-speciesList[i]
  #print(new.AlphaCode)

  #build new unmarked frame
  new.umf<-NULL
  new.umf<-makePOccu_umf_with_covs(DataIn=callback2020_DE_NJ_covs, SpeciesIn=new.AlphaCode)

  #fit model with Veg covariates
  # new.mod<-occu(~Visit ~High_Marsh+Low_Marsh+Mudflat+Phragmites+Pool_Pan+Stream+Upland+Terrestrial_Border+Patch_Area+Dist_to_ocean+Marsh_Ditch+Tidal_Restriction, data=new.umf)
  new.mod<-NULL
  new.mod<-occu(~Visit ~High_Marsh+Low_Marsh+Mudflat+Phragmites+Pool_Pan+Stream+Upland+Marsh_Ditch+Tidal_Restriction, data=new.umf)
  
  # get predicted estimates of avian occupancy as a function of vegetation cover on map
  occ_pred <- predict(new.mod, 
                    newdata = as.data.frame(pred_surface), 
                    type = "state",appendData=TRUE)

# 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)

#add species name
pred_occ$AlphaCode<-new.AlphaCode

occu.pred.save<-rbind(occu.pred.save, pred_occ)

  })

}

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

#make AlphaCode a factor
pred_occ.melt$AlphaCode<-factor(pred_occ.melt$AlphaCode,levels=speciesList)

Figures 8A-F. Predicted probability of species occupancy (mean ± SE) around DE Bay.

letterList<-c("A","B","C","D","E","F")

for(i in 1:length(speciesList)){
  
  pred_occ.melt.sub<-subset(pred_occ.melt, AlphaCode==speciesList[i])

  map.occu.plot<-ggmap(basemap,darken = c(0.4, "black"))+
    geom_tile(data=subset(pred_occ.melt.sub, Measure=="occ_prob"), aes(x=Longitude, y=Latitude, fill=`p(Occupancy)`,color=`p(Occupancy)`))+
    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(paste("Figure 8",letterList[i],". ",sep=""),speciesList[i]," - Predicted probability of occupancy around Delaware Bay",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=15, st.dist=0.03, height = 0.015, anchor=c(x=-74.5,y=38.55), border.size=0.5,st.size=2)

  map.occu.plot.facet<-map.occu.plot+facet_grid(~Measure)
  
  map.occu.plot.facet.save<-map.occu.plot+facet_grid(~Measure)
  
  #create filename with timestamp
  new.filename<-paste(paste(speciesList[i],"_Map_Predicted_Occupancy_DE_bay_",sep=""),Sys.time(),".png",sep="")
    
  ggsave(map.occu.plot.facet.save, file=paste("./Results/Figures/",new.filename, sep=""),width=14, height=8, dpi=600)
  
print(map.occu.plot)
  
}

Plot and save NOAA 0.3-m sea level rise map.

#read in shape file of SLR DE Bay 1ft
slr1ft<-readOGR(dsn="./Data/SLR_data/", layer="DE_Bay_slr_1ft")

#str(slr1ft)

# #plot and save
# png(filename="./Results/Figures/DE_Bay_1ft_slr_vector.png", width=10, height=10)
# 
# plot(slr1ft, col= "seagreen3")
# 
# dev.off()

#fortify
slr1ft_df <- fortify(slr1ft)
#names(slr1ft_df)

#reorder
slr1ft_df_sort<-slr1ft_df[order(slr1ft_df$order, decreasing=FALSE),]

#make factors
slr1ft_df_sort$id<-factor(slr1ft_df_sort$id, levels=unique(slr1ft_df_sort$id))
slr1ft_df_sort$group<-factor(slr1ft_df_sort$group, levels=unique(slr1ft_df_sort$group))

slrPlot<-ggmap(basemap,darken = c(0.4, "black"))+
  geom_polygon(data=slr1ft_df_sort, aes(x=long,  y=lat, group=1, fill="red"),alpha=0.6,inherit.aes = FALSE)+
  theme(panel.background=element_rect(fill="white",color="black"),
        panel.border = element_rect(fill="transparent",color="black"))

#save plot
ggsave(slrPlot, file="./Results/Figures/DE_Bay_1ft_slr.png", width=10, height=10)

Figure 9. NOAA 0.3-m (1 foot) sea level rise data to estimate how much coastal habitat is lost for breeding birds.

#Load the saved image of raster veg layer
#include_graphics("./Results/Figures/DE_Bay_1ft_slr.png")

include_graphics("./Results/Figures/SLR_1ft_NOAA.png")

Figures 10A-F. Predicted probability of species occupancy (mean ± SE) around DE Bay with 60% reduction in habitat due to 0.3-m (1 foot) sea level rise.

slr1ft<-readOGR(dsn="./Data/SLR_data/", layer="DE_Bay_slr_1ft")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/zach/ZachTeam Dropbox/Zach Ladin/Projects/USFWS_DE_Bay_Shorebirds/Rmarkdown/Data/SLR_data", layer: "DE_Bay_slr_1ft"
## with 53617 features
## It has 7 fields
## Integer64 fields read as strings:  Id gridcode

for(i in 1:length(speciesList)){
  
  pred_occ.melt.sub<-subset(pred_occ.melt, c(AlphaCode==speciesList[i] & Measure=="occ_prob"))
  
  xyz.df<-pred_occ.melt.sub[,c("Longitude","Latitude","p(Occupancy)")]
  
  #convert to raster
  pred_occ.melt.sub.ras<-rasterFromXYZ(xyz.df)
  
  #save
  writeRaster(pred_occ.melt.sub.ras, filename = paste("./Results/Predicted_rasters",paste(speciesList[i],"_prob_occupancy_1ft.tif",sep=""),sep=""), options=c('TFW=YES'),overwrite=TRUE)

  #how many cells with values
  
  #now mask predicted map with slr 1 ft vector
  pred_occ.melt.sub.ras.mask<-mask(pred_occ.melt.sub.ras,slr1ft, inverse=TRUE)
  
  writeRaster(pred_occ.melt.sub.ras.mask, filename = paste("./Results/Predicted_rasters",paste(speciesList[i],"_prob_occupancy_1ft_SLR.tif",sep=""),sep=""), options=c('TFW=YES'), overwrite=TRUE)

  #plot(pred_occ.melt.sub.ras.mask, col=parula(100))
  
  #convert to data.frame
  pred_occ.melt.sub.ras.mask.df<-na.omit(as.data.frame(pred_occ.melt.sub.ras.mask,xy=TRUE))
  names(pred_occ.melt.sub.ras.mask.df)<-c("Longitude","Latitude","p(Occupancy)")
  
  #compute percent change in number of cells
  nCellsFull<-nrow(pred_occ.melt.sub)
  nCellsMask<-nrow(pred_occ.melt.sub.ras.mask.df)
  pctChange<- (nCellsFull - nCellsMask)/abs(nCellsFull)*100
  
  map.occu.plot.slr<-ggmap(basemap,darken = c(0.4, "black"))+
    geom_tile(data=pred_occ.melt.sub.ras.mask.df, aes(x=Longitude, y=Latitude, fill=`p(Occupancy)`,color=`p(Occupancy)`))+
    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(paste("Figure 13",letterList[i],". ",sep=""),speciesList[i]," - Predicted probability of occupancy around Delaware Bay with 0.3-m sea level rise",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=15, st.dist=0.03, height = 0.015, anchor=c(x=-74.5,y=38.55), border.size=0.5,st.size=2)
  
print(map.occu.plot.slr)

}

References

Correll, M.D., Hantson, W., Hodgman, T.P., Cline, B.B., Elphick, C.S., Gregory Shriver, W., Tymkiw, E.L. and Olsen, B.J., 2019. Fine-scale mapping of coastal plant communities in the northeastern USA. Wetlands, 39, pp.17-28.

Fiske, I. and Chandler, R., 2011. Unmarked: an R package for fitting hierarchical models of wildlife occurrence and abundance. Journal of statistical software, 43, pp.1-23.

Ladin, Z.S., Wiest, W.A., Correll, M.D., Tymkiw, E.L., Conway, M., Olsen, B.J., Elphick, C.S., Thompson, W.L. and Shriver, W.G., 2020. Detection of local-scale population declines through optimized tidal marsh bird monitoring design. Global Ecology and Conservation, 23, p.e01128.

Office for Coastal Management, 2020: NOAA Office for Coastal Management Sea Level Rise Data: 1-10 ft Sea Level Rise Inundation Extent, https://www.fisheries.noaa.gov/inport/item/48106.

SHARP 2017. “Bird survey database: 2011-2014”. Saltmarsh Habitat and Avian Research Program. https://www.tidalmarshbirds.org.

Wiest, WA, MD Correll, BJ Olsen, CS Elphick, TP Hodgman, DR Curson, WG Shriver. 2016. Population estimates for tidal marsh birds of high conservation concern in the northeastern USA from a design-based survey. Condor: Ornithological Applications 118(2):274-288. https://doi.org/10.1650/CONDOR-15-30.1