#' ## Secretive marshbird response to invasive wetland plant management #' #' ### Nina Hill, hillx725@umn.edu #' David E. Andersen, dea@umn.edu #' #' University of Minnesota, Minnesota Cooperative Fish & Wildlife Research Unit #' #' 135 Skok Hall, Upper Buford Circle, Saint Paul, Minnesota 55108 #' format(Sys.Date(), "%Y %B %d") # Code related to Chapter 1 - Estimating avian densities in heterogenous habitats # ## Prepare marshbird and wetland data, calculate habitat-informed detection probabilities # ## Apply habitat-informed detection probabilities to maximum counts per year to estimate density of birds per point # Set up workspace, echo=FALSE------------------------- setwd("F:/forDRUM") library(plyr); library(dplyr); library(tidyr);library(ggplot2) rm(list=ls(all=TRUE)) # Clear out any objects in memory right = function(text, num_char) { substr(text, nchar(text) - (num_char-1), nchar(text)) } # Part A: Format data and calculate habitat-informed detection probabilities---------------------------- # 1) Import data MBIRD and MSURV, combine into MDATA and reshape into MPOINT-------- # create list of Buffer distances for binning detection counts buffdist <- data.frame(BuffLow=seq(0,600,by=23)) buffdist$Mid <- buffdist$BuffLow + 11.5 buffdist$Max <- buffdist$BuffLow + 23 BuffLow=buffdist$BuffLow BuffHigh=buffdist$Max[1:26] # Make a list of the species of interest bird.broadcast <- c("AMBI","LEBI","PBGR","SORA","VIRA") # bird.secondary=c("YHBL", "WISN", "AMCO", "BLTE", "UPSA", "NO BIRDS DETECTED") birdcolors <- c("orange", "turquoise","purple","tomato","dodgerblue") # Load BIRD detection data----------------------------------------------- # for ALL detections (see next section for subsetting to only during-survey detections). # I did a lot of data wrangling already in Excel, including creating many columns # of summarized data, and creating extra identifiers (ex. OutsideTimeONLY) # 1783 obs in full dataset # 1514 obs reduced by removing repeated and secondary species in Excel # NOTE 6/15/17: # I imported Morris_data and Morris_DistanceBins into ArcMap to determine IN/OUT # I updated the IN/OUT data for MSURV and added it to this script so it can select MBIRD.SURV data MBIRD <- read.csv("data/01_MBIRD.csv") %>% filter(Species %in% bird.broadcast) %>% select("SurveyID", "DetectionRecID", "PointID", "Species", "Distance", "OutsideTimeONLY", "OutsideTarget", "Previous", "IN_NWI16", "IN_Treat") %>% # OutsideTarget column was recorded by field observer mutate(BinLow = cut(Distance, breaks=BuffLow, include.lowest=TRUE, right=FALSE)) %>% # cut the distance values into 23-m distance bins mutate(Count = 1) levels(MBIRD$BinLow) <- BuffLow # Load SURVEY details--------------------------------------------------- MSURV.orig <- read.csv("data/04_MSURV.csv") MSURV <- MSURV.orig %>% select("SurveyID", "Year", "Round", "PointID", "Treatment", "WetlandID", "OtherWet", "OutTarget", "WetCONF") %>% # RouteName mutate(Treatment = factor(Treatment, levels =c("Low","Moderate","High","Extra High"))) %>% mutate(Year = factor(Year)) %>% mutate(SurveyID = as.integer(SurveyID)) %>% mutate(YRround = paste(Year, Round, sep="_")) # RouteList <- levels(MSURV$RouteName) # MSURV$RouteName <- factor(MSURV$RouteName, RouteList) MSURVvisit <- MSURV %>% dplyr::count(PointID, Round, Year , name = "visits", .drop=FALSE) # %>% # pivot_wider(PointID, names_from = c(Year, Round), # values_from = Year, # values_fn = list(Year = length), # values_fill=0, # names_sort = TRUE) # number of survey visits each round-year aggregate(list(nvisits=MSURVvisit$visits), list(Round=MSURVvisit$Round, Year=MSURVvisit$Year), sum) # number of surveys we did each year-round novisit15 <- MSURVvisit %>% filter(Year == 2015 & visits == 0) %>% dplyr::count(PointID, Year, name="novisit15") %>% filter(novisit15 == 2) novisit16 <- MSURVvisit %>% filter(Year == 2016 & visits == 0) %>% dplyr::count(PointID, Year, name="novisit16") %>% filter(novisit16 == 2) # Notes--- # MSURV = each record represents a survey conducted during 2015 & 2016. # Each survey location(PointID) has up to 4 survey records. # Some have fewer: # a) 2 or 3 points dropped after 2015 R1, # b) 12 points added 2015 R2, # c) 16 points added new 2016 R1, # note: no new points added 2016 R2. # 467 unique surveys conducted (in Excel, confirmed here) # 348 survey with broadcast species detected (119 surveys conducted when no broadcast species detected, about 25% (0.254818)) # 42 surveys when only secondary, but NO broadcast spp detected # 75 surveys when NO BIRDS DETECTED = No broadcast& No secondary=NONE # Sum ALL detection counts of broadcast species into MDATA.ALL---------------- # Leave all the detections in the datset incl OUTside target wetland, # or OUTside survey time for this analysis of distance/density calculations birdcounts <- MBIRD %>% filter(Species %in% bird.broadcast) %>% pivot_wider(id_cols = c(SurveyID, PointID), names_from = Species, values_from = DetectionRecID, values_fn = list(DetectionRecID = length), values_fill = 0) MSURV <- left_join(MSURV, birdcounts, by=c("SurveyID", "PointID"), all=T) # merge MSURV survey and MBIRD bird detection data, long format MDATA.ALL <- left_join(MSURV, MBIRD, by=c("SurveyID", "PointID")) %>% droplevels() # total counts going into habitat-informed detection probability analysis print(table(MDATA.ALL$Species, MDATA.ALL$Year)) # quick summary of species counts # some detections are missing distance estimates from field. I estimated several using notes and imagery in GIS # the rest I could fil in with median value -- but I have chosen to leave out print(table(MDATA.ALL$Species[is.na(MDATA.ALL$Distance)])) # 185 AMBI detections, but 177 detections with distance estimations # 8 obs w/out distance estimates @ 802E, 24E, 77H, 2028H, 90H, # MDATA.ALL <- MDATA.ALL %>% # group_by(Species) %>% # mutate(Distance = replace(Distance, is.na(Distance), median(Distance, na.rm = T))) rm(MBIRD, MSURV, birdcounts) # 2) Prepare Wetland data----------------------------------------- # Import and format GIS ALL (landscape) wetland data # MN NWI2016 and manually updated wetland polygons clipped into 23m-rings around survey point # 6/8/17 I made a single shapefile with wetland polys for ALL points. I ran all the tools through # model builder in GIS and exported the .dbf. Use package library(foreign) to use read.dbf into R. # Load NWI wetland summaries by distance bins (all wetland cover types) wet_by_dist_ALL <- read.csv("data/wet_by_dist_ALL.csv", header=T) # LOAD PBGR WETLAND (deep and open water wetland cover types) wet_by_dist_ALL_deep <- read.csv("data/wet_by_dist_ALL_deep.csv", header=T) colnames(wet_by_dist_ALL_deep) <- c("BuffOuter", "WetHA") # Load new NWI calcs at Mpoint level--- # MPOINT.ALL_nwi <- read.csv("F:/Morris/data/MPOINT.ALL_NWI.csv", header=T) # MPOINT.ALL_deep <- read.csv("F:/Morris/data/MPOINT.ALL_deep.csv", header=T) # 3) Calculate habitat-informed detection probabilities---------------------------------------------- # MORRIS STUDY AREA - Adjustment values # Use ALL bird detections to calculate adjustment values for each species, each point, # each 23-m distance bin: # adj value for dist bin = (sum down ALL raw CTS in distance bin / sum down ALL wetHA in distance bin). # # For each species, scale the adjustment values by the first bin, so that g(0)=1. # **this is a problem for the bittern species # **SEE 03b_Morris_AMBI_LowDet0.R for more details; I copied the important code into this script. # a) Sum detections within 23m DISTANCE BINS--- # Reformat dataframe from long to wide, tally counts by species-distance mpointall <- MDATA.ALL %>% filter(!is.na(Distance)) %>% pivot_wider(id_cols = c(PointID), values_from = Count, values_fn = sum, names_from = c(Species, BinLow), names_sort = TRUE, values_fill=0) # Plots! Species counts by distance bins (unadjusted)--- for (ii in 1:length(bird.broadcast)){ aa <- mpointall[,grep(bird.broadcast[ii], colnames(mpointall))] bb <- colSums(aa[,!colnames(aa)=="PointID"], na.rm=T) myplot<-barplot(bb[1:27], names.arg = BuffLow, col = birdcolors[ii], las=2, xlab="(Lower) Bin Distance", ylab="Count", ylim=c(0,1.1*max(bb[1:27], na.rm=T)), main=c(bird.broadcast[ii],"ALL Unadjusted Counts by Distance"),space = 0) text(x=myplot, y=bb[1:27], label = round(bb[1:27],digits = 2), pos = 3, cex = 0.8, col = "grey") plotname <- paste(bird.broadcast[ii], "ALLctsbydistance_withOUTmedian.png", sep="") dev.copy(png, file=plotname) dev.off() } rm(ii, aa, bb, plotname, myplot) bird.broadcast <- bird.broadcast[c(1,3,4,5)] # not enough LEBI detections to calc adjustments, drop birdcolors <- birdcolors[c(1,3,4,5)] # b) Divide counts within distance bins by summed wetland cover area within distance bins--- # Calculate adjustment value (count/wetHA) per distance bin for each species: # the numbers here will be the adjustment value that I apply to the count data ALL.distctadjust <- data.frame(Distance = BuffLow, wet_by_dist_ALL = wet_by_dist_ALL) for (ii in 1:length(bird.broadcast)){ thissp <- bird.broadcast[ii] thisdat <- mpointall[,grep(thissp, colnames(mpointall))] colsumz <- colSums(thisdat[,grep(paste(thissp,"0", sep="_"), colnames(thisdat)):ncol(thisdat)], na.rm=T) # sum species counts within distance bins colsumz <- c(colsumz, rep(0, nrow(ALL.distctadjust)-length(colsumz))) # and fill NA distances with 0 ALL.distctadjust[,paste(thissp, "sumct", sep="")] <- colsumz # save new column for species count # calculate and save adjustment adjvalue calc by (ct/HA wet) per distance bin ALL.distctadjust[,paste(thissp, "firstadjstep", sep="")] <- colsumz / ALL.distctadjust$wetdistALL } rm(ii, thissp, thisdat, colsumz) # add the deep/open water areas to use in PBGR calcs # the adjustment adjvalue calc by (ct/HA wet) per distance bin ALL.distctadjust$wetpbgr <- c(wet_by_dist_ALL_deep$WetHA, NA) thisdat <- mpointall[ ,grep("PBGR", colnames(mpointall))] colsumz <- colSums(thisdat[ ,grep("PBGR_0", colnames(thisdat)):ncol(thisdat)], na.rm=T) colsumz <- c(colsumz, rep(0, nrow(ALL.distctadjust)-length(colsumz))) # and fill NA distances with 0 ALL.distctadjust[, "PBGR2sumct"] <- colsumz # save the sum ct by dist bin ALL.distctadjust[, "PBGR2firstadjstep"] <- colsumz / ALL.distctadjust$wetpbgr rm(thisdat, colsumz) # ALLdistct <- MDATA.ALL%>% # select(BinLow) %>% # dplyr::full_join(MDATA.ALL, buffdist, by=c("BinLow" = "BuffLow")) # filter(!is.na(Distance)) %>% # dplyr::count(Species, BinLow, .drop =F) %>% # tidyr::complete(Species, BinLow, fill=list(n=0)) %>% # pivot_wider(id_cols = c(BinLow), names_from=Species, values_from=n, values_fill=0) # c) Scale the adjustment value by the value of the first column, g(0)=1--- # Loop through species to make data.frame of summed counts, adjust values, and scaled adjustment values for (ii in 1:length(bird.broadcast)){ thissp <- bird.broadcast[ii] thiscol <- paste(thissp,"firstadjstep",sep="") # use the value in the first bin to scale the values of other bins ALL.distctadjust[,paste(thissp, "ScalAdj", sep="")] <- ALL.distctadjust[,thiscol]/ ALL.distctadjust[1,thiscol] } rm(ii, thiscol, thissp) # update ScaledAdj for PBGR2 deep/open water # use the value in the first bin to scale the values of other bins ALL.distctadjust[,"PBGR2ScalAdj"] <- ALL.distctadjust[,"PBGR2firstadjstep"] / ALL.distctadjust[1,"PBGR2firstadjstep"] # d) Determine cut-off distance estimates--- # Determine the effective survey area, or the cut-off (truncation) distance for detecting species and calc densities. # We tested 3 methods used to find cutoff distance: # 1) distance that includes 90% of detections, I also did 95% # 2) distance = 2*mean detection distance, I also did 3*mean # 3) distance where detection probability <0.05 (calculate detection probability function and find dist where <0.05) # Resulting densities were highly correlated, so we decided to use the simplest, method 1 # Method 1: where are 90% of detections?--- # build dataframe to hold Cutoff estimates for all the species bird.broadcast <- c(bird.broadcast, "PBGR2") birdcolors <- c(birdcolors, "purple") Cutoff <- data.frame(Species=c(bird.broadcast)) for (ii in 1:length(bird.broadcast)){ dat <- ALL.distctadjust[, paste(bird.broadcast[ii], "sumct", sep="")] # DistBySpp Cutoff$CTtot[ii] <- sum(dat) # 177 detections total Cutoff$CT90obs[ii] <- 0.90*Cutoff$CTtot[ii] # so the 159.3, or rather 159th detection represents 90% Cutoff$CT90distfield[ii] <- sort(MDATA.ALL$Distance[which(MDATA.ALL$Species == bird.broadcast[ii])])[Cutoff$CT90obs[ii]] Cutoff$CT90middist[ii] <- buffdist$Mid[min(which(cumsum(dat) > Cutoff$CT90obs[ii]))] # the 159th detection has a distance of Cutoff$CT90distbin[ii] <- buffdist$Max[min(which(cumsum(dat) > Cutoff$CT90obs[ii]))] # 207, 218.5m Cutoff$CT95val[ii] <- 0.95*Cutoff$CTtot[ii] # so the 168.15, or rather 168th detection represents 95% Cutoff$CT95middist[ii] <- buffdist$Mid[min(which(cumsum(dat) > Cutoff$CT95val[ii]))] # so all detections from distance bins 1-13, or 0-276, or midpoint 287.5m } rm(ii, dat) # PBGR and PBGR2 deep/open water are same here since cutoff tallies counts and we haven't added in wetland areas yet adj <- ALL.distctadjust[,c("Distance", grep("ScalAdj", colnames(ALL.distctadjust), value=T))] adj[which(adj$Distance > Cutoff$CT90distbin[which(Cutoff$Species == "AMBI")]), grep("AMBI", colnames(adj))] <- NA adj[which(adj$Distance > Cutoff$CT90distbin[which(Cutoff$Species == "PBGR")]), grep("PBGR", colnames(adj))] <- NA adj[which(adj$Distance > Cutoff$CT90distbin[which(Cutoff$Species == "SORA")]), grep("SORA", colnames(adj))] <- NA adj[which(adj$Distance > Cutoff$CT90distbin[which(Cutoff$Species == "VIRA")]), grep("VIRA", colnames(adj))] <- NA # # check # aa <- sort(MDATA.ALL$Distance[which(MDATA.ALL$Species == "AMBI")]) # (length(aa))*0.90 # 185 AMBI observations (.9=166.5), but 177 AMBI with distance estimates (.9=159.3) 250-300 # aa[160] # =250m # pp <- sort(MDATA.ALL$Distance[which(MDATA.ALL$Species == "PBGR")]) # (length(pp))*0.90 # 314*.9=282.6; 305*.9=274.5 # pp[275] # =345m # ss <- sort(MDATA.ALL$Distance[which(MDATA.ALL$Species == "SORA")]) # (length(ss))*0.90 # 358*.9=322.2; 350*.9=315 # ss[315] # =200m # vv <- sort(MDATA.ALL$Distance[which(MDATA.ALL$Species == "VIRA")]) # (length(vv))*0.90 # 157*.9=141.3; 154*.9=138.6 # vv[139] # =100m # # # cutoff distances similar both calculations # # 159.3=250m, removing NA; 166.5=240m w median # # 274.5=345m; 282.6=325m # # 315=200m; 322.2=200m # # 138.6=100m; 141.3=100m # # ggplot(MDATA.ALL %>% filter(Species == "AMBI"), aes(x=Distance))+ # geom_histogram(binwidth = 23, fill=birdcolors[1], color="grey20")+ # geom_vline(xintercept = aa[160]) # ggplot(MDATA.ALL %>% filter(Species == "PBGR"), aes(x=Distance))+ # geom_histogram(binwidth = 23, fill=birdcolors[2], color="grey20")+ # geom_vline(xintercept = pp[275]) # # ggplot(MDATA.ALL %>% filter(Species == "SORA"), aes(x=Distance))+ # geom_histogram(binwidth = 23, fill=birdcolors[3], color="grey20")+ # geom_vline(xintercept = ss[315]) # ggplot(MDATA.ALL %>% filter(Species == "VIRA"), aes(x=Distance))+ # geom_histogram(binwidth = 23, fill=birdcolors[4], color="grey20")+ # geom_vline(xintercept = vv[139]) # # rm(aa,pp,ss,vv) rm(MSURV.orig, MSURVvisit, wet_by_dist_ALL, wet_by_dist_ALL_deep) # DONE calculating adjustment values*********************************************************** # Part B: apply habitat-informed detection probabilities to adjust counts into densities--------------------------- # 1) Import and format TARGET wetland area---- # formated GIS wetland data (NWI and manually updated wetland polygons clipped into 23m buffers around survey point) # then summed across all points into total wetland area within distance bins surrounding all survey points # LOAD Target WETLAND of all NWI and the deep/open water wetland for pbgr MPOINT.targetwet <- read.csv("data/MPOINT.Target_NWI.csv", header=T) MPOINT.targetpbgr <- read.csv("data/MPOINT.TARG_deep.csv", header=T) # 80 wetlands with deep/open water cover types (remaining target wetlands other wetland types) # sum area of wetland at survey location up to cutoff distance for each species for(ss in 1:length(bird.broadcast)){ thissp <- bird.broadcast[ss] thisdist <- Cutoff[which(Cutoff$Species == thissp), "CT90distbin"] thiscol <- paste("WetHA", thisdist, sep="_") MPOINT.targetwet[,paste(thissp, "WetHA", sep="_")] <- rowSums(MPOINT.targetwet[,c(3:grep(thiscol, colnames(MPOINT.targetwet)))], na.rm = T) } # sum wetlandHA of deep/open water for PBGR and add it to same df as other wetha MPOINT.targetpbgr$PBGR2_WetHA <- rowSums(MPOINT.targetpbgr[,c(3:16)], na.rm=T) # "WetHA_BuffOuter23":"WetHA_BuffOuter322" MPOINT.targetwet$PBGR2_WetHA <- MPOINT.targetpbgr$PBGR2_WetHA[match(MPOINT.targetwet$PointID, MPOINT.targetpbgr$PointID)] MPOINT.targetwet$PBGR2_WetHA[is.na(MPOINT.targetwet$PBGR2_WetHA)] <- 0 rm(MPOINT.targetpbgr) # 2) Sum SURVEY-only detection counts for broadcast species---------------------- # subset detections that were made DURING SURVEY, IN target wetland, and individuals NOT PREVIOUSly detected # I used GIS to verify if detections were IN or OUT of target wetland ## NOTE: I used IN NWI polygons to determine IN/OUT. Verified in or out as indicated from field survey MDATA.SURV <- MDATA.ALL %>% filter(!OutsideTimeONLY == "Y") %>% filter(!Previous == "Y") %>% filter(IN_NWI16 == "IN") %>% # n=600 # checked IN/OUT in GIS with NWI2016 update select(DetectionRecID, SurveyID, PointID, Treatment, Year, Round, YRround, Species, Count, Distance, BinLow) # filter(!OutTarget == "Y") %>% # n=192 # filter(!OutsideTarget == "Y") %>% # n=650; from field observer # filter(IN_Treat == "IN") %>% # n=783 # checked in GIS is bird in/out of treatment polygon rm(MDATA.ALL) # mpointSURV <- MDATA.SURV %>% # pivot_wider(id_cols = c(PointID, Year), # names_from = c(Species, BinLow), names_sort = TRUE, # values_from = Count, values_fn = sum) %>% # complete(PointID, Year) # Notes--- # down from 1033 ALL obs to 604 SURV obs, 57% #3.17.19 to 603 # 2104 records = 5 species for the times I did surveys (131 PointIDs * 2 years * 2 rounds) # 5 POINTS when NO BIRDS DETECTED ever during the 4 surveys (2x2015, 2x2016) # 68 POINTS where NO BROADCAST birds were ever recorded # 50 POINTS or surveys where NO Broadcast and NO secondary? # MPOINT.ALL data good to go. Each row represents a point, with total count of ALL detections for each species. # 6/15/17 Add IN/OUT data that I corrected in GIS to the MBIRD data # Used GIS to verify that detections were IN or OUT of target wetland # 3) Subset the max counts (between the 2 survey visits) for each year--------------------------- # Subset MDATA.SURV survey data for max birds detected in each year survYEAR<-c("2015","2016") MDATA.SURVmax<-data.frame(matrix(ncol=14, nrow=0)) # Loop that selects max detection of species for each year (between R1 and R2) for (hh in 1:length(bird.broadcast)){ sppdat <- MDATA.SURV[MDATA.SURV$Species== bird.broadcast[hh],] # subset bird species for (ii in 1:length(unique(sppdat$PointID))){ ptdat <- sppdat[sppdat$PointID == unique(sppdat$PointID)[ii],] # subset PointID if(is.na(ptdat$SurveyID[1])){ next } # if there aren't species det at that Point, save blank record for (jj in 1:length(survYEAR)){ yrdat <- ptdat[ptdat$Year == unique(ptdat$Year)[jj],] # subset the Year if(is.na(yrdat$SurveyID[1])){ next } # if there's 0 det that year, save blank record curround1<-yrdat[yrdat$Round == 1,] # subset Round1 curround2<-yrdat[yrdat$Round == 2,] # subset Round2 ifelse( nrow(curround1)>nrow(curround2), MDATA.SURVmax <- rbind(MDATA.SURVmax, curround1), # if R1 has more, save those records ifelse( nrow(curround1)% complete(PointID, Year) %>% pivot_wider(id_cols = c(PointID, Year), names_from = c(Species, BinLow), names_sort = TRUE, values_from = Count, values_fn = sum, values_fill = 0) mpointSURVmax[which(mpointSURVmax$PointID %in% novisit15$PointID & mpointSURVmax$Year == 2015), !colnames(mpointSURVmax) %in% c("PointID", "Year")] <- NA mpointSURVmax[which(mpointSURVmax$PointID %in% novisit16$PointID & mpointSURVmax$Year == 2016), !colnames(mpointSURVmax) %in% c("PointID", "Year")] <- NA # observations at target wetlands for analysis table(MDATA.SURVmax$Species, MDATA.SURVmax$Year) rm(MDATA.SURV, MDATA.SURVmax, novisit15, novisit16) # 5) Adjust MPOINT.SURVmax counts in distance bins, UP TO CUTOFF DISTANCE--------------- # 6) Calculate densities--- # Create MPOINT.SURV of MAX counts in distance bins UP TO CUTOFF DISTANCE # Adjust the SURVmax counts using the scaled detection function within distance bins # Calculate bird density per point by (adj bird count / total wetland area at point). # need to remove calcs past cutoff adjusted <- data.frame(PointID = mpointSURVmax$PointID, Year = mpointSURVmax$Year) densities <- data.frame(PointID = mpointSURVmax$PointID, Year = mpointSURVmax$Year)%>% left_join(MPOINT.targetwet[,c("PointID","AMBI_WetHA", "PBGR_WetHA", "SORA_WetHA", "VIRA_WetHA")], by="PointID") %>% left_join(MPOINT.targetwet[,c("PointID", "PBGR2_WetHA")], by="PointID") for(ss in 1:length(bird.broadcast)){ thissp <- bird.broadcast[ss] for(jj in 1:length(BuffLow)){ thiscol <- paste(thissp, BuffLow[jj], sep="_") kk <- mpointSURVmax[,grep(paste("^",thiscol,"$", sep=""), colnames(mpointSURVmax))] adjusted[,thiscol] <- kk / adj[, paste(thissp, "ScalAdj", sep="")][jj] } # Sum across distances for the adjusted count per survey point densities[,paste(thissp, "adjct", sep="_")] <- rowSums(adjusted[,grep(thissp, colnames(adjusted))], na.rm = T) # Divide by wetlandHA at survey point to get densities densities[,paste(thissp, "dens", sep="_")] <- densities[,paste(thissp, "adjct", sep="_")] / densities[,paste(thissp, "WetHA", sep="_")] } rm(ss, jj, kk, thissp, thiscol) # add deep/open water for PBGR thissp <- "PBGR2" for(jj in 1:length(BuffLow)){ thiscol <- paste("PBGR", BuffLow[jj], sep="_") kk <- mpointSURVmax[,grep(paste("^",thiscol,"$", sep=""), colnames(mpointSURVmax))] # get the same PBGR max counts adjusted[, paste("PBGR2", BuffLow[jj], sep="_")] <- kk / adj[,"PBGR2ScalAdj"][jj] } densities[,paste(thissp, "adjct", sep="_")] <- rowSums(adjusted[,grep(thissp, colnames(adjusted))], na.rm = T) densities[,paste(thissp, "dens", sep="_")] <- densities[,paste(thissp, "adjct", sep="_")] / densities[,paste(thissp, "WetHA", sep="_")] rm(jj, kk, thissp, thiscol) write.csv(densities, "data/MPOINT.densities.csv", row.names = F) # Plots of counts and adjustment values by distance for each species--------------------------- for (ii in 1:length(bird.broadcast)){ aa <- as.numeric(ALL.distctadjust[, grep(paste(bird.broadcast[ii], "sumct", sep=""), colnames(ALL.distctadjust))]) myplot<-barplot(aa, # names.arg = BuffLow, col = birdcolors[ii], las=2, xlab="Lower Bin Distance", ylab="Summed counts by distance bin", ylim=c(0,1.1*max(aa[!is.infinite(aa) & !is.nan(aa)])), main=c(bird.broadcast[ii],paste(""," Summed counts")), space = 0) text(x=myplot, y=aa, label = round(aa, digits = 2), pos = 3, cex = 0.8, col = "dimgrey") } rm(ii, aa, myplot, pbgr2) par(mfrow=c(1,1), mar=c(4.5,3,3,2)) # Plot detection probability adjustment values, scaled to first bin--- for (ii in 1:length(bird.broadcast)){ aa <- as.numeric(ALL.distctadjust[, grep(paste(bird.broadcast[ii], "ScalAdj", sep=""), colnames(ALL.distctadjust))]) myplot <- barplot(aa, # names.arg = BuffLow, col = birdcolors[ii], las=2, xlab="Lower Bin Distance", ylab="Scaled Adjustment Values by distance bin", ylim=c(0,1.1*max(aa[!is.infinite(aa) & !is.na(aa)])), main=c(bird.broadcast[ii],paste(""," Scaled adjustment value")), space = 0) text(x=myplot, y=aa, label = round(aa, digits = 2), pos = 3, cex = 0.8, col = "dimgrey") } rm(ii, aa, myplot) par(mfrow=c(1,1), mar=c(4.5,4,3,2)) # similar plot but xy instead of bar for (ii in 1:length(bird.broadcast)){ spdat<-as.numeric(ALL.distctadjust[,paste(bird.broadcast[ii], "ScalAdj", sep="")]) myplot <- plot(spdat, col = birdcolors[ii], las=2, pch=19, xlab="Midpoint Distance interval (m)", ylab="Scaled Adjustment Values", ylim=c(0,1.1), main=c(bird.broadcast[ii], paste(""," Scaled adjustment values by distance interval")) ) lines(spdat, col=birdcolors[ii]) text(x=c(0:30), y=spdat, label = round(spdat, digits = 2), pos = 3, cex = 0.8, col = "dimgrey") # curve() } rm(ii, spdat) # END transforming detection counts to densities******************************************************************