{stochLAB}
is a tool to run Collision Risk Models (CRMs)
for seabirds on offshore wind farms.
The {stochLAB}
package is an adaptation of the R
code developed by Masden
(2015) to incorporate variability and uncertainty in the avian
collision risk model originally developed by Band
(2012).
Code developed under {stochLAB}
substantially
re-factored and re-structured Masden’s (heavily script-based)
implementation into a user-friendly, streamlined, well documented and
easily distributed tool. Furthermore, the package lays down the code
infrastructure for easier incorporation of new functionality, e.g. extra
parameter sampling features, model expansions, etc.
In addition, previous code underpinning core calculations for the extended model has been replaced by an alternative approach, resulting in significant gains in computational speed over Masden’s code. This optimization is particularly beneficial under a stochastic context, when core calculations are called repeatedly during simulations.
For a more detailed overview type ?stochLAB
, once
installed!
You can install the released version of stochLAB from CRAN with:
install.packages("stochLAB")
You can install the development version with:
# install.packages("devtools")
::install_github("HiDef-Aerial-Surveying/stochLAB") devtools
This is a basic example of running the stochastic collision model for one seabird species and one turbine/wind-farm scenario, with fictional input parameter data.
library(stochLAB)
# ------------------------------------------------------
# Setting some of the required inputs upfront
<- data.frame(
b_dens month = month.abb,
mean = runif(12, 0.8, 1.5),
sd = runif(12, 0.2, 0.3))
# Generic FHD bootstraps for one species, from Johnson et al (2014)
<- generic_fhd_bootstraps[[1]]
fhd_boots
# wind speed vs rotation speed vs pitch
<- data.frame(
wind_rtn_ptch wind_speed = seq_len(30),
rtn_speed = 10/(30:1),
bld_pitch = c(rep(90, 4), rep(0, 8), 5:22))
# wind availability
<- data.frame(
windavb month = month.abb,
pctg = runif(12, 85, 98))
# maintenance downtime
<- data.frame(
dwntm month = month.abb,
mean = runif(12, 6, 10),
sd = rep(2, 12))
# seasons specification
<- data.frame(
seas_dt season_id = c("a", "b", "c"),
start_month = c("Jan", "May", "Oct"), end_month = c("Apr", "Sep", "Dec"))
# ----------------------------------------------------------
# Run stochastic CRM, treating rotor radius, air gap and
# blade width as fixed parameters (i.e. not stochastic)
stoch_crm(
model_options = c(1, 2, 3),
n_iter = 1000,
flt_speed_pars = data.frame(mean = 7.26, sd = 1.5),
body_lt_pars = data.frame(mean = 0.39, sd = 0.005),
wing_span_pars = data.frame(mean = 1.08, sd = 0.04),
avoid_bsc_pars = data.frame(mean = 0.99, sd = 0.001),
avoid_ext_pars = data.frame(mean = 0.96, sd = 0.002),
noct_act_pars = data.frame(mean = 0.033, sd = 0.005),
prop_crh_pars = data.frame(mean = 0.06, sd = 0.009),
bird_dens_opt = "tnorm",
bird_dens_dt = b_dens,
flight_type = "flapping",
prop_upwind = 0.5,
gen_fhd_boots = fhd_boots,
n_blades = 3,
rtr_radius_pars = data.frame(mean = 80, sd = 0), # sd = 0, rotor radius is fixed
air_gap_pars = data.frame(mean = 36, sd = 0), # sd = 0, air gap is fixed
bld_width_pars = data.frame(mean = 8, sd = 0), # sd = 0, blade width is fixed
rtn_pitch_opt = "windSpeedReltn",
windspd_pars = data.frame(mean = 7.74, sd = 3),
rtn_pitch_windspd_dt = wind_rtn_ptch,
trb_wind_avbl = windavb,
trb_downtime_pars = dwntm,
wf_n_trbs = 200,
wf_width = 15,
wf_latitude = 56.9,
tidal_offset = 2.5,
lrg_arr_corr = TRUE,
verbose = TRUE,
seed = 1234,
out_format = "summaries",
out_sampled_pars = TRUE,
out_period = "seasons",
season_specs = seas_dt,
log_file = paste0(getwd(), "scrm_example.log")
)
This is an example usage of stoch_crm()
for multiple
scenarios. The aim is two-fold:
Suggest how input parameter datasets used in the previous
implementation can be reshaped to fit stoch_crm()
’s
interface. Suggested code is also relevant in the context of multiple
scenarios applications, since the wide tabular structure of these
datasets is likely the favoured format for users to compile input
parameters under different scenarios.
Propose a functional programming framework to run
stoch_crm()
for multiple species and wind-farm/turbines
features.
Please note the example runs on fictional data.
# --------------------------------------------------------- #
# ---- Reshaping into list-column data frames ----
# --------------------------------------------------------- #
#
# Here embracing list-columns tibbles, but lists could be used instead
# --- bird features
<- bird_pars_wide_example %>%
bird_pars ::relocate(Flight, .after = dplyr::last_col()) %>%
dplyr::pivot_longer(AvoidanceBasic:Prop_CRH_ObsSD) %>%
tidyr::mutate(
dplyrpar = dplyr::if_else(grepl("SD|sd|Sd", name), "sd", "mean"),
feature = gsub("SD|sd|Sd","", name)) %>%
::select(-name) %>%
dplyr::pivot_wider(names_from = par, values_from = value) %>%
tidyr::nest(pars = c(mean, sd)) %>%
tidyr::pivot_wider(names_from = feature, values_from = pars) %>%
tidyr::add_column(prop_upwind = 0.5)
tibble
# --- bird densities: provided as mean and sd Parameters for Truncated Normal lower
# bounded at 0
<- dens_tnorm_wide_example %>%
dens_pars ::add_column(
tibbledens_opt = rep("tnorm", nrow(.)),
.after = 1) %>%
::pivot_longer(Jan:DecSD) %>%
tidyr::mutate(
dplyrpar = dplyr::if_else(grepl("SD|sd|Sd", name), "sd", "mean"),
month = gsub("SD|sd|Sd","", name)) %>%
::select(-name) %>%
dplyr::pivot_wider(names_from = par, values_from = value) %>%
tidyr::nest(mth_dens = c(month, mean, sd))
tidyr
# --- FHD data from Johnson et al (2014) for the species under analysis
<- generic_fhd_bootstraps[bird_pars$Species]
gen_fhd_boots
# --- seasons definitions (made up)
<- list(
season_dt Arctic_Tern = data.frame(
season_id = c("breeding", "feeding", "migrating"),
start_month = c("May", "Sep", "Jan"),
end_month = c("Aug", "Dec", "Apr")),
Black_headed_Gull = data.frame(
season_id = c("breeding", "feeding", "migrating"),
start_month = c("Jan", "May", "Oct"),
end_month = c("Apr", "Sep", "Dec")),
Black_legged_Kittiwake = data.frame(
season_id = c("breeding", "feeding", "migrating"),
start_month = c("Dec", "Mar", "Sep"),
end_month = c("Feb", "Aug", "Nov")))
# --- turbine parameters
## address operation parameters first
<- turb_pars_wide_example %>%
trb_opr_pars ::select(TurbineModel, JanOp:DecOpSD) %>%
dplyr::pivot_longer(JanOp:DecOpSD) %>%
tidyr::mutate(
dplyrmonth = substr(name, 1, 3),
par = dplyr::case_when(
grepl("SD|sd|Sd", name) ~ "sd",
grepl("Mean|MEAN|mean", name) ~ "mean",
TRUE ~ "pctg"
%>%
)) ::select(-name) %>%
dplyr::pivot_wider(names_from = par, values_from = value) %>%
tidyr::nest(
tidyrwind_avbl = c(month, pctg),
trb_dwntm = c(month, mean, sd))
## address turbine features and subsequently merge operation parameters
<- turb_pars_wide_example %>%
trb_pars ::select(TurbineModel:windSpeedSD ) %>%
dplyr::relocate(RotorSpeedAndPitch_SimOption, .after = 1) %>%
dplyr::pivot_longer(RotorRadius:windSpeedSD) %>%
tidyr::mutate(
dplyrpar = dplyr::if_else(grepl("SD|sd|Sd", name), "sd", "mean"),
feature = gsub("(SD|sd|Sd)|(Mean|MEAN|mean)","", name)
%>%
) ::select(-name) %>%
dplyr::pivot_wider(names_from = par, values_from = value) %>%
tidyr::nest(pars = c(mean, sd)) %>%
tidyr::pivot_wider(names_from = feature, values_from = pars) %>%
tidyr::left_join(., trb_opr_pars)
dplyr
# --- windspeed, rotation speed and blade pitch relationship
wndspd_rtn_ptch_example
# --- windfarm parameters
<- data.frame(
wf_pars wf_id = c("wf_1", "wf_2"),
n_turbs = c(200, 400),
wf_width = c(4, 10),
wf_lat = c(55.8, 55.0),
td_off = c(2.5, 2),
large_array_corr = c(FALSE, TRUE)
)
# -------------------------------------------------------------- #
# ---- Run stoch_crm() for multiple scenarios ----
# -------------------------------------------------------------- #
# --- Set up scenario combinations
<- tidyr::expand_grid(
scenarios_specs spp = bird_pars$Species,
turb_id = trb_pars$TurbineModel,
wf_id = wf_pars$wf_id) %>%
::add_column(
tibblescenario_id = paste0("scenario_", 1:nrow(.)),
.before = 1)
# --- Set up progress bar for the upcoming iterative mapping step
<- progress::progress_bar$new(
pb format = "Running Scenario: :what [:bar] :percent eta: :eta",
width = 100,
total = nrow(scenarios_specs))
# --- Map stoch_crm() to each scenario specification via purrr::pmap
<- scenarios_specs %>%
outputs ::pmap(function(scenario_id, spp, turb_id, wf_id, ...){
purrr
$tick(tokens = list(what = scenario_id))
pb
# params for current species
<- bird_pars %>%
c_spec ::filter(Species == {{spp}}) # {{}} to avoid issues with data masking
dplyr
# density for current species
<- dens_pars %>%
c_dens ::filter(Species == {{spp}})
dplyr
# params for current turbine scenario
<- trb_pars %>%
c_turb ::filter(TurbineModel == {{turb_id}})
dplyr
# params for current windfarm scenario
<- wf_pars %>%
c_wf ::filter(wf_id == {{wf_id}})
dplyr
# inputs in list-columns need to be unlisted, either via `unlist()` or
# indexing `[[1]]`
# switching off `verbose`, otherwise console will be # cramped with log
messagesstoch_crm(
model_options = c(1, 2, 3),
n_iter = 1000,
flt_speed_pars = c_spec$Flight_Speed[[1]],
body_lt_pars = c_spec$Body_Length[[1]],
wing_span_pars = c_spec$Wingspan[[1]],
avoid_bsc_pars = c_spec$AvoidanceBasic[[1]],
avoid_ext_pars = c_spec$AvoidanceExtended[[1]],
noct_act_pars = c_spec$Nocturnal_Activity[[1]],
prop_crh_pars = c_spec$Prop_CRH_Obs[[1]],
bird_dens_opt = c_dens$dens_opt,
bird_dens_dt = c_dens$mth_dens[[1]],
flight_type = c_spec$Flight,
prop_upwind = c_spec$prop_upwind,
gen_fhd_boots = gen_fhd_boots[[spp]],
n_blades = c_turb$Blades,
rtr_radius_pars = c_turb$RotorRadius[[1]],
air_gap_pars = c_turb$HubHeightAdd[[1]],
bld_width_pars = c_turb$BladeWidth[[1]],
rtn_pitch_opt = c_turb$RotorSpeedAndPitch_SimOption,
bld_pitch_pars = c_turb$Pitch[[1]],
rtn_speed_pars = c_turb$RotationSpeed[[1]],
windspd_pars = c_turb$windSpeed[[1]],
rtn_pitch_windspd_dt = wndspd_rtn_ptch_example,
trb_wind_avbl = c_turb$wind_avbl[[1]],
trb_downtime_pars = c_turb$trb_dwntm[[1]],
wf_n_trbs = c_wf$n_turbs,
wf_width = c_wf$wf_width,
wf_latitude = c_wf$wf_lat,
tidal_offset = c_wf$td_off,
lrg_arr_corr = c_wf$large_array_corr,
verbose = FALSE,
seed = 1234,
out_format = "summaries",
out_sampled_pars = FALSE,
out_period = "seasons",
season_specs = season_dt[[spp]],
log_file = NULL
)
})
# --- close progress bar
$terminate()
pb
# --- identify elements of output list
names(outputs) <- scenarios_specs$scenario_id
outputs
This is an example usage of band_crm()
. This is for a
single species and single set of turbine parameters. This replicates the
Band (2012) worksheet. The stoch_crm()
function wraps
around this function, where band_crm()
acts in essence as a
single draw of stoch_crm()
.
Please note the example runs on fictional data.
# ------------------------------------------------------
# Run with arbitrary parameter values, for illustration
# ------------------------------------------------------
# Setting a dataframe of parameters to draw from
<- data.frame(
params flight_speed = 13.1, # Flight speed in m/s
body_lt = 0.85, # Body length in m
wing_span = 1.01, # Wing span in m
flight_type = "flapping", # flapping or gliding flight
avoid_rt_basic = 0.989, # avoidance rate for option 1 and 2
avoid_rt_ext = 0.981, # extended avoidance rate for option 3 and 4
noct_activity = 0.5, # proportion of day birds are inactive
prop_crh_surv = 0.13, # proportion of birds at collision risk height (option 1 only)
prop_upwind = 0.5, # proportion of flights that are upwind
rotor_speed = 15, # rotor speed in m/s
rotor_radius = 120, # radius of turbine in m
blade_width = 5, # width of turbine blades at thickest point in m
blade_pitch = 15, # mean radius pitch in Radians
n_blades = 3, # total number of blades per turbine
hub_height = 150, # height of hub in m above HAT
n_turbines = 100, # number of turbines in the wind farm
wf_width = 52, # width across longest section of wind farm
wf_latitude = 56, # latitude of centroid of wind farm
tidal_offset = 2.5, # mean tidal offset from HAT of the wind farm
lrg_arr_corr = TRUE # apply a large array correction?
)
# Monthly bird densities
<- data.frame(
b_dens month = month.abb,
dens = runif(12, 0.8, 1.5)
)
# flight height distribution from Johnston et al
<- data.frame(
gen_fhd_dat height = Johnston_Flight_heights_SOSS$Height,
prop = Johnston_Flight_heights_SOSS$Gannet.est
)
# monthly operational time of the wind farm
<- data.frame(
turb_oper month = month.abb,
prop_oper = runif(12,0.5,0.8)
)
::band_crm(
stochLABmodel_options = c(1,2,3),
flight_speed = params$flight_speed,
body_lt = params$body_lt,
wing_span = params$wing_span,
flight_type = params$flight_type,
avoid_rt_basic = params$avoid_rt_basic,
avoid_rt_ext = params$avoid_rt_ext,
noct_activity = params$noct_activity,
prop_crh_surv = params$prop_crh_surv,
dens_month = b_dens,
prop_upwind = params$prop_upwind,
gen_fhd = gen_fhd_dat,
site_fhd = NULL, # Option 4 only
rotor_speed = params$rotor_speed,
rotor_radius = params$rotor_radius,
blade_width = params$blade_width,
blade_pitch = params$blade_pitch,
n_blades = params$n_blades,
hub_height = params$hub_height,
chord_prof = chord_prof_5MW,
n_turbines = params$n_turbines,
turb_oper_month = turb_oper,
wf_width = params$wf_width,
wf_latitude = params$wf_latitude,
tidal_offset = params$tidal_offset,
lrg_arr_corr = params$lrg_arr_corr
)