Simple simulation example • collision

library(collision)
library(Distance) # for distance modelling
#> Loading required package: mrds
#> This is mrds 3.0.1
#> Built: R 4.6.0; ; 2026-04-24 22:38:47 UTC; unix
#> 
#> Attaching package: 'Distance'
#> The following object is masked from 'package:mrds':
#> 
#>     create.bins
library(mc2d)  # for the PERT distribution
#> Loading required package: mvtnorm
#> 
#> Attaching package: 'mc2d'
#> The following objects are masked from 'package:base':
#> 
#>     pmax, pmin
library(ggplot2)
library(data.table) # for easy summarising and joins
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%
library(MASS) # for multivariate normal

Stochastic Example

This is a good example to start with stochastic simulation. It models a (fictional) site with four proposed turbines, of two different models, and a (fictional) point count survey of a raptor species.

For this example we assume flight density is constant across the site.

We are going to calculate a simple stochastic model output here.

In the examples below we are just outlining how to use the model but there’s a lot of analysis required beforehand to fit a distance curve and extract the effective detection radius, analyse flight heights, and understand the spatial and temporal variability in the flight patterns. Please refer to other vignettes for more information.

The inner functions of the package are designed to take single value inputs. The input functions define_bird and define_turbine can define single or stochastic inputs. So we need to:

Define stochastic rules for inputs
- bird(s): define_bird()
- turbine(s): define_turbine()
- proportion of flights at and below rotor swept height
- effective detection radius (EDR): obtained from distance modelling (Buckland et al. 2001)
- encounter rate: encounter_rate()
Sample from the input distributions into the collision calculations to determine the distribution of possible collision outcomes.

This allows the simulation to more accurately reflect the uncertainty due to both natural variance (e.g. bird body length) and measurement uncertainty (e.g. estimated flight heights).

Define inputs

Turbine model

Step one is to set up a turbine - here’s some examples of two generic turbines, set up using the define_turbine function. These parameters can be defined as probability distributions or as single numbers.


turbine_model1 <- define_turbine(
  model_id = "Turbine 180",
  blade_length = 88,
  blade_thickness_narrow = 0.5,
  blade_thickness_wide = 3,
  d_base = 8.7,
  d_rotormin = 6.57,
  d_top = 3.5,
  hh = 149,
  max_chord = 5,
  min_chord = 2,
  max_nac_h = 3,
  max_nac_l = 20,
  max_width_nacelle = 3,
  rpm = set_random("rpert", min = 0, max = 9, mode = 6.5, shape = 8),
  rotor_diam = 180,
  tilt_deg = 6,
  prop_operational = set_random("rbeta", shape1 = 20, shape2 = 0.5)
)

turbine_model2 <- define_turbine(
  model_id = "Turbine 150",
  blade_length = 73.7,
  blade_thickness_narrow = 0.4,
  blade_thickness_wide = 3,
  d_base = 7,
  d_rotormin = 5.6,
  d_top = 3,
  hh = 112,
  max_chord = 4.2,
  min_chord = 1.5,
  max_nac_h = 3,
  max_nac_l = 13,
  max_width_nacelle = 3,
  rpm = set_random("rpert", min = 0, max = 10.8, mode = 9.1, shape = 8),
  rotor_diam = 150,
  tilt_deg = 6,
  prop_operational = set_random("rbeta", shape1 = 20, shape2 = 0.5)
)

Bird species

We make a ‘bird’ object - a list of key parameters. These parameters can be defined as probability distributions or as single numbers.


wte <- define_bird(
  species = "Wedge-tailed Eagle",
  bird_length = set_random("rpert", min = 0.85, max = 1.05,
                           mean = 0.945, shape = 2) ,
  bird_speed = set_random("rlnorm", meanlog = 2.8, sdlog = 0.2),
  prop_day = set_random("runif", min = 0.48, max = 0.52),
  prop_year = 1,
  avoidance_dynamic = set_random("rpert", min = 0.88, max = 0.95, mode = 0.92,
                                 shape = 4), # shape = 4 is default
  avoidance_static = 0.9999
)

Note: the Beta PERT distribution (rpert) is a good distribution to use if you only have the mode (or mean) and the minimum and maximum parameters for a species (eg. mean body length = 0.95m, range: 0.6-1.2m). See the distributions vignette for guidance on choosing and fitting the distributions.

Define a Set of Turbines

