R/CollisionModelv4.R

In josephstatwick/EagleModelCalc: Eagle fatality calculator for wind farms

#' @title Eagle Fatality model
#' @author Joe Statwick (Arcadis)
#' @author Mark Otto (FWS)
#'
#' @description This function calculates and plots the annual predicted fatalities at a location.
#' @param cProject The name of the project, e.g., Big Farm 2
#' @param Name The name of the modeler, e.g., Jane Doe
#' @param nTurbine the number of turbines
#' @param HazDiam diameter of turbine blades (in m) (optional)
#' @param HazRadm radius of turbine blades (in m) (optional)
#' @param HazRadKm radius of turbine blades (in km) (optional)
#' @param HzKM2 Hazardous area around each turbine (in km^2) (can be supplied directly or calculated via the optional arguments HazDiam, HazRadm, or HazRadKm)
#' @param CntMin Duration of one pre-construction observation, in minutes (optional)
#' @param CntHr Duration of one pre-construction observation, in hours (can be supplied directly or calculated via the optional argument CntMin)
#' @param EOutMin Number of Eagle carcasses observed during PCMM surveys. (Defaults to -1 if PCMM not performed)
#' @param PCMMsurveys Number of turbines surveyed during each PCMM survey event. Can be prorated for partially-surveyed turbines. (Defaults to 0 if PCMM not performed)
#' @param PCMMhrs Number of daylight hours in PCMM survey period (Defaults to 0 if PCMM not performed)
#' @param Oprate Average annual operational rate of turbines. Defaults to 1
#' @param CprHrKM2 Total surveyed hours * hazardous area surveyed in post-construction mortality monitoring (Can be supplied directly or calculated via the optional arguments PCMMsurveys and PCMMhrs)
#' @param ExpSvy A data frame containing strata-wise (seasonal) data for eagle minutes (EMin), number of surveys (nCnt), size of a survey (CntKM2), and daylight hours (DayLtHr). Can be created manually or generated automatically by EagleWizard(). EMin, nCnt, and CntKM2 default to -1, 0, and 0 if pre-construction surveys were not performed, but DayLtHr must be supplied. Use DayLen() to calculate annual or seasonal daylight hours.
#' @param CPrPriors Collision probability priors, either USFWS or Bay. Defaults to USFWS. If neither, custom priors must be provided.
#' @param aPriCpr Alpha parameter for the exposure prior. Must be supplied if neither USFWS nor Bay priors are selected.
#' @param bPriCpr Beta parameter for the exposure prior. Must be supplied if neither USFWS nor Bay priors are selected.
#' @param AddTot Add strata for total (TRUE) or not (FALSE)
#' @export
#' @import rv


### Draft USFWS Collision Fatality Model Code version 4.1 (23 Apr 2013)               ###
### This code is a working version only. It is not intended for general distribution. ###
### Check back often for updates to the model/code                                    ###

## Please review Model Code v4.1 README for general instructions
## requires the rv and maptools package for R

## Analysis Inputs ##
EagleModel<-function(cProject, Name, nTurbine, HazDiam, HazRadm = HazDiam/2, HazRadKm = (HazRadm/1000), HzKM2 = (nTurbine*pi*HazRadKm^2), CntMin = NULL, CntHr = (CntMin/60), EOutMin = -1, PCMMsurveys = 0, PCMMhrs = 0, Oprate = 1, CprHrKM2 = (PCMMsurveys*PCMMhrs)*HzKM2/nTurbine*Oprate, ExpSvy=data.frame(row.names=c("Annual"),EMin=-1, nCnt=0, CntKM2= 0, DayLtHr=NULL), CPrPriors = "USFWS", aPriCPr = NULL, bPriCPr = NULL, AddTot = FALSE){

  UCI<-c(0.5,0.8,0.9,0.95)
  nSim<-100000
  setnsims(nSim)

### Survey Inputs ###

 nSvy<-nrow(ExpSvy)
 cSvy<-(rownames(ExpSvy))

 SmpHrKM2<-with(ExpSvy,nCnt*CntHr*CntKM2)
 ExpFac<-ExpSvy$DayLtHr*HzKM2*Oprate

# Calculate the fatalities and store as a temporary object.
 if(!is.null(aPriCPr)){
 tmp<-with(ExpSvy,mapply(simFatalCPr,EMin=EMin,EOutMin=EOutMin,SmpHrKM2=SmpHrKM2,CprHrKM2=CprHrKM2, ExpFac=ExpFac, CPrPriors=CPrPriors, aPriCPr=aPriCPr, bPriCPr=bPriCPr, SIMPLIFY=FALSE))
 } else {
   tmp<-with(ExpSvy,mapply(simFatalCPr,EMin=EMin,EOutMin=EOutMin,SmpHrKM2=SmpHrKM2,CprHrKM2=CprHrKM2, ExpFac=ExpFac, CPrPriors=CPrPriors, SIMPLIFY=FALSE))
 }



# R code to get the survey specific simulations in an rv vector.
 Fatalities<-rvnorm(nSvy)
 Exp<-data.frame(Mean=rep(NA,nSvy),SD=NA,row.names=cSvy)
 for(i in 1:nSvy){
  Fatalities[i]<-tmp[[i]]
  Exp[i,]<-attr(tmp[[i]],"Exp")
 }
 rm(tmp)
 names(Fatalities)<-cSvy

# Summarize the results, including a total if needed.
 nSvy<-length(Fatalities)
 if(is.null(nSvy))nSvy<-1
 FatalStats<-RVSmry(cSvy,Fatalities,probs=UCI)
 if(AddTot){
  FatalStats<-rbind(
   FatalStats,
   RVSmry("Total",sum(Fatalities),probs=UCI)
  )
 }

# Look at the results
 cat(cProject,"\n")
 cat(paste(Name,", ",date(),"\n",sep=""))
 #Number of Turbines
 print(nTurbine)
 #Hazardous Area Per Turbine (km^2)
 #print(HzKM2PT)
 print(ExpSvy)
 #Exposure rate
 print(Exp,digits=3)
 #Annual Collision Fatalities
 print(FatalStats,digits=2)
 print(c(EMin=ExpSvy$EMin, EOutMin=EOutMin, SmpHrKM2=SmpHrKM2, CprHrKM2=CprHrKM2, ExpFac=ExpFac))

# Plots
 nPlot<-nSvy+as.integer(AddTot)
 nCol<-floor(sqrt(nPlot))
 nRow<-ceiling(nPlot/nCol)
 xlim<-range(rvrange(Fatalities))

 par(mfrow=c(nRow,nCol))
 for(iPlot in 1:nSvy){
  plotFatal(Fatalities[iPlot],probs=UCI,
#  xlim=xlim,add=FALSE,  # uncomment this line to put the graphs for all of the strata on the same scale
  main=cSvy[iPlot])
 }
 if(AddTot)plotFatal(sum(Fatalities),main="Total")


}