--- title: "Ecological Applications: Acoustic Exposure reveals spatial and temporal variation in curtailment effectiveness at reducing bat fatality at wind turbines" author: "Trevor Peterson, Adam Rusk, Caroline Byrne, and Seta Aghababian" date: "Submitted UPDATE 2024" output: html_document --- ```{r Chunk Setup, echo=FALSE} knitr::opts_chunk$set(echo = FALSE, message=FALSE, error=FALSE) ``` # Front Matter Script adapted from original analysis script by Adam Rusk (adam.rusk@westwoodps.com) May 2021 and modified for publication by Trevor Peterson (trevor.peterson@stantec.com) and Seta Aghababian (seta.aghababian@stantec.com), November 2024 ```{r Function:Install Packages} #Function to check and install a package checkPackage <- function(package_name, install=F, location=NULL, version=NULL){ #Check if package is installed if (!requireNamespace(package_name, quietly = TRUE, lib.loc=location)) { cat(paste("Package",package_name,"not installed.\n\n")) status <- FALSE #not installed #Install if specified if(isTRUE(install)){ install_version(package_name, lib=location, version=version) status <- TRUE #installed } #end install if specified } else { status <- TRUE #installed } #end check if package is installed if(isTRUE(status)){ library(package_name, lib.loc=location, character.only = TRUE) cat(paste(package_name,packageVersion(package_name,lib.loc = location), "installed and loaded.\n\n")) } #load pacakge if package is installed } ``` ```{r Function:Linear Regression Plotting} plotReg <- function(Data, x, y, xaxis, yaxis, xlims=c(0,100), ylims=c(0,100), TFont=TitleFont, KFont=TickFont, AFont=AxisFont) { #Linear model model <- lm(get(y)~get(x),Data) sum <- summary(model) print(sum) #Generate model predictions newdata <- data.frame(ins=seq(xlims[1]-50,xlims[2]+50,0.25)) newdata[[x]] <- newdata$ins newdata$ins <- NULL newdata <- cbind(newdata, predict(model,newdata=newdata,type="response",se=TRUE)) %>% mutate(LL=fit-1.96*se.fit, UL=fit+1.96*se.fit) #Pull model outputs r2 <- paste("R^2", ifelse(round(sum$r.squared)>0.001,"==","<"), #Override rounding to 1 as round to 0.99 ifelse(round(sum$r.squared,2)==1, 0.99,round(sum$r.squared,2))) p <- paste("italic('P')", ifelse(round(sum$coefficients[2,4])>0.001,"==","<"), ifelse(round(sum$coefficients[2,4],3)<0.001,0.001, round(sum$coefficients[2,4],3))) #Plot plot <- ggplot()+ coord_cartesian()+ #Error lines geom_ribbon(data=newdata,aes(x=!!sym(x),ymin = LL, ymax = UL), alpha=1,linetype=2,fill=NA,color="black",size=0.65)+ #Regression line geom_line(data=newdata,aes(x=!!sym(x),y=fit),size=1.2) + #Data points (on top of the line) geom_point(data=Data,aes(x=!!sym(x),y=!!sym(y),color=Site),size=2.5)+ #Automated x limits of +5/-5 from min/max values scale_x_continuous(limits=xlims, expand=c(0,0)) + #Automated y limits of +5/-5 from min/max values scale_y_continuous(limits=ylims, expand=c(0,0), oob=squish_infinite) + #Fix color scaling and legend title scale_color_manual(name=NULL,values=c("red","blue"), labels=c("Kings Point","North Fork Ridge")) + #Fix axis labels xlab(xaxis) + ylab(yaxis) + #R^2 annotation with scaled placement annotate("text",label=r2, x=xlims[2]*0.35, y=ylims[2]*0.85, parse=TRUE,size=5)+ #P annotation with scaled placement annotate("text",label=p, x=xlims[2]*0.35, y=ylims[2]*0.75, parse=TRUE,size=5)+ #Adjusting theme theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), plot.title = element_text(size=TFont,face="bold",hjust=0.5), axis.title.x=element_text(size=AFont, vjust = 0), axis.title.y=element_text(size=AFont, vjust=1), axis.text.x=element_text(colour = "black", size=KFont), axis.text.y=element_text(colour = "black", size=KFont), legend.title=element_text(size=AFont), legend.text=element_text(size=AFont), legend.background = element_blank(), legend.box.background=element_rect(color="black")) return(plot) } ``` ```{r Project Setup} ### Original R Studio session ran with R Version 4.3.3 ### ### Set Randomization Seed #### set.seed(333) #### Set Pathways #### wd <- getwd() #working directory Exposure.data <- list.files(wd,pattern="EchoPitch.csv", full.names=T,recursive=F) #Acoustic exposure data file Fatality.data <- list.files(wd,pattern="Carcass_Acoustics.csv", full.names=T,recursive=F) #Fatality/Exposure data file #### Libraries #### require(remotes) #Create a local library directory for package versions used by this script #Utilizing a local library with specific versions used during analysis to ensure #the code is stable for end users libraries <- paste(wd,"PetersonEtAl2025_Libraries",sep="/") #Create directory if needed if(!dir.exists(libraries)){ dir.create(libraries) cat(paste("Created directory for code-specific library versions in",libraries)) } #Libraries needed to install and load needed <- data.frame(name=c("dplyr", #for data manipulation "tidyr", #for data manipulation "tibble", #for data manipulation "ggplot2", #for plotting "lubridate", #for data manipulation "cowplot", #for plotting "mgcv", #for modeling "mgcViz",#for plotting "DHARMa", #for modeling "ggeffects", #for modeling "ggnewscale", #for plotting "scales" #for plotting ), version=c("1.1.4", "1.3.1", "3.2.1", "3.5.1", "1.9.3", "1.1.3", "1.9.1", "0.1.11", "0.4.6", "1.5.0", "0.4.10", "1.3.0")) #Install and load packages #Preference is for pacakges used during analysis to be loaded in their original #version, but if that's not possible default to already installed versions or #most recent versions published on CRAN invisible(sapply(needed$name,function(package){ tryCatch({#Check install specified version to the local package library if needed #and load library checkPackage(package,install=T,location=libraries, version=needed$version[needed$name==package])}, #If something errors out... error=function(e){ cat(paste("Error in installing/loading",package, ", loading an already installed version or", "installing/loading most recent version on CRAN.\n\n")) checkPackage(package,install=T) } ) })) #### Clean Up Setup #### remove(checkPackage,needed,libraries) ``` ```{r Data Import} #Biweekly Bat Acoustic Data by Turbine, Measured and Simulated Treatment Exposure, #and Exposure Under Additional Simulated Curtailment Strategies EchoPitch <- read.csv(Exposure.data) %>% #Fix data import formats mutate(across(c(Site,Turbine,Year),function(x){factor(x,levels=sort(unique(x)))}), DetectorPosition=factor(DetectorPosition,levels=c("Nacelle","Mid-tower")), across(Biweek:End,function(x){as.POSIXct(x,format="%Y-%m-%d",tz="Etc/GMT+6")}) ) #Bat Acoustic Exposure and Adjusted Carcass Recovery Biweekly Carcasses <- read.csv(Fatality.data) %>% #Fix data import formats mutate(across(c(Site,Year),function(x){factor(x,levels=sort(unique(x)))}), DetectorPosition=factor(DetectorPosition,levels=c("Nacelle","Mid-tower")), across(Biweek,function(x){as.POSIXct(x,format="%Y-%m-%d",tz="Etc/GMT+6")}) ) ``` # Tables ```{r Table 1} #Summary of total bat activity and exposed bat activty Table1 <- EchoPitch %>% select(-Passes_Ops) %>% #Not reporting number of passes with WS/Temp group_by(Site,Year,Treatment,DetectorPosition) %>% summarise(Start=min(Start), End=max(End), nTurb=length(unique(Turbine)), across(c(DN,Bat_Passes:Simulated),function(x){sum(x,na.rm=T)}), .groups="drop") %>% mutate(across(Bat_Passes:Simulated,list(Rate=function(x){round(x/DN,1)})), Treatment=factor(Treatment,levels=c("8.0 m/s Interim", "5.0 m/s Curtail", "Control")), Dates=paste0(format(Start,"%m/%d"),"-",format(End,"%m/%d"))) %>% select(-Start,-End) %>% relocate(Dates,.after=DetectorPosition) %>% arrange(Year,Site,Treatment) #Combine Pass and Rate Columns invisible(sapply(colnames(Table1)[11:13],function(rate.col){ pass.col <- gsub("_Rate","",rate.col) Table1[[pass.col]] <<- paste0(Table1[[pass.col]]," (",Table1[[rate.col]],")") Table1[[rate.col]] <<- NULL })) colnames(Table1) <- c("Site","Year","Operational Group","Detector Position", "Start-End Date","Number of Turbines","Detector-Nights", "Total Bat Passes (Pass Rate)", "Exposed Bat Passes (Exposed Rate)", "Simulated Exposed Bat Passes (Simulated Rate)") Table1#return table ``` ```{r Table 2} #Mean bi-weekly simulated exposure at varying cut-ins Table2 <- EchoPitch %>% pivot_longer(cols=Curt_3:Curt_8, names_to="Simulated_Treatment", values_to="Exposed_Passes") %>% mutate(Simulated_Treatment=recode(Simulated_Treatment, "Curt_3"="3.0 m/s", "Curt_4"="4.0 m/s", "Curt_4.5"="4.5 m/s", "Curt_5"="5.0 m/s", "Curt_5.5"="5.5 m/s", "Curt_6"="6.0 m/s", "Curt_6.5"="6.5 m/s", "Curt_7"="7.0 m/s", "Curt_7.5"="7.5 m/s", "Curt_8"="8.0 m/s"), Percent_Exposed=Exposed_Passes/Passes_Ops*100) %>% #Clipping to associated date range as shown in figure filter(Biweek>=as.POSIXct("2021-02-19",tz="Etc/GMT+6") & Biweek<=as.POSIXct("2021-11-26",tz="Etc/GMT+6")) %>% group_by(Simulated_Treatment,DetectorPosition) %>% summarise(Percent_Total=round(sum(Exposed_Passes)/sum(Passes_Ops)*100,1), Percent_Average=round(mean(Percent_Exposed,na.rm=T),1), .groups="drop") %>% arrange(DetectorPosition,Simulated_Treatment) Table2 ``` # Plots ```{r Establish Plot Format} #Standard Plot Format out <- getwd() Font <- "sans" TitleFont <- 16 TickFont <- 12 AxisFont <- 14 plotformat <- theme(plot.title=element_text(family = Font, face="bold", size=TitleFont, vjust = 1))+ theme(axis.title.x=element_text(family = Font, face="bold", size=AxisFont, vjust = 0)) + theme(axis.title.y=element_text(family = Font, face="bold", size=AxisFont, vjust=0.9, hjust=1))+ theme(axis.text.x=element_text(angle=45, family = Font, size=TickFont, colour = "black", vjust=.5)) + theme(axis.text.y=element_text(family = Font, size=TickFont, colour = "black")) + theme(legend.title=element_text(family = Font, size=AxisFont)) + theme(legend.text=element_text(family = Font, size=TickFont), strip.text = element_text(size=TickFont)) ``` ```{r Figure 3} #Seasonal Bat Activity Rate Modeling ####Data and Model Preparation#### modelData <- EchoPitch %>% #biweekly data #Regrouping original data to ignore changes in treatment group_by(Site,Turbine,DetectorPosition,Year,Biweek) %>% summarise(across(DN:Bat_Passes,sum), .groups="drop") %>% mutate(Pass_Rate=Bat_Passes/DN, #bat pass rate DoY=yday(Biweek), #numerical day of year for start of 14-day period Turbine_ID=factor(paste(Site,Turbine)), #Unique ID for same number turbines Biweek=as.Date(Biweek) #adjusting for plot ) %>% select(Site,Year,Turbine,Turbine_ID,DetectorPosition,Pass_Rate,DoY,Biweek,DN,Bat_Passes) #No need to run pearson correlations because there aren't multiple numeric inputs. #For purposes of this analysis we are only interested in the #seasonal pattern and not underlying variation that day of year captures like #wind speed, temperature, etc. #Strong data bias near zero ggplot(modelData,aes(x=Pass_Rate)) + geom_histogram() + facet_grid(Site~DetectorPosition, scales="free") ####Gam Model#### model <- gam(Pass_Rate~ s(DoY) + #regular spline is valid in this analysis to capture #variation as data is already aggregated on a biweekly #basis which has a similar effect to a penalized spline #that is often necessary for smooth day of year variation DetectorPosition + Site+ Year+ s(Turbine_ID,bs="re") #added in random effect which smoothed the error , data=modelData, method="REML", family="tw" #Using Tweedie with auto-fit to capture biased #distribution near zero of continuous response variable ) ####Model Fit Checks#### FitCheck <- check(getViz(model)) FitCheck #QQ: OK #Residuals v Predictor: No overly concerning patterns #Residuals Histogram: Mostly normal distribution #edf: Non-significant knots #Response v Fitted: Fairly good linear estimations on par with values #Create Supplement of Model Fit Plots Supplement <- plot_grid(FitCheck[[1]]+theme_bw(), FitCheck[[2]]+theme_bw(), FitCheck[[3]]+theme_bw(), FitCheck[[4]]+theme_bw(), labels = c("A", "B", "C", "D") ) #Create a supplement plot from gam.check figures resid <- residuals(model,type="deviance") linpred <- napredict(model$na.action,model$linear.predictors) observed <- napredict(model$na.action,model$y) qq.gam(model, rep = 0, level = 0.9, type = "deviance", rl.col = 2, rep.col = "gray80") concurvity(model) #Autocorrelation not present for day of year #Quick DHARMa overview check resid <- simulateResiduals(model) plot(resid) #Outliers are not of concern based on other outputs and predictions #Data is not overdispersed #No significant quantiles detected but quickly making sure quantiles look ok per variable plotResiduals(resid,modelData$DoY) #slightly off lower quantile but fine plotResiduals(resid,modelData$Turbine_ID) #some low turbines but ok plotResiduals(resid,modelData$Site) #normal plotResiduals(resid,modelData$DetectorPosition) #normal #Taken together, nothing here indicates that predictions will be heavily influenced #and the model in general has an acceptable fit to interpret predictions from ####Model Results#### summary(model) #Mid-tower has significantly higher pass rates (expected) #No significant difference between sites (interesting, but they're also nearby) #Significant difference between years (unexpected, but impact is low 0.1 difference) #Significant seasonal variation in pass rate (expected) #Significant variation between turbines (interesting) ggplot(modelData %>% mutate(Turbine_ID=gsub(" ","\n",Turbine_ID)) %>% filter(DetectorPosition=="Nacelle"), aes(x=DoY,y=Pass_Rate,color=Year)) + geom_line(aes(group=Year)) + facet_wrap(~Turbine_ID) + ggtitle("Nacelle") #Makes sense, some turbines have increased activity at them, based on my experience #this is relatively normal at wind facilities though I'm not sure what the underlying #reasoning is (e.g. land cover, topography, etc.) ####Model predictions#### #Generate predictions predictions <- ggpredict(model, #Modeling only the period from February to November terms=c("DoY [50:330]", "Site","DetectorPosition" )) %>% data.frame() %>% `colnames<-`(c("DoY","Pass_Rate","Std.Error", "Conf.Low","Conf.High","Site", "DetectorPosition" )) %>% mutate(Biweek=as.Date(DoY,origin="2020-12-31", tz="Etc/GMT+6")) #Plot Figure3 <- ggplot(data=modelData,aes(x=Biweek,y=Pass_Rate,fill=Site)) + #Original Data Points geom_point(colour="black",pch=21, size=1.5) + #Facet facet_grid(rows=vars(DetectorPosition), #cols=vars(Year), scales="free_y") + #Data Point Fill scale_fill_manual(values=c("Kings Point"="red","North Fork Ridge"="blue", limits=force),"Recorded") + #Set New Fill for Predicted Values ggnewscale::new_scale_fill()+ #Fitted Values geom_line(size=0.75, data=predictions,aes(x=Biweek,y=Pass_Rate,color=Site)) + #Confidence Interval geom_ribbon(data=predictions,aes(ymin=Conf.Low,ymax=Conf.High, fill=Site),alpha=.3,color=NA)+ #Predicted Color and Fill scale_color_manual(values=c("Kings Point"="red","North Fork Ridge"="blue"), "Predicted")+ scale_fill_manual(values=c("Kings Point"="red","North Fork Ridge"="blue"), "Predicted") + #Axis Labels ylab("Bat Passes per Detector Night") + xlab("Date")+ #Scale X-Axis to Same Values as Seen in Later Histograms scale_x_date(breaks=sort(unique(modelData$Biweek))[1:21],date_labels="%b %d", limits=c(min(predictions$Biweek),max(predictions$Biweek))) + #Set Themes theme_bw() + plotformat + theme(axis.text.x = element_text(angle = 45, hjust=1, vjust=1.1), axis.title.y = element_text(hjust=.7)) Figure3 ggsave(plot = Figure3, filename = paste(out,"Figure3.pdf",sep="/"), width = 17, height = 9, units = "cm", dpi = 300) ggsave(plot=Supplement, filename = paste(out,"Supplement.pdf",sep="/"), width = 18, height = 17, units = "cm", dpi = 300) ``` ```{r Figure 4} #Number of searches conducted Carcasses %>% select(Site,Year,Biweek,N_Searches) %>% unique() %>% #both detector positions share the same estimates group_by(Site) %>% summarise(n=sum(N_Searches)) #Number of carcasses recovered Carcasses %>% select(Site,Year,Biweek,N_Carcasses) %>% unique() %>% #both detector positions share the same estimates group_by(Site) %>% summarise(n=sum(N_Carcasses)) #Average Exposure EchoPitch %>% #per height group_by(DetectorPosition) %>% summarise(N_M=sum(Measured), N_S=sum(Simulated), Rate=sum(Measured)/sum(DN)) EchoPitch %>% #overall summarise(N_M=sum(Measured), N_S=sum(Simulated), Rate=sum(Measured)/sum(DN)) #Comparison of measured exposure and fatality recovery #Generate Figures lapply(c("Nacelle","Mid-tower"),function(p){ #Plot Data plotData <- Carcasses %>% filter(DetectorPosition==p) #Adjust extended axis limits so figures are even limits <- data.frame(position=c("Nacelle","Mid-tower"), xmin=c(0,0), ymin=c(0,0), xmax=c(12,80), ymax=c(19,19) ) #Make Plot plot <- plotReg(Data=plotData, x="Exposed_Rate", y="Estimate_PerTurbine", xaxis="Biweekly Acoustic Exposure Rate", yaxis="Biweekly Fatality Estimate\nper Turbine", xlims=c(limits$xmin[limits$position==p], limits$xmax[limits$position==p]), ylims=c(limits$ymin[limits$position==p], limits$ymax[limits$position==p]), TFont=14, KFont=10, AFont=12 ) #Add Title plot <- plot + ggtitle(paste(p,"Detectors")) #Environment object name obj <- paste("plot",ifelse(p=="Nacelle","nac","mid"),sep="_") #Assign plot to environment assign(obj,plot,envir=.GlobalEnv) return() }) ####Bind Figures#### #Get legend legend <- get_legend(plot_nac+ theme(legend.direction = "horizontal")) #Join both plots without legends prow <- plot_grid( plot_nac + theme(legend.position = "none", plot.margin = margin(r = 20,t=1,b=1,l=1)), plot_mid + theme(legend.position = "none", plot.margin = margin(r = 20,t=1,b=1,l=1)), nrow=1, rel_widths=c(0.98,1) ) #Join combined plot to legend Figure4 <- plot_grid(prow, legend, nrow=2, rel_heights = c(1,0.2)) Figure4 ggsave(plot=Figure4, filename = paste(out,"Figure4.pdf",sep="/"), width = 17, height = 9, units = "cm", dpi = 300) remove(plot_nac,plot_mid,prow,legend) ``` ```{r Figure 5} #Seasonal Exposure Ranges EchoPitch %>% #per height group_by(Site,Year,Turbine,DetectorPosition) %>% summarise(across(DN:Simulated,sum), .groups="drop") %>% mutate(Simulated_Rate=Simulated/DN, Measured_Rate=Measured/DN) %>% group_by(DetectorPosition) %>% summarise(Min_Simualted=min(Simulated_Rate), Max_Simulated=max(Simulated_Rate), Min_Measured=min(Measured_Rate), Max_Measured=max(Measured_Rate)) #Generate Figures lapply(c("Nacelle","Mid-tower"),function(p){ #Plot Data plotData <- EchoPitch %>% filter(DetectorPosition==p) %>% group_by(Site,Year,Turbine,DetectorPosition) %>% summarise(across(DN:Simulated,sum), .groups="drop") %>% mutate(Measured_Rate=Measured/DN, Simulated_Rate=Simulated/DN) #Adjust extended axis limits so figures are even limits <- data.frame(position=c("Nacelle","Mid-tower"), xmin=c(0,10), ymin=c(0,10), xmax=c(5,60), ymax=c(5,60) ) #Make Plot plot <- plotReg(Data=plotData, x="Measured_Rate", y="Simulated_Rate", xaxis="Acoustic Exposure Based on\nRotor Speed", yaxis="Acoustic Exposure Based on\nCurtailment Criteria", xlims=c(limits$xmin[limits$position==p], limits$xmax[limits$position==p]), ylims=c(limits$ymin[limits$position==p], limits$ymax[limits$position==p]), TFont=14, KFont=10, AFont=12 ) #Add Title plot <- plot + ggtitle(paste(p,"Detectors")) #Environment object name obj <- paste("plot",ifelse(p=="Nacelle","nac","mid"),sep="_") #Assign plot to environment assign(obj,plot,envir=.GlobalEnv) return() }) ####Bind Figures#### #Get legend legend <- get_legend(plot_nac+ theme(legend.direction = "horizontal")) #Join both plots without legends prow <- plot_grid( plot_nac + theme(legend.position = "none"), plot_mid + theme(legend.position = "none", plot.margin=margin(5,10,5,5)), nrow=1, rel_widths=c(0.93,1) ) #Join combined plot to legend Figure5 <- plot_grid(prow, legend, nrow=2, rel_heights = c(1,0.2)) Figure5 ggsave(plot=Figure5, filename = paste(out,"Figure5.pdf",sep="/"), width = 17, height = 9.5, units = "cm", dpi = 300) remove(plot_nac,plot_mid,prow,legend) ``` ```{r Figure 6} #Simulated exposure rate under varying cut-in speeds Figure6_Data <- EchoPitch %>% pivot_longer(cols=Curt_3:Curt_8, names_to="Simulated_Treatment", values_to="Exposed_Passes") %>% mutate(Simulated_Treatment=recode(Simulated_Treatment, "Curt_3"="3.0 m/s", "Curt_4"="4.0 m/s", "Curt_4.5"="4.5 m/s", "Curt_5"="5.0 m/s", "Curt_5.5"="5.5 m/s", "Curt_6"="6.0 m/s", "Curt_6.5"="6.5 m/s", "Curt_7"="7.0 m/s", "Curt_7.5"="7.5 m/s", "Curt_8"="8.0 m/s"), Rate_Exposed=Exposed_Passes/DN) %>% #Clipping treatments to display filter(!Simulated_Treatment %in% c("3.0 m/s","8.0 m/s","7.5 m/s")) %>% #Clipping to associated date range as shown in figure1 filter(Biweek<=as.POSIXct("2021-12-01",tz="Etc/GMT+6")) %>% group_by(Simulated_Treatment,DetectorPosition,Biweek) %>% summarise(Rate_Total=round(sum(Exposed_Passes)/sum(DN)*100,1), #seasonal percent Rate_Average=round(mean(Rate_Exposed,na.rm=T),1), #mean N=length(Turbine[!is.na(Rate_Exposed)]), #sample size SD=sd(Rate_Exposed,na.rm=T), #standard deviation Margin=qt(0.975,df=N-1)*SD/sqrt(N), #error margin UL=Rate_Average+Margin,#upper limit LL=ifelse(Rate_Average-Margin<0,0,Rate_Average-Margin), #lower limit .groups="drop") %>% arrange(DetectorPosition,Simulated_Treatment) %>% mutate(Biweek=factor(format(Biweek,"%b %d"), levels=c(format(sort(unique(Biweek)),"%b %d")) ) ) #Available intervals Figure6_Data %>% select(Biweek,DetectorPosition,N) %>% unique() %>% summarise(N=sum(N)) #Biweekly simulated exposure Figure6 <- ggplot(Figure6_Data,aes(x=Biweek,y=Rate_Average,fill=Simulated_Treatment)) + geom_col(position=position_dodge2()) + geom_errorbar(aes(ymin=LL,ymax=UL),position=position_dodge2(),size=.35)+ facet_grid(DetectorPosition~., scales="free_y")+ scale_y_continuous(expand=expansion(mult=c(0,0.05)), #limits=c(0,NA) )+ labs(x="Date", y="Mean Acoustic Exposure Rate")+ theme_bw() + plotformat + theme(axis.text.x=element_text(angle=45,hjust=1,vjust=1.1), axis.title.y=element_text(hjust=.5), strip.text = element_text(size=12), panel.spacing = unit(1, "lines")) + scale_fill_manual(values=c("#EB1E2C", "#FD6F30", #"#EF8A0C", "#F9D23C", #"#4E9F50", "#87D180", "#3896C4", "#91DCEA", "#A26DC2", "#C799BC"), name="Cut-In Speed") Figure6 ggsave(plot=Figure6, filename = paste(out,"Figure6.pdf",sep="/"), width = 17, height = 9, units = "cm", dpi = 300) ``` ```{r Figure 7} #Simulated percent exposure under varying cut-in speeds Figure7_Data <- EchoPitch %>% pivot_longer(cols=Curt_3:Curt_8, names_to="Simulated_Treatment", values_to="Exposed_Passes") %>% mutate(Simulated_Treatment=recode(Simulated_Treatment, "Curt_3"="3.0 m/s", "Curt_4"="4.0 m/s", "Curt_4.5"="4.5 m/s", "Curt_5"="5.0 m/s", "Curt_5.5"="5.5 m/s", "Curt_6"="6.0 m/s", "Curt_6.5"="6.5 m/s", "Curt_7"="7.0 m/s", "Curt_7.5"="7.5 m/s", "Curt_8"="8.0 m/s"), #Overwrite when DN>0 and Passes=0 Percent_Exposed=ifelse(DN>0 & Passes_Ops==0,0,Exposed_Passes/Passes_Ops*100)) %>% #Clipping treatments to display filter(!Simulated_Treatment %in% c("3.0 m/s","8.0 m/s","7.5 m/s")) %>% #Clipping to associated date range as shown in figure1 filter(Biweek<=as.POSIXct("2021-12-01",tz="Etc/GMT+6")) %>% group_by(Simulated_Treatment,DetectorPosition,Biweek) %>% summarise(Percent_Total=round(sum(Exposed_Passes)/sum(Passes_Ops)*100,1), #seasonal rate Percent_Average=round(mean(Percent_Exposed,na.rm=T),1), #mean N=length(Turbine[!is.na(Percent_Exposed)]), #sample size SD=sd(Percent_Exposed,na.rm=T), #standard deviation Margin=qt(0.975,df=N-1)*SD/sqrt(N), #error margin UL=ifelse(Percent_Average+Margin>100,100,Percent_Average+Margin), #upper limit LL=ifelse(Percent_Average-Margin<0,0,Percent_Average-Margin), #lower limit .groups="drop") %>% arrange(DetectorPosition,Simulated_Treatment) %>% mutate(Biweek=factor(format(Biweek,"%b %d"), levels=c(format(sort(unique(Biweek)),"%b %d")) ) ) #Number of available intervals Figure7_Data %>% select(Biweek,DetectorPosition,N) %>% unique() %>% summarise(N=sum(N)) #Biweekly simulated exposure Figure7 <- ggplot(Figure7_Data,aes(x=Biweek,y=Percent_Average,fill=Simulated_Treatment)) + geom_col(position=position_dodge2()) + geom_errorbar(aes(ymin=LL,ymax=UL),position=position_dodge2(),size=.3)+ facet_grid(DetectorPosition~., scales="free_y")+ scale_y_continuous(expand=expansion(mult=c(0,0.02)), limits=c(0,100) )+ labs(x="Date", y="Mean Acoustic Exposure (%)")+ theme_bw() + plotformat+ theme(axis.text.x=element_text(angle=45,hjust=1,vjust=1.1), axis.title.y=element_text(hjust=.5), strip.text = element_text(size=12), panel.spacing = unit(1, "lines")) + scale_fill_manual(values=c("#EB1E2C", "#FD6F30", #"#EF8A0C", "#F9D23C", #"#4E9F50", "#87D180", "#3896C4", "#91DCEA", "#A26DC2", "#C799BC"), name="Cut-In Speed") Figure7 ggsave(plot=Figure7, filename = paste(out,"Figure7.pdf",sep="/"), width = 17, height = 10, units = "cm", dpi = 300) ```