The basic template for the model outputs is a data.frame of turbine inputs. These can be read in from a csv, or you can set them up with a few basic pieces of information and the turbine definition above.

For each turbine we need a turbine ID, and a location. Location can be NA if you are not including any spatial modelling, although you will need the average minimum distance between turbines for cluster_correction_l() to obtain an approximate correction for turbine clustering.

df_turbines <- data.frame(
  turbine_id = c("T01", "T02", "T03", "T04"),
  model = rep(c("Turbine 180", "Turbine 150"), each = 2),
  lat = c(-32.505, -32.521, -32.523, -32.516),
  lon = c(143.441, 143.442, 143.425, 143.457)
)

Survey data

Now we need to determine the input values from the observations. In order to determine the error on the estimate of the parameter we use the bootstrap method to repeatedly resample the distribution. This is just an example of how you can obtain the distribution, it is up to the analyst to determine how best to model these parameters.

The package observations dataset:

summary(df_obs)
#>     distance           size              type         height      
#>  Min.   : 136.0   Min.   :1.000   Length   :120   Min.   :  2.00  
#>  1st Qu.: 519.5   1st Qu.:1.000   N.unique :  1   1st Qu.: 33.50  
#>  Median : 818.5   Median :1.000   N.blank  :  0   Median : 77.00  
#>  Mean   : 852.6   Mean   :1.283   Min.nchar:  6   Mean   : 97.63  
#>  3rd Qu.:1134.2   3rd Qu.:2.000   Max.nchar:  6   3rd Qu.:122.00  
#>  Max.   :2334.0   Max.   :2.000                   Max.   :518.00  
#>    survey_id          object      
#>  Min.   :  2.00   Min.   :  1.00  
#>  1st Qu.: 27.25   1st Qu.: 30.75  
#>  Median : 54.50   Median : 60.50  
#>  Mean   : 53.56   Mean   : 60.50  
#>  3rd Qu.: 77.25   3rd Qu.: 90.25  
#>  Max.   :100.00   Max.   :120.00
# converting to data.table for ease of joining and summarising
dt_obs <- setDT(copy(df_obs))

is a dataset of 120 observations from 100 point count surveys, including the count (size) of raptors and associated distance at first observation (distance).

And the package survey dataset:

summary(df_survey)
#>    survey_id      survey_duration    survey_type 
#>  Min.   :  1.00   Min.   :45      Length   :100  
#>  1st Qu.: 25.75   1st Qu.:45      N.unique :  1  
#>  Median : 50.50   Median :45      N.blank  :  0  
#>  Mean   : 50.50   Mean   :45      Min.nchar:  5  
#>  3rd Qu.: 75.25   3rd Qu.:45      Max.nchar:  5  
#>  Max.   :100.00   Max.   :45

dt_survey <- setDT(copy(df_survey))

is a dataset of the metadata for the 100 point count surveys, including the duration of each survey.

It’s very important that the field data include any surveys with no sightings as well as positive detections.

Between the two of them these tables must include:

Distance to the bird at first sighting in metres. If transect sampling, this is usually the perpendicular (right-angle) distance from the transect line.
Duration of the surveys in minutes
Count of individuals in each survey (size)

Flight heights

In order to determine the flight flux we need an effective detection height, but unlike with the EDR we can’t fit a distance model to the heights because distance modelling relies on the assumption that the birds are uniformly distributed at all distances, which we know is not the case in vertical space (i.e. the density of flights tends to drop off with increasing height). The simplest way to avoid biasing the estimate by artificially inflating or deflating the flux through the turbine is to desktop truncate the observations to the maximum tip height of the turbine and use that as the effective detection height (all the observations can still be used to fit the observer’s detection function since it is just used for the horizontal distance correction). Other methods for estimating an effective detection height can be used, but this method is the most straightforward and will work well in almost all cases. Truncating the the maximum turbine height also means that the proportion of flights at rotor swept height (prop_at_height) and below rotor swept height (prop_below_height) sum to 1.

The proportion of flights at and below rotor swept height account for the amount of flights at risk of being struck by the blades of the turbine. We only need to calculate a distribution for prop_below_height since prop_at_height = 1 - prop_below_height. Since the proportion should be bounded at 0 and 1, a beta distribution is a good option.

# see distributions vignette for guidance on (one way) to fit the prop_below_height distribution
# turbine1
prop_below_height1 <- set_random("rbeta", shape1 = 6.86, shape2 = 35.3)

