--- title: "Calculating avoidance rates for collision risk modelling" author: "Luke Ozsanlav-Harris" date: "23/02/2022" output: html_document --- ```{r eval=FALSE, include=FALSE} #Tips for me to make this R markdown #This web page has lots of good tips: https://www.earthdatascience.org/courses/earth-analytics/document-your-science/intro-to-the-rmarkdown-format-a#nd-knitr/ ``` ## Project aim: Calculate avoidance rates for seabirds for use within collison risk modelling This is for key species: northern gannet, black-legged kittiwake, lesser black-backed gull, great black-backed gull, herring gull, Sandwich tern, and appropriate species groups such as all terns, small gulls, large gulls and all gulls Avoidance rates will be calculated using the basic and extended band models, using the stochastic version of both of these to calculate more robust CIs around estimates ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ### Load the required packages ```{r Load-packages, message=FALSE, warning=FALSE} library(lubridate) library(truncnorm) library(foreach) library(doParallel) library(ggplot2) library(sf) library(tidyverse) library(ggspatial) library(png) library(gghalves) library(htmltools) library(htmlwidgets) library(webshot) library(sparkline) library(formattable) library(magrittr) library(knitr) ``` ## 1. Read in the data sets This includes the observed passage and collision data as well as the modeled flight height distributions from Jonstone *et al* (2014) At the end it would be good to show the Avoidance rates by species compared to those calculated by Cook (2021) ```{r Read-in-data, warning=FALSE} ## Set working directory where data & outputs will be stored options("scipen"=100, "digits"=8) ## read in data coll.dat = read.table("Ozsanlav-Harris et al (2022) Collision Data.csv", header = T, sep = ",") ### collated collision & survey data flightheight = read.table("Flight Height Data.csv", header = T, sep = ",") ### modeled proportion of birds at collision risk height BH.boot = read.table("BH_Boot.csv", header = T, sep = ",") ### Black-headed Gull Flight height bootstraps CG.boot = read.table("CG_Boot.csv", header = T, sep = ",") ### Common Gull Flight height bootstraps CT.boot = read.table("CT_Boot.csv", header = T, sep = ",") ### Common Tern Flight height bootstraps GB.boot = read.table("GB_Boot.csv", header = T, sep = ",") ### Great Black-backed Gull Flight height bootstraps HG.boot = read.table("BH_Boot.csv", header = T, sep = ",") ### Herring Gull Flight height bootstraps KI.boot = read.table("KI_Boot.csv", header = T, sep = ",") ### Kittiwake Flight height bootstraps LB.boot = read.table("LB_Boot.csv", header = T, sep = ",") ### Lesser Black-backed Gull Flight height bootstraps NE.boot = read.table("NE_Boot.csv", header = T, sep = ",") ### Sandwich Tern Flight height bootstraps LG.boot = read.table("LG_Boot.csv", header = T, sep = ",") ### Little Gull Flight height bootstraps ``` ------ ## 2. Create the functions used in collsion risk modelling The code in the functions will not be given to minimize space but each of the main functions and there use will be given #### Function 1- probability of collisions **pcoll()** : Calculate the probability that a bird collides with a rotor if it flies through the rotor sweep zone for the basic band model The flux used alongside the collision probability from this function is the percentage of birds passing through the rotor sweep zone derived from the observed percentage of birds at rotor height data and the area of the rotor sweep zone compared to the total area at collision risk height #### Function 2- collision integral **coll.int()**: Calculate the average collision integral/probability, this allows the probability of collision to change along the length of a rotor blade The flux used alongside the collision probability from this function is just the total flux, but the proportion of flux at rotor height and at each point long the blade is calculated using the flight height distributions in Jonstone *et al* (2014) #### Function 3 - daylength calcualtor **nhoursday()**: Calculates the length of day and night in hours given a start date, end data and latitude #### Function 4 - standard devation around a ratio **std.dev()**: Delta method for estimating SD around a ratio ```{r Create-functions-used-in-avoidance-calcualtions, include=FALSE} ## #### 2.1 Function for probability of collision from Basic Band model #### ## ## use k =1 & 3 blades as default ## Explaining the variables- ## Aspect: bird length/wingspan ## max.chord: width of turbine blade ## chord: ## rotor.period: rotor speed in rpm ## pitch: pitch of the rotor balds ## prop.flights.upwind: ## radius: ## rotor.rad: radius of the rotor ## prop.flights.upwind: ## This function calculate the probability that a bird collides with a rotor if it flies through the rotor sweep zone for the basic band model ## This is dependent on the turbines itself and the morphology and flight speed of the bird. pcoll = function(k =1, no.blades = 3, max.chord, pitch, rotor.rad, rotor.period,wingspan, bird.length, aspect, flap.glide, bird.speed,prop.flights.upwind,radius,chord) { ## p(collision) data alpha = bird.speed * rotor.period / (radius * 2 * rotor.rad * pi) alphaLTaspect = which(alpha 1] = 1 overall_up_p_collision = 2 * (sum(radius[1:19] * up_p_collision[1:19]) + up_p_collision[20]/2) * 0.05 ## downwind p(collision) down_col_length = abs(-max.chord * chord * sin(pitch*pi/180) + (alpha * max.chord * chord * cos(pitch*pi/180))) down_col_length2 = (wingspan * alpha) + down_col_length down_col_length2[alphaLTaspect] = bird.length + down_col_length[alphaLTaspect] down_p_collision = (no.blades/rotor.period)*down_col_length2/bird.speed down_p_collision[down_p_collision > 1] = 1 overall_down_p_collision = 2 * (sum(radius[1:19] * down_p_collision[1:19]) + down_p_collision[20]/2) * 0.05 return(((1-prop.flights.upwind) * overall_down_p_collision) + (prop.flights.upwind * overall_up_p_collision)) } ## ## #### 2.2 Estimate the number of hours daylight & night during each study #### ## ## ## Assuming it uses trigonometry and latitude to calculate sun angle and the day length ## There are R package that clauclate this for you nhoursday = function(start, end, latitude) { day.seq = seq(start, end, by = "day") day.seq = yday(day.seq) DayLength = vector() NightLength = vector() for(zz in 1:length(day.seq)){ P = asin(0.39795 * cos (0.2163108 + 2 * atan (0.9671396 * tan(0.0086 * (day.seq[zz] - 186))))) DayLength[zz] = 24 - (24/pi) * acos((sin(0.8333*pi/180) + sin(latitude * pi / 180) * sin (P)) / (cos(latitude *pi /180 ) * cos(P))) NightLength[zz] = 24 - DayLength[zz] } DayLength = sum(DayLength) NightLength = sum(NightLength) return(list(DayLength, NightLength)) } ## #### 2.3 Delta method for estimating SD around a ratio #### ## std.dev = function(b,c){ var.b <- var(b) var.c <- var(c) cov.bc <- cov(b,c) a <- sum(b)/sum(c) pd.a.b <- 1/sum(c) pd.a.c <- -sum(b)/sum(c)^2 pd.mat <- matrix(c(pd.a.b, pd.a.c), nrow=1) vcv.mat <- matrix(c(var.b, cov.bc, cov.bc, var.c), nrow = 2, byrow = TRUE) mat1 <- pd.mat %*% vcv.mat var.a <- mat1 %*% t(pd.mat) return(sqrt(var.a)) } ## ## #### 2.4 Functions for extended Band Collision Risk Model #### ## ## ## Calculate collision probability at any given point along the rotor blade pcollex = function (r,phi, updown, Rotor.Diameter, Blade.Width, Rotor.Speed, Pitch, flapglide, wingspan, birdlength, speed, Blades,coverC = coverC) { if(r>1) cell_val =21 if(r<=1) which(coverC$rad == abs(round(ceiling(r*100)/5) * 5)/100) -> cell_val upper = coverC[cell_val,1] lower = coverC[cell_val - 1, 1] p = (r - lower) / (upper-lower) c = coverC[cell_val-1, 2] + p * (coverC[cell_val,2] - coverC[cell_val-1,2]) radius = Rotor.Diameter/2 * r if(radius ==0) return(1) if(radius !=1) ({ chord = Blade.Width * c omega = Rotor.Speed * 2 * pi /60 pitch = Pitch * pi /180 phi = phi * pi /180 if(updown == "up") direction = 1 if(updown != "up") direction = -1 if(flapglide == "flap") wingspan-> Wingspan2 if(flapglide != "flap") wingspan * abs(cos(phi))-> Wingspan2 multiplier = Blades * omega /(2 * pi * speed) alpha = speed/(radius * omega) CollideLength_1 = abs(direction * chord * sin(pitch) + alpha * chord * cos(pitch)) if(birdlength > Wingspan2 * alpha) birdlength -> CollideLength2 if(birdlength <= Wingspan2 * alpha) Wingspan2 * alpha -> CollideLength2 return(multiplier * (CollideLength_1 + CollideLength2)) }) } ## calculate collision integral to account for risk and birds both varying with height ## Function 1 that feeds into collision integral pcollxy = function(x,y,updown, Rotor.Diameter, Blade.Width, Rotor.Speed, Pitch, flapglide, wingspan, birdlength, speed, Blades,coverC = coverC) { r = (x * x + y * y) ^ 0.5 if( y == 0 & x >= 0) pi /2 -> phi if( y == 0 & x < 0) -pi/2 -> phi if(y != 0) atan (x/y) -> phi if( y <0) (phi + pi)*180/pi -> phi2 if(y >= 0) phi*180/pi -> phi2 pcollex(r,phi2,updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) } ## Function 2 that feeds into collision integral xrisksum2 = function (y, xinc,updown, Rotor.Diameter, Blade.Width, Rotor.Speed, Pitch, flapglide, wingspan, birdlength, speed, Blades,coverC = coverC) { xmax = (1- y * y) ^ 0.5 imax = as.integer (xmax/xinc) risk = (pcollxy (imax * xinc, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 + pcollxy(xmax, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC)/2) * (xmax - imax * xinc) if(imax<=0) return(2*risk) if(imax>0)({ risk2 = risk + (pcollxy(0,y,updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 + pcollxy(imax * xinc, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC)/2) * xinc for (i in 1: (imax - 1)) { risk2 = risk2 + pcollxy(i * xinc, y, updown, Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) * xinc } return(2*risk2) }) } ## collison integral function coll.int = function(HD.y,HD.d.y, FluxInt, Blades, HubHeight, Rotor.Diameter, Blade.Width, Rotor.Speed,Pitch,speed,flapglide,birdlength,wingspan, FH.dat, Prop_Upwind,coverC = coverC){ CollMinUP = HD.d.y[1] * xrisksum2(HD.y[1], 0.05, "up", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at lowest point on rotor blade CollIntUP = CollMinUP + HD.d.y[41] * xrisksum2(HD.y[41], 0.05, "up", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at highest point on rotor blade for (i in 2:40) {CollIntUP = CollIntUP + HD.d.y[i] * xrisksum2(HD.y[i], 0.05, "up", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) } #### Fill in intermediate heights CollIntUP = CollIntUP * 0.05 * (2/pi) ## Down wind CollMinDown = HD.d.y[1] * xrisksum2(HD.y[1], 0.05, "down", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at lowest point on rotor blade CollIntDown = CollMinDown + HD.d.y[41] * xrisksum2(HD.y[41], 0.05, "down", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) / 2 ### risk at highest point on rotor blade for (i in 2:40) { CollIntDown = CollIntDown + HD.d.y[i] * xrisksum2(HD.y[i], 0.05, "down", Rotor.Diameter=Rotor.Diameter, Blade.Width=Blade.Width, Rotor.Speed=Rotor.Speed, Pitch = Pitch, flapglide = flapglide, wingspan=wingspan, birdlength=birdlength, speed=speed, Blades=Blades,coverC = coverC) } #### Fill in intermediate heights CollIntDown = CollIntDown * 0.05 * (2/pi) ## Average Collision Integral CollInt = (Prop_Upwind * CollIntUP) + ((1-Prop_Upwind) * CollIntDown) return(CollInt) } ## ## #### 2.5 Functions for Christie & Urquhart extension to Band model #### ## ## #### Step 1 PrC = function(Radius, Blades, pitch, chord.width, angvel, length, wingspan, delta, BF10G11, WGamma, WBSpeed, RelRad, RelChord, theta, c, s, Tarc, RotorV) { ## vector velocity of bird vbx = -WBSpeed * sin(delta) vby = 0 vbz = WBSpeed * cos(delta) ## vector of chord at theta = 0 cx0 = RelChord[1] * chord.width * cos(pitch) cy0 = 0 cz0 = RelChord[1] * chord.width * sin(pitch) ## vector length of bird blx = -length * sin (WGamma) bly = 0 blz = length * cos(WGamma) ## vector wingspan of gliding bird bwgx = wingspan * cos(WGamma) bwgy = 0 bwgz = wingspan * sin (WGamma) ## vector velocity of rotor vrx = RotorV[1] * c vry = -RotorV[1] * s vrz = 0 ## vector velocity of bird relative to blade u = vbx-vrx v = vby-vry w= vbz - vrz if (w ==0 | is.na(w)) w = 0.000000001 ## so we don't divide by 0 ## vector of chord rotated through theta cx = cx0 * c cy = -s * cx0 cz = cz0 l = cx+blx m = cy + bly n = cz + blz L1 = abs(((w*l-u*n)*c+(v*n-w*m)*s)/w) l2 = cx-blx m2= cy-bly n2 = cz-blz L2 = abs(((w*l2-u*n2)*c+(v*n2-w*m2)*s)/w) l3 = cx + bwgx m3 = cy+bwgy n3 = cz + bwgz L3 = abs(((w*l3-u*n3)*c+(v*n3-w*m3)*s)/w) l4 = cx-bwgx m4 = cy-bwgy n4 = cz-bwgz L4 = abs(((w*l4-u*n4)*c+(v*n4-w*m4)*s)/w) ## wingspan of a flapping bird rotated p = c* cos(WGamma) q = -s * cos(WGamma) r = c*sin(WGamma) k = wingspan/(p^2+q^2+r^2)^0.5 bwfx = k*p bwfy = k*q bwfz = k*r l5 = cx+bwfx m5 = cy+bwfy n5 = cz+bwfz L5 = abs(((w*l5-u*n5)*c+(v*n5-w*m5)*s)/w) l6 = cx-bwfx m6 = cy-bwfy n6 = cz-bwfz L6 = abs(((w*l6-u*n6)*c+(v*n6-w*m6)*s)/w) ## is bird flapping or gliding? # 0 = flapping # 1 = gliding if(BF10G11 == 0) ColLength = pmax(L1, L2, L3, L4, L5, L6) if(BF10G11 == 1) ColLength = pmax(L1, L2, L3, L4) Ptot = min(ColLength[1]/Tarc[1], 1) + (2*min(ColLength[2]/Tarc[1], 1)) + (2*min(ColLength[3]/Tarc[1], 1)) + (2*min(ColLength[4]/Tarc[1], 1)) + (2*min(ColLength[5]/Tarc[1], 1)) + (2*min(ColLength[6]/Tarc[1], 1)) + (min(ColLength[7]/Tarc[1], 1)) PrC = Ptot/12 return(PrC) } ################ Step 2 mean.coll.risk = function(RelRad, RelChord, simpsons.weight, theta, c,s,Radius, Blades, pitch, chord.width, period, length, wingspan, speed, BF10G11, min.angle, max.angle, wind.speed){ angle.head = seq(min.angle,max.angle,5) cor.pcoll = vector() for(abc in 1:length(angle.head)) { ## transform parameters for model delta = angle.head[abc]*pi/180 WGamma = delta - asin(wind.speed/speed*sin(delta)) WBSpeed = (wind.speed^2 + speed^2 - 2*wind.speed * speed *cos(WGamma))^0.5 angvel = 2*pi/period Tarc = 2* pi * RelRad * Radius/Blades RotorV = RelRad * Radius * angvel col.vec = vector() for(i in 1:20){ col.vec[i] = 2 * 0.05 *RelRad[i] *PrC(Radius, Blades, pitch, chord.width,angvel, length, wingspan,delta,BF10G11, WGamma, WBSpeed, RelRad[i],RelChord[i],theta,c, s, Tarc[i], RotorV[i])} ######## average collision risk along blade, weighted by Simpson's Rule to better approximate ara under the curve col.entry = sum(simpsons.weight*col.vec)/3 ######## Correct to reflect that area visible to bird as it approaches ######## is smaller than area visible to bird approaching perpendicularly ## calculate frontal area AFront = pi * Radius^2 ## calculate area covered by side of the turbine ASide = (mean(RelChord) * chord.width * sin(pitch) + wingspan) * 2 * Radius ## Calculate the area visible AVis = AFront * abs(cos(delta)) + ASide * abs(sin(delta)) ## now correct probability of collision cor.pcoll[abc] = col.entry * AVis/AFront } return(cor.pcoll) } #### End of functions #### ``` ------ ## 3. Calcualte length of day and night for all wind farms during the study period This allows the total flux through the windfarm to be calculated during the study period based off of the hourly passage rate ```{r Lenghth-of-day-and-night} ## format start date and end date into dates coll.dat$start.date = as.Date(coll.dat$start.date, format = "%d/%m/%Y") coll.dat$end.date = as.Date(coll.dat$end.date, format = "%d/%m/%Y") ##### create columns for number of ours daylight & night in each survey coll.dat$HoursDay = 0 coll.dat$HoursNight = 0 ### now use this to estimate n hours daylight & night for(i in 1:nrow(coll.dat)) { if(!is.na(coll.dat$end.date[i]) )({ hours = nhoursday(coll.dat$start.date[i], coll.dat$end.date[i],coll.dat$Latitude[i]) coll.dat$HoursDay[i] = hours[1] coll.dat$HoursNight[i] = hours[2] }) } ``` ------ ## 4. Calculate predicted number of collsions for Basic Band model (2012) This section of code calculates the number of collisions in the absence of avoidance for the basic band model ```{r Basic-band-predicted-collions} ## need to get area of survey window at collision risk height ## calculate area of survey window that is at CRH coll.dat$area.survey.window = coll.dat$width.of.survey.area * coll.dat$rotor.diameter ## calculate the total rotor sweep zone for all turbines in the study coll.dat$total.rotor.frontal.area = coll.dat$number.of.turbines * pi * (coll.dat$rotor.diameter/2)^2 ## For sites where we do not have a proportion of birds at risk height from site-specific survey data, we need to estimate this from the ## Johnston et al. modeled distributions ## initiate column coll.dat$FH.modelled = 0 ## loop through each row in the main data set with passage & collions data for(i in 1:nrow( coll.dat)) { ## extract the modeled flight height data for the appropriate species ## This data set is in 1 meter increments, therefore the row number just correspond to a number of meters sp.fh.dat = flightheight[ , grepl( paste(coll.dat$Species.for.Flight.Height.Distribution[i], ".est", sep = ""), names( flightheight ) ) ] min.rotor.height = coll.dat$hub.height[i] - (0.5 * coll.dat$rotor.diameter[i]) max.rotor.height = coll.dat$hub.height[i] + (0.5 * coll.dat$rotor.diameter[i]) ## calculate the proportion of birds within the rotor sweep zone based off of the modeled flight height values coll.dat$FH.modelled[i] = sum(sp.fh.dat[min.rotor.height:max.rotor.height]) } #### get proportion of birds at risk height for sites with no site-specific data coll.dat$X..flights.at.rotor.height[is.na(coll.dat$X..flights.at.rotor.height)] = coll.dat$FH.modelled[is.na(coll.dat$X..flights.at.rotor.height)] #### Now work out the probability of collision at each site ## vectors needed for the Band model radius = seq(0.05, 1, 0.05) chord = c(0.73, 0.79, 0.88, 0.96, 1, 0.98, 0.92, 0.85, 0.80, 0.75,0.7, 0.64,0.58,0.52,0.47,0.41,0.37,0.30,0.24, 0) ## this applies function 2.1 to each row in coll.dat ## calculate the probability that a bird collides with a rotor if it flies through the rotor sweep zone for(i in 1:nrow(coll.dat)){ coll.dat$Pcoll[i] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = coll.dat$rotor.pitch[i], rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/coll.dat$rotor.speed[i],wingspan = coll.dat$wingspan[i], bird.length = coll.dat$BirdLength[i], aspect = coll.dat$BirdLength[i]/coll.dat$wingspan[i], flap.glide = 0, bird.speed = coll.dat$flight.speed[i],prop.flights.upwind = 0.5,radius,chord) } ########## Now need to add column to main data set shown total number of birds passing at collision risk height coll.dat$HoursDay = as.numeric(coll.dat$HoursDay) coll.dat$HoursNight = as.numeric(coll.dat$HoursNight) ## re-calculate passage rate for the rows that already had a passage rate entered in excel coll.dat$passage.rate[!is.na(coll.dat$Number.of.birds.recorded)] = coll.dat$Number.of.birds.recorded[!is.na(coll.dat$Number.of.birds.recorded)]/coll.dat$N.Hours.monitoring[!is.na(coll.dat$Number.of.birds.recorded)] ### re-calculate calculate total flux by adding passage rates*number of hours during the day and night (had a correction of 0.25 or 0 depending on species) coll.dat$total.flux[!is.na(coll.dat$passage.rate)] = (coll.dat$passage.rate[!is.na(coll.dat$passage.rate)] * coll.dat$HoursNight[!is.na(coll.dat$passage.rate)] * coll.dat$Nocturnal.Correction[!is.na(coll.dat$passage.rate)]) + (coll.dat$passage.rate[!is.na(coll.dat$passage.rate)] * coll.dat$HoursDay[!is.na(coll.dat$passage.rate)] ) ## for rows with density estimates calculate total flux with survey area/density/flight speed/Number of hours daylight coll.dat$total.flux[!is.na(coll.dat$density)] = (3600/(coll.dat$width.of.survey.area[!is.na(coll.dat$density)]/coll.dat$flight.speed[!is.na(coll.dat$density)])) * coll.dat$HoursDay[!is.na(coll.dat$density)] * (coll.dat$density[!is.na(coll.dat$density)] * coll.dat$Survey.area[!is.na(coll.dat$density)]) ## NOTE: i think this is wrong to multiply by * (coll.dat$total.rotor.frontal.area/coll.dat$area.survey.window) ## this is a ratio between the width of the survey window and the toal frontal area of the rotor sweep zones ## In a turbines array the width of the survey area may only accommodate a portion of the turbines as they are in rows ## But to calcualte total rotor frontal area all the turbines in the array are used ## I think this then makes the ratio calculate too large, i.e. the proportion of the survey width that is in the sweep zone is too high ## Calculate total flux within the rotor sweep zone coll.dat$flux.CRH = coll.dat$total.flux * coll.dat$X..flights.at.rotor.height * (coll.dat$total.rotor.frontal.area/coll.dat$area.survey.window) ## calculte the number of collsions in the absence of avoidance by multiplying the total flux through the rotor sweep zone by the probability of colising with the rotor coll.dat$PredColl.NoAvoid.Opt1 = coll.dat$Pcoll * coll.dat$flux.CRH ``` ## 5. Summarise basic band model avoidance rates for species and species groups These plots show the avoidance rates at the species level and some key groupings for the basic band model The confidence intervals here just represent the variation between data set rows in avoidance rates (a row can be one wind farm over multiple years, one wind farm over one season or 1 turbine within a wind farm) ```{r Species-level-avoidance-basic-band, echo=FALSE, message=FALSE, warning=FALSE} ### get avoidance rates summarized by species & group spp = unique(coll.dat$Species) Basic.Band.Spp.Avoid = vector() Basic.Band.Spp.AvoidSD = vector() Basic.Band.Spp.AvoidUpCI = vector() Basic.Band.Spp.AvoidLowCI = vector() ## now for each unique value in the species column calculate ratio of estimated collision to predicted collisions (calculated above) for(i in 1:length(spp)) { Basic.Band.Spp.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == spp[i]])) Basic.Band.Spp.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$Species == spp[i]]) Basic.Band.Spp.AvoidUpCI[i] = Basic.Band.Spp.Avoid[i] + (1.96 * Basic.Band.Spp.AvoidSD[i]) Basic.Band.Spp.AvoidLowCI[i] = Basic.Band.Spp.Avoid[i] - (1.96 * Basic.Band.Spp.AvoidSD[i]) } ## create a table of the output Basic.Band.Avoid.Spp.Table = data.frame(Species = spp, Avoidance = Basic.Band.Spp.Avoid, AvoidanceSD = Basic.Band.Spp.AvoidSD, AvoidanceLCL = Basic.Band.Spp.AvoidLowCI, AvoidanceUCL = Basic.Band.Spp.AvoidUpCI) ``` ```{r Group-level-avoidance-basic-band, echo=FALSE, message=FALSE, warning=FALSE} ## now do the same as above but for species groups group = unique(coll.dat$group) Basic.Band.Group.Avoid = vector() Basic.Band.Group.AvoidSD = vector() Basic.Band.Group.AvoidUpCI = vector() Basic.Band.Group.AvoidLowCI = vector() for(i in 1:length(group)) { Basic.Band.Group.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == group[i]])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == group[i]])) Basic.Band.Group.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group == group[i]], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group == group[i]]) Basic.Band.Group.AvoidUpCI[i] = Basic.Band.Group.Avoid[i] + (1.96 * Basic.Band.Group.AvoidSD[i]) Basic.Band.Group.AvoidLowCI[i] = Basic.Band.Group.Avoid[i] - (1.96 * Basic.Band.Group.AvoidSD[i]) } ## think this might be fixing some kind of error Basic.Band.Group.Avoid[4] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group %in% c("small gull", "large gull", "gull")])) Basic.Band.Group.AvoidSD[4] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")], coll.dat$PredColl.NoAvoid.Opt1[coll.dat$group %in% c("small gull", "large gull", "gull")]) Basic.Band.Group.AvoidUpCI[4] = Basic.Band.Group.Avoid[4] + (1.96 * Basic.Band.Group.AvoidSD[4]) Basic.Band.Group.AvoidLowCI[4] = Basic.Band.Group.Avoid[4] - (1.96 * Basic.Band.Group.AvoidSD[4]) ## create a table of results for species groups Basic.Band.Avoid.Group.Table = data.frame(Species = group, Avoidance = Basic.Band.Group.Avoid, AvoidanceSD = Basic.Band.Group.AvoidSD, AvoidanceLCL = Basic.Band.Group.AvoidLowCI, AvoidanceUCL = Basic.Band.Group.AvoidUpCI) ## Add data for all gulls and terns BBAvoid <- 1-(sum(coll.dat$Estimated.Collisions)/sum(coll.dat$PredColl.NoAvoid.Opt1)) BBAvoidSD <- std.dev(coll.dat$Estimated.Collisions, coll.dat$PredColl.NoAvoid.Opt1) BBAvoidUpCI <- BBAvoid + (1.96 *BBAvoidSD) BBAvoidLowCI <- BBAvoid - (1.96 * BBAvoidSD) Basic.Band.Avoid.Group.Table$Species <- as.character(Basic.Band.Avoid.Group.Table$Species) Basic.Band.Avoid.Group.Table <- rbind(Basic.Band.Avoid.Group.Table, c("Gulls & terns", BBAvoid, BBAvoidSD, BBAvoidUpCI, BBAvoidLowCI)) Basic.Band.Avoid.Group.Table = rbind(Basic.Band.Avoid.Spp.Table,Basic.Band.Avoid.Group.Table) Basic.Band.Avoid.Group.Table$Avoidance <- as.numeric(Basic.Band.Avoid.Group.Table$Avoidance) Basic.Band.Avoid.Group.Table$AvoidanceSD <- as.numeric(Basic.Band.Avoid.Group.Table$AvoidanceSD) Basic.Band.Avoid.Group.Table$AvoidanceLCL <- as.numeric(Basic.Band.Avoid.Group.Table$AvoidanceLCL) Basic.Band.Avoid.Group.Table$AvoidanceUCL <- as.numeric(Basic.Band.Avoid.Group.Table$AvoidanceUCL) Basic.Band.Avoid.Group.Table$ARs = paste(round(Basic.Band.Avoid.Group.Table[,2],4), " (",round(Basic.Band.Avoid.Group.Table[,3],4), "; ", round(Basic.Band.Avoid.Group.Table[,4],4), " - ", round(Basic.Band.Avoid.Group.Table[,5],4),")", sep = "") ## sometimes keep NAs in so going to remove thee row is species is NA Basic.Band.Avoid.Group.Table <- filter(Basic.Band.Avoid.Group.Table, Species %in% c("Gulls & terns", "tern", "little tern", "common tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull", "common gull" ,"kittiwake")) #### Now plot the avoidance rates and their CIs ## filter out the species I want for this report BasicBand_group_avoidance <- Basic.Band.Avoid.Group.Table BasicBand_group_avoidance$Species <- as.character(BasicBand_group_avoidance$Species) BasicBand_group_avoidance$Species <- as.factor(BasicBand_group_avoidance$Species) ## order the factor so the group plot in the order I want BasicBand_group_avoidance$Species <- factor(BasicBand_group_avoidance$Species, levels = c("Gulls & terns", "tern", "little tern", "common tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull", "common gull" ,"kittiwake")) ## This bit of code was to save the outputs for the Cook (2021) data sets #BasicBand_group_avoidance$Source <- "Cook (2021)" #write_csv(BasicBand_group_avoidance, file = "~/Turbine CRM analysis/Cook (2021) outputs/BasicBand_ARs_Cook.csv") ``` ```{r Group-level-avoidance-basic-band-plotting, echo=FALSE, message=FALSE, warning=FALSE} ## Now read in the Cook 2021 data Cook_BasicBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/BasicBand_ARs_Cook.csv") ## Bind together my data set and the cook data set BasicBand_group_avoidance$Source <- "Ozsanlav-Harris et al (2022)" AllBasicBand_ARs <- rbind(BasicBand_group_avoidance, Cook_BasicBandAR) ## Filter out the species I want and then alter the species names AllBasicBand_ARs <- filter(AllBasicBand_ARs, Species %in% c("Gulls & terns", "tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull" ,"kittiwake")) AllBasicBand_ARs$Species <- as.character(AllBasicBand_ARs$Species) AllBasicBand_ARs$Species <- paste(toupper(substr(AllBasicBand_ARs$Species, 1, 1)), substr(AllBasicBand_ARs$Species, 2, nchar(AllBasicBand_ARs$Species)), sep="") AllBasicBand_ARs$Species <- factor(AllBasicBand_ARs$Species, levels = c("Gulls & terns", "Tern", "Sandwich tern", "Small gull", "Large gull", "Gull", "Great black-backed gull", "Lesser black-backed gull", "Herring gull", "Black-headed gull" ,"Kittiwake")) pd<- position_dodge(0.7) ggplot(AllBasicBand_ARs, aes(x= Species, y = Avoidance, colour = Source))+ geom_point(aes(color = Source), position = pd)+ geom_errorbar(aes(ymin=AvoidanceLCL, ymax=AvoidanceUCL, colour= Source), position=pd, width=0.5) + coord_flip() + theme_bw() + scale_colour_manual(values=c("#0072B2", "#D55E00")) + ylab("Basic Band avoidance rate") + xlab("Species/Species group") + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), axis.text=element_text(size=10), axis.title=element_text(size=14, face = "bold"), plot.title = element_text(size=12, face="bold"), panel.grid.minor.x = element_blank()) ## Read in or save the plots from Cook 2021 setwd("~/Turbine CRM analysis/OH et al (2022) plots") ## Save a plot #ggsave("BasicBand_SpeciesGroup.png", # width = 28, height = 15, units = "cm") ## Now read in the Cook 2021 data Cook_BasicBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/BasicBand_ARs_Cook.csv") Cook_BasicBandAR$ARs = paste(round(Cook_BasicBandAR$Avoidance,4), " (",round(Cook_BasicBandAR$AvoidanceSD,4), "; ", round(Cook_BasicBandAR$AvoidanceLCL,4), " - ", round(Cook_BasicBandAR$AvoidanceUCL,4),")", sep = "") ## Create a table for the ARs Cook_BasicBandAR$AR_2020 <- BasicBand_group_avoidance$ARs Cook_BasicBandAR$dif <- ((BasicBand_group_avoidance$Avoidance-Cook_BasicBandAR$Avoidance)/Cook_BasicBandAR$Avoidance)*100 colnames(Cook_BasicBandAR) TableBasicARs <- subset(Cook_BasicBandAR, select = c("Species", "ARs", "AR_2020", "dif")) TableBasicARs <- filter(TableBasicARs, Species %in% c("Gulls & terns", "tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull" ,"kittiwake")) data.table::setnames(TableBasicARs, old = c("Species", "ARs", "AR_2020", "dif"), new = c("Species/Species group", "Cook (2021)", "Ozsanlav-Harris et al (2022)", "% change in rate")) TableBasicARs$`Species/Species group` <- as.character(TableBasicARs$`Species/Species group`) TableBasicARs$`Species/Species group` <- paste(toupper(substr(TableBasicARs$`Species/Species group`, 1, 1)), substr(TableBasicARs$`Species/Species group`, 2, nchar(TableBasicARs$`Species/Species group`)), sep="") TableBasicARs$`Species/Species group` <- factor(TableBasicARs$`Species/Species group`, levels = c("Gulls & terns", "Tern", "Sandwich tern", "Small gull", "Large gull", "Gull", "Great black-backed gull", "Lesser black-backed gull", "Herring gull", "Black-headed gull" ,"Kittiwake")) ## Create the table suing Formattable TableBasicARs <- TableBasicARs[c(5,1,3,2,4,10,8,7,6,9,11),] TableBasicARs$`% change in rate` <- round(TableBasicARs$`% change in rate`, digits = 3) customGreen = "#71CA97" customRed = "#ff7f7f" improvement_formatter <- formatter("span", style = x ~ style(font.weight = "bold", color = ifelse(x > 0, customGreen, ifelse(x < 0, customRed, "black"))), x ~ icontext(ifelse(x>0, "arrow-up", ifelse(x < 0, "arrow-down", "equal")), x) ) BasicTable <- formattable(TableBasicARs, align =c("l","c","c", "c"), list( `Species/Species group` = formatter("span", style = ~ style(color = "grey", font.weight = "bold")), `% change in rate`= improvement_formatter )) w <- as.htmlwidget(BasicTable, width = 1250, height = 500) path <- html_print(w, background = "white", viewer = NULL) url <- paste0("file:///", gsub("\\\\", "/", normalizePath(path))) setwd("~/Turbine CRM analysis/OH et al (2022) plots") outputpath <- paste0("BasicBandAR_", " table", ".png") #webshot::webshot(url = url, file = outputpath, # vwidth = 1250, vheight = 500) ``` ------ ## 6. Calculate predicted number of collsions for the Extended Band model This section of code calculates the number of collisions in the absence of avoidance for the extended band model It calculates the proportion of the total flux at each point along the rotor diameter using the flight height distributions in Jonstone *et al* (2014) ```{r Extended-band-predicted-collions} ## I think that this approach only used modeled flight height data in Johnston et al 2014 ## try paralleling to improve speed comb <- function(...) { mapply('cbind', ..., SIMPLIFY=FALSE) } if(nrow(coll.dat) >= detectCores()) cl <- makeCluster(detectCores()) if(nrow(coll.dat) < detectCores()) cl = makeCluster(niter) registerDoParallel(cl) ## HD.y = round(seq(-1,1,0.05),2) rad = round(seq(0,1,0.05),2) ###### very convoluted, but necessary due to the way R treats seq... c = c(0.69,0.73,0.79,0.88,0.96,1,0.98,0.92,0.85,0.8,0.75,0.7,0.64,0.58,0.52,0.47,0.41,0.37,0.3,0.24,0) coverC = data.frame(cbind(rad,c)) coll.int.out = foreach(i = 1:nrow(coll.dat), .combine = "comb") %dopar% { ## get modelled flight height distributions for each species FH.dat = flightheight[ , grepl( paste(coll.dat$Species.for.Flight.Height.Distribution[i], ".est", sep = ""), names( flightheight ) ) ] ## Calcualte probabiltiy of collsion at 41 points along the rotor blade, number of points determined by values in HD.y height = coll.dat$hub.height[i] + (HD.y * (coll.dat$rotor.diameter[i]/2)) HD.d.y = (((height-floor(height)) * (FH.dat[ceiling(height)+1] - FH.dat[floor(height)+1])) + FH.dat[floor(height)+1]) * coll.dat$rotor.diameter[i]/2 FluxMin = HD.d.y[1] * (2 * ((1 - HD.y[1] * HD.y[1])^0.5)) / 2 ### risk at lowest point on rotor blade FluxInt = FluxMin + HD.d.y[41] * (2 * ((1 - HD.y[41] * HD.y[41])^0.5)) / 2 ### risk at highest point on rotor blade for (zz in 2:40) { FluxInt = FluxInt + HD.d.y[zz] * (2 * ((1 - HD.y[zz] * HD.y[zz])^0.5)) } #### Fill in intermediate heights FluxInt = FluxInt * 0.05 * (2/pi) ## now calculate the probability of collision if a bird passes through the rotor sweep zone using the collision integral out =coll.int( HD.y = HD.y, HD.d.y = HD.d.y, FluxInt = FluxInt, Blades = 3, HubHeight = coll.dat$hub.height[i], Rotor.Diameter = coll.dat$rotor.diameter[i], Blade.Width = coll.dat$blade.width[i], Rotor.Speed = coll.dat$rotor.speed[i], Pitch = coll.dat$rotor.pitch[i], speed = coll.dat$flight.speed[i], flapglide = "flap", birdlength = coll.dat$BirdLength[i], wingspan = coll.dat$wingspan[i], FH.dat = FH.dat, Prop_Upwind = 0.5, coverC = coverC) list(out) } stopCluster(cl) coll.dat$coll.int = coll.int.out[[1]][1,] ## calculate the total flux through the rotor sweep zone and immediately above and below the rotor sweep zone coll.dat$flux.All = coll.dat$total.flux * (coll.dat$total.rotor.frontal.area/coll.dat$area.survey.window) ## calculate the predicted number of collisions coll.dat$PredColl.NoAvoid.Opt3 = coll.dat$flux.All * coll.dat$coll.int ``` ## 7. Summarise Extended band model avoidance rates for species and species groups These plots show the avoidance rates at the species level and some key groupings for the extended band model The confidence intervals here just represent the variation between data set rows in avoidance rates (a row can be one wind farm over multiple years, one wind farm over one season or 1 turbine within a wind farm) ```{r Species-level-avoidance-extended-band, message=FALSE, warning=FALSE, include=FALSE} ## Now calculate the avoidance rates for each unique value in the species column Extended.Band.Spp.Avoid = vector() Extended.Band.Spp.AvoidSD = vector() Extended.Band.Spp.AvoidUpCI = vector() Extended.Band.Spp.AvoidLowCI = vector() for(i in 1:length(spp)) { Extended.Band.Spp.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == spp[i]])) Extended.Band.Spp.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$Species == spp[i]], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$Species == spp[i]]) Extended.Band.Spp.AvoidUpCI[i] = Extended.Band.Spp.Avoid[i] + (1.96 * Extended.Band.Spp.AvoidSD[i]) Extended.Band.Spp.AvoidLowCI[i] = Extended.Band.Spp.Avoid[i] - (1.96 * Extended.Band.Spp.AvoidSD[i]) } Extended.Band.Avoid.Spp.Table = data.frame(Species = spp, Avoidance = Extended.Band.Spp.Avoid, AvoidanceSD = Extended.Band.Spp.AvoidSD, AvoidanceLCL = Extended.Band.Spp.AvoidLowCI, AvoidanceUCL = Extended.Band.Spp.AvoidUpCI) ``` ```{r Group-level-avoidance-extended-band, echo=FALSE, message=FALSE, warning=FALSE} ## Now calculate the avoidance rates for each unique value in the species group column Extended.Band.Group.Avoid = vector() Extended.Band.Group.AvoidSD = vector() Extended.Band.Group.AvoidUpCI = vector() Extended.Band.Group.AvoidLowCI = vector() for(i in 1:length(group)) { Extended.Band.Group.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == group[i]])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group == group[i]])) Extended.Band.Group.AvoidSD[i] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group == group[i]], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group == group[i]]) Extended.Band.Group.AvoidUpCI[i] = Extended.Band.Group.Avoid[i] + (1.96 * Extended.Band.Group.AvoidSD[i]) Extended.Band.Group.AvoidLowCI[i] = Extended.Band.Group.Avoid[i] - (1.96 * Extended.Band.Group.AvoidSD[i]) } Extended.Band.Group.Avoid[4] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")])/sum(coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group %in% c("small gull", "large gull", "gull")])) Extended.Band.Group.AvoidSD[4] = std.dev(coll.dat$Estimated.Collisions[coll.dat$group %in% c("small gull", "large gull", "gull")], coll.dat$PredColl.NoAvoid.Opt3[coll.dat$group %in% c("small gull", "large gull", "gull")]) Extended.Band.Group.AvoidUpCI[4] = Extended.Band.Group.Avoid[4] + (1.96 * Extended.Band.Group.AvoidSD[4]) Extended.Band.Group.AvoidLowCI[4] = Extended.Band.Group.Avoid[4] - (1.96 * Extended.Band.Group.AvoidSD[4]) Extended.Band.Avoid.Group.Table = data.frame(Species = group, Avoidance = Extended.Band.Group.Avoid, AvoidanceSD = Extended.Band.Group.AvoidSD, AvoidanceLCL = Extended.Band.Group.AvoidLowCI, AvoidanceUCL = Extended.Band.Group.AvoidUpCI) ## Add data for all gulls and terns ExtBAvoid <- 1-(sum(coll.dat$Estimated.Collisions)/sum(coll.dat$PredColl.NoAvoid.Opt3)) ExtBAvoidSD <- std.dev(coll.dat$Estimated.Collisions, coll.dat$PredColl.NoAvoid.Opt3) ExtBAvoidUpCI <- ExtBAvoid + (1.96 *ExtBAvoidSD) ExtBAvoidLowCI <- ExtBAvoid - (1.96 * ExtBAvoidSD) Extended.Band.Avoid.Group.Table$Species <- as.character(Extended.Band.Avoid.Group.Table$Species) Extended.Band.Avoid.Group.Table <- rbind(Extended.Band.Avoid.Group.Table, c("Gulls & terns", ExtBAvoid, ExtBAvoidSD, ExtBAvoidUpCI, ExtBAvoidLowCI)) Extended.Band.Avoid.Group.Table = rbind(Extended.Band.Avoid.Spp.Table, Extended.Band.Avoid.Group.Table) Extended.Band.Avoid.Group.Table$Avoidance <- as.numeric(Extended.Band.Avoid.Group.Table$Avoidance) Extended.Band.Avoid.Group.Table$AvoidanceSD <- as.numeric(Extended.Band.Avoid.Group.Table$AvoidanceSD) Extended.Band.Avoid.Group.Table$AvoidanceLCL <- as.numeric(Extended.Band.Avoid.Group.Table$AvoidanceLCL) Extended.Band.Avoid.Group.Table$AvoidanceUCL <- as.numeric(Extended.Band.Avoid.Group.Table$AvoidanceUCL) Extended.Band.Avoid.Group.Table$ARs = paste(round(Extended.Band.Avoid.Group.Table[,2],4), " (",round(Extended.Band.Avoid.Group.Table[,3],4), "; ", round(Extended.Band.Avoid.Group.Table[,4],4), " - ", round(Extended.Band.Avoid.Group.Table[,5],4),")", sep = "") ## sometimes keep NAs in so going to remove thee row is species is NA Extended.Band.Avoid.Group.Table <- filter(Extended.Band.Avoid.Group.Table, Species %in% c("Gulls & terns", "tern", "little tern", "common tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull", "common gull" ,"kittiwake")) ## Now plot the avoidance rates and their CIs ## filter out the species I want for this report ExtBand_group_avoidance <- Extended.Band.Avoid.Group.Table ExtBand_group_avoidance$Species <- as.character(ExtBand_group_avoidance$Species) ExtBand_group_avoidance$Species <- as.factor(ExtBand_group_avoidance$Species) ## order the factgor so the group plot in the order I want ExtBand_group_avoidance$Species <- factor(ExtBand_group_avoidance$Species, levels = c("Gulls & terns", "tern", "little tern", "common tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull", "common gull" ,"kittiwake")) ## Save the data set for this plot #ExtBand_group_avoidance$Source <- "Cook (2021)" #write_csv(ExtBand_group_avoidance, file = "~/Turbine CRM analysis/Cook (2021) outputs/ExtBand_ARs_Cook.csv") ## Read in the data set Cook_ExtBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/ExtBand_ARs_Cook.csv") Cook_ExtBandAR$ARs = paste(round(Cook_ExtBandAR$Avoidance,4), " (",round(Cook_ExtBandAR$AvoidanceSD,4), "; ", round(Cook_ExtBandAR$AvoidanceLCL,4), " - ", round(Cook_ExtBandAR$AvoidanceUCL,4),")", sep = "") ``` ```{r Group-level-avoidance-extended-band-plotting, echo=FALSE, message=FALSE, warning=FALSE} ## Bind together my data set and the cook data set ExtBand_group_avoidance$Source <- "Ozsanlav-Harris et al (2022)" AllExtBand_ARs <- rbind(ExtBand_group_avoidance, Cook_ExtBandAR) ## Filter out the species I want and then alter the species names AllExtBand_ARs <- filter(AllExtBand_ARs, Species %in% c("Gulls & terns", "tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull" ,"kittiwake")) AllExtBand_ARs$Species <- as.character(AllExtBand_ARs$Species) AllExtBand_ARs$Species <- paste(toupper(substr(AllExtBand_ARs$Species, 1, 1)), substr(AllExtBand_ARs$Species, 2, nchar(AllExtBand_ARs$Species)), sep="") AllExtBand_ARs$Species <- factor(AllExtBand_ARs$Species, levels = c("Gulls & terns", "Tern", "Sandwich tern", "Small gull", "Large gull", "Gull", "Great black-backed gull", "Lesser black-backed gull", "Herring gull", "Black-headed gull" ,"Kittiwake")) pd<- position_dodge(0.7) ggplot(AllExtBand_ARs, aes(x= Species, y = Avoidance, colour = Source))+ geom_point(aes(color = Source), position = pd)+ geom_errorbar(aes(ymin=AvoidanceLCL, ymax=AvoidanceUCL, colour= Source), position=pd, width=0.5) + coord_flip() + theme_bw() + scale_colour_manual(values=c("#0072B2", "#D55E00")) + ylab("Extended Band avoidance rate") + xlab("Species/Species group") + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), axis.text=element_text(size=10), axis.title=element_text(size=14, face = "bold"), plot.title = element_text(size=12, face="bold"), panel.grid.minor.x = element_blank()) ## Read in or save the plots from Cook 2021 setwd("~/Turbine CRM analysis/OH et al (2022) plots") ## Save a plot #ggsave("ExtBand_SpeciesGroup.png", # width = 28, height = 15, units = "cm") ## Now read in the Cook 2021 data Cook_ExtBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/ExtBand_ARs_Cook.csv") Cook_ExtBandAR$ARs = paste(round(Cook_ExtBandAR$Avoidance,4), " (",round(Cook_ExtBandAR$AvoidanceSD,4), "; ", round(Cook_ExtBandAR$AvoidanceLCL,4), " - ", round(Cook_ExtBandAR$AvoidanceUCL,4),")", sep = "") ## Create a table for the ARs Cook_ExtBandAR$AR_2020 <- ExtBand_group_avoidance$ARs Cook_ExtBandAR$dif <- ((ExtBand_group_avoidance$Avoidance-Cook_ExtBandAR$Avoidance)/Cook_ExtBandAR$Avoidance)*100 colnames(Cook_ExtBandAR) TableExtARs <- subset(Cook_ExtBandAR, select = c("Species", "ARs", "AR_2020", "dif")) TableExtARs <- filter(TableExtARs, Species %in% c("Gulls & terns", "tern", "sandwich tern", "small gull", "large gull", "gull", "great black-backed gull", "lesser black-backed gull", "herring gull", "black-headed gull" ,"kittiwake")) data.table::setnames(TableExtARs, old = c("Species", "ARs", "AR_2020", "dif"), new = c("Species/Species group", "Cook (2021)", "Ozsanlav-Harris et al (2022)", "% change in rate")) TableExtARs$`Species/Species group` <- as.character(TableExtARs$`Species/Species group`) TableExtARs$`Species/Species group` <- paste(toupper(substr(TableExtARs$`Species/Species group`, 1, 1)), substr(TableExtARs$`Species/Species group`, 2, nchar(TableExtARs$`Species/Species group`)), sep="") TableExtARs$`Species/Species group` <- factor(TableExtARs$`Species/Species group`, levels = c("Gulls & terns", "Tern", "Sandwich tern", "Small gull", "Large gull", "Gull", "Great black-backed gull", "Lesser black-backed gull", "Herring gull", "Black-headed gull" ,"Kittiwake")) ## Create the table suing Formattable TableExtARs <- TableExtARs[c(5,1,3,2,4,10,8,7,6,9,11),] TableExtARs$`% change in rate` <- round(TableExtARs$`% change in rate`, digits = 3) customGreen = "#71CA97" customRed = "#ff7f7f" improvement_formatter <- formatter("span", style = x ~ style(font.weight = "bold", color = ifelse(x > 0, customGreen, ifelse(x < 0, customRed, "black"))), x ~ icontext(ifelse(x>0, "arrow-up", ifelse(x < 0, "arrow-down", "equal")), x) ) ExtTable <- formattable(TableExtARs, align =c("l","c","c", "c"), list( `Species/Species group` = formatter("span", style = ~ style(color = "grey", font.weight = "bold")), `% change in rate`= improvement_formatter )) w <- as.htmlwidget(ExtTable, width = 1250, height = 500) path <- html_print(w, background = "white", viewer = NULL) url <- paste0("file:///", gsub("\\\\", "/", normalizePath(path))) setwd("~/Turbine CRM analysis/OH et al (2022) plots") outputpath <- paste0("ExtBandAR_", " table", ".png") #webshot::webshot(url = url, file = outputpath, # vwidth = 1250, vheight = 500) ``` ------ ## 8. Calculate predicted number of collisions for the stochastic basic Band model This gives confidence intervals for the avoidance rates by allowing some of the input parameters to vary for the basic band model ```{r Stochastic-band-predicted-collions} ## Now for the stochastic avoidance rates ## First step is to simulate the number of birds passing through the rotor swept areas ## Now need to work out the number of birds passing through each wind farm during the study period nsim = 1000 flux.dat = matrix(nrow = nrow(coll.dat), ncol = nsim) ### for sites where we already have total flux (e.g. Oosterbierum) flux.dat[!is.na(coll.dat$total.flux),] = coll.dat$total.flux[!is.na(coll.dat$total.flux)] ### For sites where we have a number of birds recorded and a number of hours monitoring coll.dat$passage.rate[!is.na(coll.dat$Number.of.birds.recorded)] = coll.dat$Number.of.birds.recorded[!is.na(coll.dat$Number.of.birds.recorded)]/coll.dat$N.Hours.monitoring[!is.na(coll.dat$Number.of.birds.recorded)] for(i in 1:nrow(coll.dat[!is.na(coll.dat$passage.rate),])) { nbirds.pass = rpois(nsim,coll.dat$passage.rate[!is.na(coll.dat$passage.rate)][i]) flux.dat[!is.na(coll.dat$passage.rate),] [i,]= (nbirds.pass * coll.dat$HoursNight[!is.na(coll.dat$passage.rate)][[i]] * coll.dat$Nocturnal.Correction[!is.na(coll.dat$passage.rate)] [i]) + (nbirds.pass * coll.dat$HoursDay[!is.na(coll.dat$passage.rate)] [[i]]) } ## get total number of flights through rotor sweep NOT adjusted for proportion at risk height for(i in 1:ncol(flux.dat)) { flux.dat[,i] = flux.dat[,i] / coll.dat$area.survey.window * coll.dat$total.rotor.frontal.area } ### Do the same for sites where bird data are availabile as densities rather than counts for(i in 1:nrow(flux.dat[!is.na(coll.dat$density),])) { flux.dat[!is.na(coll.dat$density),][i,] = (3600/(coll.dat$width.of.survey.area[!is.na(coll.dat$density)][i]/rtruncnorm(n = nsim, a = 4, coll.dat$flight.speed[!is.na(coll.dat$density)][i]))) * coll.dat$HoursDay[!is.na(coll.dat$density)][i] * (rpois(nsim, coll.dat$density[!is.na(coll.dat$density)][i]) * coll.dat$Survey.area[!is.na(coll.dat$density)][i]) } ## We now have 1000 estimates of the number of birds passing through turbine rotors, need to correct for PCH and probability of collision ## Now work out the probability of collision at each site ## vectors needed for the Band model ##### now simulate Pcoll Pcoll.dat = matrix(nrow = nrow(coll.dat), ncol = nsim) for(i in 1:nrow(Pcoll.dat)) { for(y in 1:ncol(Pcoll.dat)){ simRotorSpeed = rnorm(1,coll.dat$rotor.speed[i],0.5) simRotorPitch = rnorm(1,coll.dat$rotor.pitch[i],0.1) simBodyLength = rnorm(1, coll.dat$BirdLength[i], coll.dat$BirdLengthSD[i]) simWingspan = rnorm(1, coll.dat$wingspan[i], coll.dat$wingspanSD[i]) simFlightSpeed = rtruncnorm(n = 1,a = 4, mean = coll.dat$flight.speed[i], sd = coll.dat$flightspeedSD[i]) Pcoll.dat[i,y] = pcoll(k =1, no.blades = 3, max.chord = coll.dat$blade.width[i], pitch = simRotorPitch, rotor.rad = coll.dat$rotor.diameter[i]/2, rotor.period = 60/simRotorSpeed,wingspan = simWingspan, bird.length = simBodyLength, aspect = simBodyLength/simWingspan, flap.glide = 0, bird.speed = simFlightSpeed,prop.flights.upwind = 0.5,radius,chord) }} ### now estimate collisions in the absence of avoidance coll.no.avoid = matrix(nrow = nrow(coll.dat), ncol = nsim) for( i in 1:ncol(coll.no.avoid)){ coll.no.avoid[,i] = Pcoll.dat[,i] * flux.dat[,i] * coll.dat$X..flights.at.rotor.height } ``` ## 9. Summarise stochastic basic band model avoidance rates for species and species groups ```{r Species-level-avoidance-stochastic-basic-band, echo=FALSE, message=FALSE, warning=FALSE} ## now bootstrapped avoidance rates HG.Basic.sCRM.Avoid = vector() LB.Basic.sCRM.Avoid = vector() GB.Basic.sCRM.Avoid = vector() KI.Basic.sCRM.Avoid = vector() BH.Basic.sCRM.Avoid = vector() CG.Basic.sCRM.Avoid = vector() LT.Basic.sCRM.Avoid = vector() ST.Basic.sCRM.Avoid = vector() CT.Basic.sCRM.Avoid = vector() SG.Basic.sCRM.Avoid = vector() LG.Basic.sCRM.Avoid = vector() AG.Basic.sCRM.Avoid = vector() AT.Basic.sCRM.Avoid = vector() ALL.Basic.sCRM.Avoid = vector() for (i in 1:ncol(coll.no.avoid)){ HG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "herring gull"])/sum(coll.no.avoid[coll.dat$Species == "herring gull", i])) LB.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "lesser black-backed gull"])/sum(coll.no.avoid[coll.dat$Species == "lesser black-backed gull", i])) GB.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "great black-backed gull"])/sum(coll.no.avoid[coll.dat$Species == "great black-backed gull", i])) KI.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "kittiwake"])/sum(coll.no.avoid[coll.dat$Species == "kittiwake", i])) BH.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "black-headed gull"])/sum(coll.no.avoid[coll.dat$Species == "black-headed gull", i])) CG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common gull"])/sum(coll.no.avoid[coll.dat$Species == "common gull", i])) LT.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern"])/sum(coll.no.avoid[coll.dat$Species == "little tern", i])) ST.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern"])/sum(coll.no.avoid[coll.dat$Species == "sandwich tern", i])) CT.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern"])/sum(coll.no.avoid[coll.dat$Species == "common tern", i])) SG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "small gull"])/sum(coll.no.avoid[coll.dat$group == "small gull", i])) LG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "large gull"])/sum(coll.no.avoid[coll.dat$group == "large gull", i])) AG.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("large gull", "small gull", "gull")])/sum(coll.no.avoid[coll.dat$group %in% c("gull","large gull", "small gull"), i])) AT.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern"])/sum(coll.no.avoid[coll.dat$group == "tern", i])) ALL.Basic.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions)/sum(coll.no.avoid[,i])) } Basic.sCRM.Avoid.Table = data.frame(Species = c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Common Gull", "Little Tern", "Sandwich Tern", "Common Tern", "Small Gull", "Large Gull", "All Gull", "All Tern", "Gulls & terns"), Avoidance = c(median(HG.Basic.sCRM.Avoid),median(LB.Basic.sCRM.Avoid),median(GB.Basic.sCRM.Avoid),median(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0]),median(BH.Basic.sCRM.Avoid), median(CG.Basic.sCRM.Avoid),median(LT.Basic.sCRM.Avoid), median(ST.Basic.sCRM.Avoid),median(CT.Basic.sCRM.Avoid),median(SG.Basic.sCRM.Avoid),median(LG.Basic.sCRM.Avoid),median(AG.Basic.sCRM.Avoid),median(AT.Basic.sCRM.Avoid), median(ALL.Basic.sCRM.Avoid)), AvoidanceSD = c(sd(HG.Basic.sCRM.Avoid),sd(LB.Basic.sCRM.Avoid),sd(GB.Basic.sCRM.Avoid),sd(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0]),sd(BH.Basic.sCRM.Avoid), sd(CG.Basic.sCRM.Avoid),sd(LT.Basic.sCRM.Avoid), sd(ST.Basic.sCRM.Avoid),sd(CT.Basic.sCRM.Avoid),sd(SG.Basic.sCRM.Avoid),sd(LG.Basic.sCRM.Avoid),sd(AG.Basic.sCRM.Avoid),sd(AT.Basic.sCRM.Avoid), sd(ALL.Basic.sCRM.Avoid, na.rm = T)), AvoidanceLCL = c(quantile(HG.Basic.sCRM.Avoid, 0.025),quantile(LB.Basic.sCRM.Avoid, 0.025),quantile(GB.Basic.sCRM.Avoid, 0.025),quantile(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0], 0.025),quantile(BH.Basic.sCRM.Avoid, 0.025),quantile(CG.Basic.sCRM.Avoid, 0.025),quantile(LT.Basic.sCRM.Avoid, 0.025),quantile(ST.Basic.sCRM.Avoid, 0.025),quantile(CT.Basic.sCRM.Avoid, 0.025),quantile(SG.Basic.sCRM.Avoid, 0.025),quantile(LG.Basic.sCRM.Avoid, 0.025),quantile(AG.Basic.sCRM.Avoid, 0.025),quantile(AT.Basic.sCRM.Avoid, 0.025), quantile(ALL.Basic.sCRM.Avoid, 0.025)), AvoidanceUCL = c(quantile(HG.Basic.sCRM.Avoid, 0.975),quantile(LB.Basic.sCRM.Avoid, 0.975),quantile(GB.Basic.sCRM.Avoid, 0.975),quantile(KI.Basic.sCRM.Avoid[KI.Basic.sCRM.Avoid >0], 0.975),quantile(BH.Basic.sCRM.Avoid, 0.975),quantile(CG.Basic.sCRM.Avoid, 0.975),quantile(LT.Basic.sCRM.Avoid, 0.975),quantile(ST.Basic.sCRM.Avoid, 0.975),quantile(CT.Basic.sCRM.Avoid, 0.975),quantile(SG.Basic.sCRM.Avoid, 0.975),quantile(LG.Basic.sCRM.Avoid, 0.975),quantile(AG.Basic.sCRM.Avoid, 0.975),quantile(AT.Basic.sCRM.Avoid, 0.975), quantile(ALL.Basic.sCRM.Avoid, 0.975))) Basic.sCRM.Avoid.Table$ARs = paste(round(Basic.sCRM.Avoid.Table[,2],4), " (",round(Basic.sCRM.Avoid.Table[,3],4), "; ", round(Basic.sCRM.Avoid.Table[,4],4), " - ", round(Basic.sCRM.Avoid.Table[,5],4),")", sep = "") ## Now plot the avoidance rates and their CIs ## filter out the species I want for this report BasicBandSt_group_avoidance <- Basic.sCRM.Avoid.Table BasicBandSt_group_avoidance$Species <- as.character(BasicBandSt_group_avoidance$Species) BasicBandSt_group_avoidance$Species <- as.factor(BasicBandSt_group_avoidance$Species) ## order the factgor so the group plot in the order I want BasicBandSt_group_avoidance$Species <- factor(BasicBandSt_group_avoidance$Species, levels = c("Gulls & terns", "All Tern", "Little Tern", "Common Tern", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "Great Black-Backed Gull", "Lesser Black-backed Gull", "Herring Gull", "Black-headed Gull", "Common Gull" ,"Kittiwake")) ## Save the data set for this plot #BasicBandSt_group_avoidance$Source <- "Cook (2021)" #write_csv(BasicBandSt_group_avoidance, file = "~/Turbine CRM analysis/Cook (2021) outputs/StBasicBand_ARs_Cook.csv") ## Read in the data set Cook_StBasicBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/StBasicBand_ARs_Cook.csv") Cook_StBasicBandAR$ARs = paste(round(Cook_StBasicBandAR$Avoidance,4), " (",round(Cook_StBasicBandAR$AvoidanceSD,4), "; ", round(Cook_StBasicBandAR$AvoidanceLCL,4), " - ", round(Cook_StBasicBandAR$AvoidanceUCL,4),")", sep = "") ``` ``` {r Group-level-avoidance-StBasic-band-plotting, echo=FALSE, message=FALSE, warning=FALSE} ## Bind together my data set and the cook data set BasicBandSt_group_avoidance$Source <- "Ozsanlav-Harris et al (2022)" AllSt_BasicBand_ARs <- rbind(BasicBandSt_group_avoidance, Cook_StBasicBandAR) unique(AllSt_BasicBand_ARs$Species) ## Filter out the species I want and then alter the species names AllSt_BasicBand_ARs <- filter(AllSt_BasicBand_ARs, Species %in% c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "All Tern", "Gulls & terns")) AllSt_BasicBand_ARs$Species <- as.character(AllSt_BasicBand_ARs$Species) AllSt_BasicBand_ARs$Species <- paste(toupper(substr(AllSt_BasicBand_ARs$Species, 1, 1)), substr(AllSt_BasicBand_ARs$Species, 2, nchar(AllSt_BasicBand_ARs$Species)), sep="") AllSt_BasicBand_ARs$Species <- factor(AllSt_BasicBand_ARs$Species, levels = c("Gulls & terns", "All Tern", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "Great Black-Backed Gull", "Lesser Black-backed Gull", "Herring Gull", "Black-headed Gull" ,"Kittiwake")) pd<- position_dodge(0.7) ggplot(AllSt_BasicBand_ARs, aes(x= Species, y = Avoidance, colour = Source))+ geom_point(aes(color = Source), position = pd)+ geom_errorbar(aes(ymin=AvoidanceLCL, ymax=AvoidanceUCL, colour= Source), position=pd, width=0.5) + coord_flip() + theme_bw() + scale_colour_manual(values=c("#0072B2", "#D55E00")) + ylab("Stochastic Basic Band avoidance rate") + xlab("Species/Species group") + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), axis.text=element_text(size=10), axis.title=element_text(size=14, face = "bold"), plot.title = element_text(size=12, face="bold"), panel.grid.minor.x = element_blank()) ## Read in or save the plots from Cook 2021 setwd("~/Turbine CRM analysis/OH et al (2022) plots") ## Save a plot #ggsave("BasicBandSt_SpeciesGroup.png", # width = 28, height = 15, units = "cm") ## Read in the data set Cook_StBasicBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/StBasicBand_ARs_Cook.csv") Cook_StBasicBandAR$ARs = paste(round(Cook_StBasicBandAR$Avoidance,4), " (",round(Cook_StBasicBandAR$AvoidanceSD,4), "; ", round(Cook_StBasicBandAR$AvoidanceLCL,4), " - ", round(Cook_StBasicBandAR$AvoidanceUCL,4),")", sep = "") ## Create a table for the ARs Cook_StBasicBandAR$AR_2020 <- BasicBandSt_group_avoidance$ARs Cook_StBasicBandAR$dif <- ((BasicBandSt_group_avoidance$Avoidance-Cook_StBasicBandAR$Avoidance)/Cook_StBasicBandAR$Avoidance)*100 colnames(Cook_StBasicBandAR) TableStBasicARs <- subset(Cook_StBasicBandAR, select = c("Species", "ARs", "AR_2020", "dif")) data.table::setnames(TableStBasicARs, old = c("Species", "ARs", "AR_2020", "dif"), new = c("Species/Species group", "Cook (2021)", "Ozsanlav-Harris et al (2022)", "% change in rate")) TableStBasicARs <- filter(TableStBasicARs, `Species/Species group` %in% c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "All Tern", "Gulls & terns")) ## Create the table using Formattable TableStBasicARs <- TableStBasicARs[c(4,5,1,2,3,9,8,7,6,10,11),] TableStBasicARs$`% change in rate` <- round(TableStBasicARs$`% change in rate`, digits = 3) customGreen = "#71CA97" customRed = "#ff7f7f" improvement_formatter <- formatter("span", style = x ~ style(font.weight = "bold", color = ifelse(x > 0, customGreen, ifelse(x < 0, customRed, "black"))), x ~ icontext(ifelse(x>0, "arrow-up", ifelse(x < 0, "arrow-down", "equal")), x) ) StBasicTable <- formattable(TableStBasicARs, align =c("l","c","c", "c"), list( `Species/Species group` = formatter("span", style = ~ style(color = "grey", font.weight = "bold")), `% change in rate`= improvement_formatter )) w <- as.htmlwidget(StBasicTable, width = 1250, height = 500) path <- html_print(w, background = "white", viewer = NULL) url <- paste0("file:///", gsub("\\\\", "/", normalizePath(path))) setwd("~/Turbine CRM analysis/OH et al (2022) plots") outputpath <- paste0("StBasicBandAR_", " table", ".png") #webshot::webshot(url = url, file = outputpath, # vwidth = 1250, vheight = 500) ``` ------ ## 10. Calculate predicted number of collisions for the stochastic extended Band model This gives confidence intervals for the avoidance rates by allowing some of the input parameters to vary for the extended band model ```{r Stochastic-extended-band-predicted-collions} if(nrow(coll.dat) >= detectCores()) cl <- makeCluster(detectCores()) if(nrow(coll.dat) < detectCores()) cl = makeCluster(niter) registerDoParallel(cl) coll.int.out2 = foreach(z = 1:nsim, .combine = "comb") %dopar% { out = vector() for(i in 1:nrow(coll.dat)){ fh.dist.sample = sample(seq(2,201,1),1) if(coll.dat$Species.for.Flight.Height.Distribution[i] == "BlackHeadedGull") FH.dat = BH.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "CommonGull") FH.dat = CG.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "LesserBlackBackedGull") FH.dat = LB.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "HerringGull") FH.dat = HG.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "GreatBlackBackedGull") FH.dat = GB.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "Kittiwake") FH.dat = KI.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "SandwichTern") FH.dat = NE.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "LittleGull") FH.dat = LG.boot[,fh.dist.sample] if(coll.dat$Species.for.Flight.Height.Distribution[i] == "CommonTern") FH.dat = CT.boot[,fh.dist.sample] height = coll.dat$hub.height[i] + (HD.y * (coll.dat$rotor.diameter[i]/2)) HD.d.y = (((height-floor(height)) * (FH.dat[ceiling(height)+1] - FH.dat[floor(height)+1])) + FH.dat[floor(height)+1]) * coll.dat$rotor.diameter[i]/2 FluxMin = HD.d.y[1] * (2 * ((1 - HD.y[1] * HD.y[1])^0.5)) / 2 ### risk at lowest point on rotor blade FluxInt = FluxMin + HD.d.y[41] * (2 * ((1 - HD.y[41] * HD.y[41])^0.5)) / 2 ### risk at highest point on rotor blade for (zz in 2:40) { FluxInt = FluxInt + HD.d.y[zz] * (2 * ((1 - HD.y[zz] * HD.y[zz])^0.5)) } #### Fill in intermediate heights FluxInt = FluxInt * 0.05 * (2/pi) out[i] =coll.int( HD.y = HD.y, HD.d.y = HD.d.y, FluxInt = FluxInt, Blades = 3, HubHeight = coll.dat$hub.height[i], Rotor.Diameter = coll.dat$rotor.diameter[i], Blade.Width = coll.dat$blade.width[i], Rotor.Speed = rnorm(1,coll.dat$rotor.speed[i],0.5), Pitch = rnorm(1, coll.dat$rotor.pitch[i],0.1), speed = rnorm(1, coll.dat$flight.speed[i],coll.dat$flightspeedSD[i]), flapglide = "flap", birdlength = rnorm(1, coll.dat$BirdLength[i],coll.dat$BirdLengthSD[i]), wingspan = rnorm(1, coll.dat$wingspan[i],coll.dat$wingspanSD[i]), FH.dat = FH.dat, Prop_Upwind = 0.5, coverC = coverC) } list(out) } write.csv(coll.int.out2[[1]], "CollIntSim.csv") ### now estimate collisions in the absence of avoidance coll.no.avoid.Opt3 = matrix(nrow = nrow(coll.dat), ncol = nsim) for( i in 1:ncol(coll.no.avoid.Opt3)){ coll.no.avoid.Opt3[,i] = coll.int.out2[[1]][,i] * flux.dat[,i] } ``` ## 11. Summarise stochastic extended band model avoidance rates for species and species groups ```{r Species-level-avoidance-stochastic-ext-band, echo=FALSE, message=FALSE, warning=FALSE} HG.Extended.sCRM.Avoid = vector() LB.Extended.sCRM.Avoid = vector() GB.Extended.sCRM.Avoid = vector() KI.Extended.sCRM.Avoid = vector() BH.Extended.sCRM.Avoid = vector() CG.Extended.sCRM.Avoid = vector() LT.Extended.sCRM.Avoid = vector() ST.Extended.sCRM.Avoid = vector() CT.Extended.sCRM.Avoid = vector() SG.Extended.sCRM.Avoid = vector() LG.Extended.sCRM.Avoid = vector() AG.Extended.sCRM.Avoid = vector() AT.Extended.sCRM.Avoid = vector() ALL.Extended.sCRM.Avoid = vector() for (i in 1:ncol(coll.no.avoid.Opt3)){ HG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "herring gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "herring gull", i])) LB.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "lesser black-backed gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "lesser black-backed gull", i])) GB.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "great black-backed gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "great black-backed gull", i])) KI.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "kittiwake"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "kittiwake", i])) BH.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "black-headed gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "black-headed gull", i])) CG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common gull"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "common gull", i])) LT.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "little tern"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "little tern", i])) ST.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "sandwich tern"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "sandwich tern", i])) CT.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$Species == "common tern"])/sum(coll.no.avoid.Opt3[coll.dat$Species == "common tern", i])) SG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "small gull"])/sum(coll.no.avoid.Opt3[coll.dat$group == "small gull", i])) LG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "large gull"])/sum(coll.no.avoid.Opt3[coll.dat$group == "large gull", i])) AG.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group %in% c("large gull", "small gull", "gull")])/sum(coll.no.avoid.Opt3[coll.dat$group %in% c("gull","large gull", "small gull"), i])) AT.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions[coll.dat$group == "tern"])/sum(coll.no.avoid.Opt3[coll.dat$group == "tern", i])) ALL.Extended.sCRM.Avoid[i] = 1-(sum(coll.dat$Estimated.Collisions)/sum(coll.no.avoid.Opt3[,i])) } Extended.sCRM.Avoid.Table = data.frame(Species = c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Common Gull", "Little Tern", "Sandwich Tern", "Common Tern", "Small Gull", "Large Gull", "All Gull", "All Tern", "Gulls & terns"), Avoidance = c(median(HG.Extended.sCRM.Avoid),median(LB.Extended.sCRM.Avoid),median(GB.Extended.sCRM.Avoid),median(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0]),median(BH.Extended.sCRM.Avoid), median(CG.Extended.sCRM.Avoid),median(LT.Extended.sCRM.Avoid), median(ST.Extended.sCRM.Avoid),median(CT.Extended.sCRM.Avoid),median(SG.Extended.sCRM.Avoid),median(LG.Extended.sCRM.Avoid),median(AG.Extended.sCRM.Avoid),median(AT.Extended.sCRM.Avoid), median(ALL.Extended.sCRM.Avoid)), AvoidanceSD = c(sd(HG.Extended.sCRM.Avoid),sd(LB.Extended.sCRM.Avoid),sd(GB.Extended.sCRM.Avoid),sd(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0]),sd(BH.Extended.sCRM.Avoid), sd(CG.Extended.sCRM.Avoid),sd(LT.Extended.sCRM.Avoid), sd(ST.Extended.sCRM.Avoid),sd(CT.Extended.sCRM.Avoid),sd(SG.Extended.sCRM.Avoid),sd(LG.Extended.sCRM.Avoid),sd(AG.Extended.sCRM.Avoid),sd(AT.Extended.sCRM.Avoid), sd(ALL.Extended.sCRM.Avoid)), AvoidanceLCL = c(quantile(HG.Extended.sCRM.Avoid, 0.025),quantile(LB.Extended.sCRM.Avoid, 0.025),quantile(GB.Extended.sCRM.Avoid, 0.025),quantile(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0], 0.025),quantile(BH.Extended.sCRM.Avoid, 0.025),quantile(CG.Extended.sCRM.Avoid, 0.025),quantile(LT.Extended.sCRM.Avoid, 0.025),quantile(ST.Extended.sCRM.Avoid, 0.025),quantile(CT.Extended.sCRM.Avoid, 0.025),quantile(SG.Extended.sCRM.Avoid, 0.025),quantile(LG.Extended.sCRM.Avoid, 0.025),quantile(AG.Extended.sCRM.Avoid, 0.025),quantile(AT.Extended.sCRM.Avoid, 0.025), quantile(ALL.Extended.sCRM.Avoid, 0.025)), AvoidanceUCL = c(quantile(HG.Extended.sCRM.Avoid, 0.975),quantile(LB.Extended.sCRM.Avoid, 0.975),quantile(GB.Extended.sCRM.Avoid, 0.975),quantile(KI.Extended.sCRM.Avoid[KI.Extended.sCRM.Avoid >0], 0.975),quantile(BH.Extended.sCRM.Avoid, 0.975),quantile(CG.Extended.sCRM.Avoid, 0.975),quantile(LT.Extended.sCRM.Avoid, 0.975),quantile(ST.Extended.sCRM.Avoid, 0.975),quantile(CT.Extended.sCRM.Avoid, 0.975),quantile(SG.Extended.sCRM.Avoid, 0.975),quantile(LG.Extended.sCRM.Avoid, 0.975),quantile(AG.Extended.sCRM.Avoid, 0.975),quantile(AT.Extended.sCRM.Avoid, 0.975), quantile(ALL.Extended.sCRM.Avoid, 0.975))) Extended.sCRM.Avoid.Table$ARs = paste(round(Extended.sCRM.Avoid.Table[,2],4), " (",round(Extended.sCRM.Avoid.Table[,3],4), "; ", round(Extended.sCRM.Avoid.Table[,4],4), " - ", round(Extended.sCRM.Avoid.Table[,5],4),")", sep = "") ## Now plot the avoidance rates and their CIs ## filter out the species I want for this report ExtBandSt_group_avoidance <- Extended.sCRM.Avoid.Table ExtBandSt_group_avoidance$Species <- as.character(ExtBandSt_group_avoidance$Species) ExtBandSt_group_avoidance$Species <- as.factor(ExtBandSt_group_avoidance$Species) ## Save the data set for this plot #ExtBandSt_group_avoidance$Source <- "Cook (2021)" #write_csv(ExtBandSt_group_avoidance, file = "~/Turbine CRM analysis/Cook (2021) outputs/StExtBand_ARs_Cook.csv") ## Read in the data set Cook_StExtBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/StExtBand_ARs_Cook.csv") Cook_StExtBandAR$ARs = paste(round(Cook_StExtBandAR$Avoidance,4), " (",round(Cook_StExtBandAR$AvoidanceSD,4), "; ", round(Cook_StExtBandAR$AvoidanceLCL,4), " - ", round(Cook_StExtBandAR$AvoidanceUCL,4),")", sep = "") ``` ```{r Group-level-avoidance-StExt-band-plotting, echo=FALSE, message=FALSE, warning=FALSE} ## order the factors so the group plot in the order I want Cook_StExtBandAR$Species <- factor(Cook_StExtBandAR$Species, levels = c("Gulls & terns", "All Tern", "Little Tern", "Common Tern", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "Great Black-Backed Gull", "Lesser Black-backed Gull", "Herring Gull", "Black-headed Gull", "Common Gull" ,"Kittiwake")) ExtBandSt_group_avoidance$Species <- factor(ExtBandSt_group_avoidance$Species, levels = c("Gulls & terns", "All Tern", "Little Tern", "Common Tern", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "Great Black-Backed Gull", "Lesser Black-backed Gull", "Herring Gull", "Black-headed Gull", "Common Gull" ,"Kittiwake")) ## Bind together my data set and the cook data set ExtBandSt_group_avoidance$Source <- "Ozsanlav-Harris et al (2022)" AllSt_ExtBand_ARs <- rbind(ExtBandSt_group_avoidance, Cook_StExtBandAR) unique(AllSt_ExtBand_ARs$Species) ## Filter out the species I want and then alter the species names AllSt_ExtBand_ARs <- filter(AllSt_ExtBand_ARs, Species %in% c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "All Tern", "Gulls & terns")) AllSt_ExtBand_ARs$Species <- as.character(AllSt_ExtBand_ARs$Species) AllSt_ExtBand_ARs$Species <- paste(toupper(substr(AllSt_ExtBand_ARs$Species, 1, 1)), substr(AllSt_ExtBand_ARs$Species, 2, nchar(AllSt_ExtBand_ARs$Species)), sep="") AllSt_ExtBand_ARs$Species <- factor(AllSt_ExtBand_ARs$Species, levels = c("Gulls & terns", "All Tern", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "Great Black-Backed Gull", "Lesser Black-backed Gull", "Herring Gull", "Black-headed Gull" ,"Kittiwake")) pd<- position_dodge(0.7) ggplot(AllSt_ExtBand_ARs, aes(x= Species, y = Avoidance, colour = Source))+ geom_point(aes(color = Source), position = pd)+ geom_errorbar(aes(ymin=AvoidanceLCL, ymax=AvoidanceUCL, colour= Source), position=pd, width=0.5) + coord_flip() + theme_bw() + ylim(0.8,1) + scale_colour_manual(values=c("#0072B2", "#D55E00")) + ylab("Stochastic Ext. Band avoidance rate") + xlab("Species/Species group") + theme(panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank(), axis.text=element_text(size=10), axis.title=element_text(size=14, face = "bold"), plot.title = element_text(size=12, face="bold"), panel.grid.minor.x = element_blank()) ## Read in or save the plots from Cook 2021 setwd("~/Turbine CRM analysis/OH et al (2022) plots") ## Save a plot #ggsave("ExtBandSt_SpeciesGroup.png", # width = 28, height = 15, units = "cm") ## Read in the data set Cook_StExtBandAR <- read_csv("~/Turbine CRM analysis/Cook (2021) outputs/StExtBand_ARs_Cook.csv") Cook_StExtBandAR$ARs = paste(round(Cook_StExtBandAR$Avoidance,4), " (",round(Cook_StExtBandAR$AvoidanceSD,4), "; ", round(Cook_StExtBandAR$AvoidanceLCL,4), " - ", round(Cook_StExtBandAR$AvoidanceUCL,4),")", sep = "") ## Create a table for the ARs Cook_StExtBandAR$AR_2020 <- ExtBandSt_group_avoidance$ARs Cook_StExtBandAR$dif <- ((ExtBandSt_group_avoidance$Avoidance-Cook_StExtBandAR$Avoidance)/Cook_StExtBandAR$Avoidance)*100 colnames(Cook_StExtBandAR) TableStExtARs <- subset(Cook_StExtBandAR, select = c("Species", "ARs", "AR_2020", "dif")) data.table::setnames(TableStExtARs, old = c("Species", "ARs", "AR_2020", "dif"), new = c("Species/Species group", "Cook (2021)", "Ozsanlav-Harris et al (2022)", "% change in rate")) TableStExtARs <- filter(TableStExtARs, `Species/Species group` %in% c("Herring Gull", "Lesser Black-backed Gull", "Great Black-Backed Gull", "Kittiwake", "Black-headed Gull", "Sandwich Tern", "Small Gull", "Large Gull", "All Gull", "All Tern", "Gulls & terns")) ## Create the table using Formattable TableStExtARs$`% change in rate` <- round(TableStExtARs$`% change in rate`, digits = 3) TableStExtARs <- TableStExtARs[c(4,5,1,2,3,9,8,7,6,10,11),] customGreen = "#71CA97" customRed = "#ff7f7f" improvement_formatter <- formatter("span", style = x ~ style(font.weight = "bold", color = ifelse(x > 0, customGreen, ifelse(x < 0, customRed, "black"))), x ~ icontext(ifelse(x>0, "arrow-up", ifelse(x < 0, "arrow-down", "equal")), x) ) StExtTable <- formattable(TableStExtARs, align =c("l","c","c", "c"), list( `Species/Species group` = formatter("span", style = ~ style(color = "grey", font.weight = "bold")), `% change in rate`= improvement_formatter )) w <- as.htmlwidget(StExtTable, width = 1250, height = 500) path <- html_print(w, background = "white", viewer = NULL) url <- paste0("file:///", gsub("\\\\", "/", normalizePath(path))) setwd("~/Turbine CRM analysis/OH et al (2022) plots") outputpath <- paste0("StExtBandAR_", " table", ".png") #webshot::webshot(url = url, file = outputpath, # vwidth = 1250, vheight = 500) ```