# turbine2
prop_below_height2 <- set_random("rbeta", shape1 = 2.65, shape2 = 42.5)

Interactions

For this example, we assume you have done the required distance correction and have arrived at a distance model.

summary(ds_raptor)
#> 
#> Summary for distance analysis 
#> Number of observations :  120 
#> Distance range         :  0  -  2334 
#> 
#> Model       : Half-normal key function with cosine adjustment terms of order 2,3 
#> 
#> Strict monotonicity constraints were enforced.
#> AIC         :  1812.536 
#> Optimisation:  mrds (slsqp) 
#> 
#> Detection function parameters
#> Scale coefficient(s):  
#>             estimate         se
#> (Intercept) 6.822439 0.08498135
#> 
#> Adjustment term coefficient(s):  
#>                 estimate        se
#> cos, order 2 -0.33444481 0.1543828
#> cos, order 3  0.09652554 0.1490168
#> 
#>                        Estimate        SE        CV
#> Average p             0.6314592  0.141781 0.2245292
#> N in covered region 190.0360410 43.949117 0.2312673

As discussed above, we are desktop truncating to the maximum rotor swept height of the turbine, so the df_obs_summary used to calculate the interaction needs to be filtered to just birds at at-risk height (see other vignettes for more detail).

# see distributions vignette for guidance on (one way) to fit the turbine flights distribution
# turbine1
n_interactions1 <- set_random("rlnorm", meanlog = -8.5, sdlog = 0.19)

# turbine2
n_interactions2 <- set_random("rlnorm", meanlog = -8.7, sdlog = 0.2)

Run the simulation

Here we set up a loop to run the simulation. Note that each iteration of the loop is independent to others, meaning that you can wrap the loop inside methods to run the iterations in parallel to save time if you wish.

Step 0 - Set up simulation parameters
Step 1 - Sample inputs for bird, turbine, prop at height, prop below height and interactions per turbine per minute.
Step 2 - Calculate turbine_flights_year() with sampled values - calculate flights through turbine per year
Step 3 - Run prob_collision_static() and prob_collision_dynamic() - calculate $P(C|I)$ for one interaction
Step 4 - Run n_collisions() - calculate number of collisions a year

For a breakdown of each step see the deterministic example.

## inputs needed:
## seed for reproducibility
## iterations - how many runs
## df_turbines  - class `data.frame`
## wte - class `birdInput`
## turbine_model - class `turbineInput`
## prop_below_height - class `randomInput`
## n_interactions - class `randomInput`


set.seed(1234)
iterations <- 1000

lst_results <- lapply(1:iterations, function(i) {
  df_turbines_results <- df_turbines

  # Step 1 - Sample inputs for bird, turbine, prop at height and prop below height
  bird_i <- sample_input(wte)
  
  turbines1_i <- sample_input(turbine_model1)
  prop_below_height1_i <- sample_input(prop_below_height1)
  prop_at_height1_i <- 1-prop_below_height1_i
  interactions1_i <- sample_input(n_interactions1)
  
  turbines2_i <- sample_input(turbine_model2)
  prop_below_height2_i <- sample_input(prop_below_height2)
  prop_at_height2_i <- 1-prop_below_height2_i
  interactions2_i <- sample_input(n_interactions2)
  
  # Step 2 calculate flux through turbine per year
  df_turbines_results[df_turbines_results$model == "Turbine 180",
                      "flight_turbine_year"] <- flights_per_year(
    flights_per_time = interactions1_i,
    time_units = "min",
    prop_day = bird_i$prop_day,
    prop_year = bird_i$prop_year # present all year
  )
  
  df_turbines_results[df_turbines_results$model == "Turbine 150",
                      "flight_turbine_year"] <- flights_per_year(
    flights_per_time = interactions2_i,
    time_units = "min",
    prop_day = bird_i$prop_day,
    prop_year = bird_i$prop_year # present all year
  )

  # Step 3 - Run `p_collisions()` - calculate $P(C|I)$ for one interaction
  df_turbines_results[df_turbines_results$model == "Turbine 180",
                      "p_coll_static"] <- prob_collision_static(
    d_base = turbines1_i$d_base,
    d_rotormin = turbines1_i$d_rotormin,
    d_top = turbines1_i$d_top,
    hh = turbines1_i$hh,
    blade_length = turbines1_i$blade_length,
    max_nac_h = turbines1_i$max_nac_h,
    max_nac_l = turbines1_i$max_nac_l,
    max_width_nacelle = turbines1_i$max_width_nacelle,
    rotor_diam = turbines1_i$rotor_diam,
    tilt_deg = turbines1_i$tilt_deg,
    max_chord = turbines1_i$max_chord,
    min_chord = turbines1_i$min_chord,
    blade_thickness_wide = turbines1_i$blade_thickness_wide,
    blade_thickness_narrow = turbines1_i$blade_thickness_narrow,
    prop_at_height = prop_at_height1_i,
    prop_below_height = prop_below_height1_i
  )
  
  df_turbines_results[df_turbines_results$model == "Turbine 150",
                      "p_coll_static"] <- prob_collision_static(
    d_base = turbines2_i$d_base,
    d_rotormin = turbines2_i$d_rotormin,
    d_top = turbines2_i$d_top,
    hh = turbines2_i$hh,
    blade_length = turbines2_i$blade_length,
    max_nac_h = turbines2_i$max_nac_h,
    max_nac_l = turbines2_i$max_nac_l,
    max_width_nacelle = turbines2_i$max_width_nacelle,
    rotor_diam = turbines2_i$rotor_diam,
    tilt_deg = turbines2_i$tilt_deg,
    max_chord = turbines2_i$max_chord,
    min_chord = turbines2_i$min_chord,
    blade_thickness_wide = turbines2_i$blade_thickness_wide,
    blade_thickness_narrow = turbines2_i$blade_thickness_narrow,
    prop_at_height = prop_at_height2_i,
    prop_below_height = prop_below_height2_i
  )

  df_turbines_results[df_turbines_results$model == "Turbine 180",
                      "p_coll_dyn"] <- prob_collision_dynamic(
    rpm = turbines1_i$rpm, # s_rot,
    blade_length = turbines1_i$blade_length,
    max_width_nacelle = turbines1_i$max_width_nacelle,
    rotor_diam = turbines1_i$rotor_diam,
    blade_thickness_wide = turbines1_i$blade_thickness_wide,
    blade_thickness_narrow = turbines1_i$blade_thickness_narrow,
    hh = turbines1_i$hh,
    bird_length = bird_i$bird_length,
    bird_speed = bird_i$bird_speed,
    prop_at_height = prop_at_height1_i,
    prop_below_height = prop_below_height1_i,
    prop_operational = turbines1_i$prop_operational
  )
  
  df_turbines_results[df_turbines_results$model == "Turbine 150",
                      "p_coll_dyn"] <- prob_collision_dynamic(
    rpm = turbines2_i$rpm, # s_rot,
    blade_length = turbines2_i$blade_length,
    max_width_nacelle = turbines2_i$max_width_nacelle,
    rotor_diam = turbines2_i$rotor_diam,
    blade_thickness_wide = turbines2_i$blade_thickness_wide,
    blade_thickness_narrow = turbines2_i$blade_thickness_narrow,
    hh = turbines2_i$hh,
    bird_length = bird_i$bird_length,
    bird_speed = bird_i$bird_speed,
    prop_at_height = prop_at_height2_i,
    prop_below_height = prop_below_height2_i,
    prop_operational = turbines2_i$prop_operational
  )


  # Step 4 - Run `n_collisions()` - calculate number of collisions a year
  df_turbines_results$n_collision <- n_collision(
    avoidance_rate_static = bird_i$avoidance_static,
    avoidance_rate_dynamic = bird_i$avoidance_dynamic,
    n_flights = df_turbines_results$flight_turbine_year,
    p_coll_static = df_turbines_results$p_coll_static,
    p_coll_dynamic = df_turbines_results$p_coll_dyn
  )

  ## for reporting - optional
  df_turbines_results$iteration <- i

  return(df_turbines_results)
})

Now we have the output of each simulation. We can use them to plot out the histogram of collision results. This shows the distribution of the long-term average annual collision rate for the wind farm.

lapply(lst_results, function(x) {
  sum(x$n_collision)
}) |>
  unlist() |>
  hist(x = _, xlab = "n_collision", main = "Histogram of collision results")

References

Buckland, S., D. Anderson, K. Burnham, Jeffrey Laake, David Borchers, and Len Thomas. 2001. Introduction to Distance Sampling: Estimating Abundance of Biological Populations. Xv. Oxford University Press. https://doi.org/10.1093/oso/9780198506492.001.0001.