Issue46a

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: collision
 Type: Package
 Title: Per-turbine Collision Risk Model for predicting avian collision with on-shore and off-shore wind turbines
-Version: 1.1.9001
+Version: 1.1.9002
 Authors@R: c(
   person("Elizabeth", "Stark", ,"estark@symbolix.com.au", role = c("aut", "cre"), 
     comment = c(ORCID = "0009-0002-6337-5021")),
@@ -22,6 +22,7 @@ Imports:
   Distance,
   methods,
   Rdpack,
+  sf,
   units
 Suggests: 
   tinytest,

NAMESPACE

@@ -4,6 +4,8 @@ S3method(sample_input,birdInput)
 S3method(sample_input,default)
 S3method(sample_input,randInput)
 S3method(sample_input,turbineInput)
+export(cluster_correction_a)
+export(cluster_correction_l)
 export(define_bird)
 export(define_turbine)
 export(edr_from_distmodel)

NEWS.md

@@ -1,7 +1,8 @@
 # collision (public beta)
 
-# collision 1.1.9000 (version 1.1 development)
+# collision 1.1 (in development)
 
+* [issue46a](https://github.com/SymbolixAU/collision/issues/46a) added spatial correction to account for turbine clustering when turbine spacing is smaller than the effective survey width/area.
 * [issue40](https://github.com/SymbolixAU/collision/issues/40) improved handling of height of observed flux window to align with the max height of the turbine.
 * [issue41](https://github.com/SymbolixAU/collision/issues/41) bug fix to ensure the proportion of time a turbine is operational is accounted for in the `p_collision_dynamic` function.
 

R/cluster_correction.R

@@ -0,0 +1,253 @@
+#' Correct for turbine clustering
+#' 
+#' The construction of the CRM requires that the flight flux through the observed 
+#' window be scaled to the size of a turbine and then applied to each turbine. If 
+#' If the (effective) detection range of the observer overlaps multiple turbines
+#' the flux needs to be scaled appropriate over the turbines in the observation
+#' range.
+#' 
+#' @details
+#' If the turbine layout is more compact and/or your effective detection width is large, then the effective detection area contains more than one turbine, in which case the observed flux needs to be weighted to account for the effective turbines in the observed area. This accounts for the fact that the observed flux is an instantaneous measure and cannot interact with multiple turbines at once.
+#' 
+#' This package provides a 1D line correction [cluster_correction_l()] and 2D aerial [cluster_correction_a()].
+#' 
+#' ##  [cluster_correction_l()]
+#' 
+#' The 1D (line) correction calculates the ratio of the effective strip width (ESW) and the mean distance between each turbine and it's nearest neighbour. It's suitable for line transects or for interim CRM early in the project life-cycle where the turbine spacing can be estimated but the layout may not be available.
+#' 
+#' ## [cluster_correction_a()]
+#' 
+#' The 2D aerial correction places a buffer of radius equal to the effective detection radius (EDR) around each turbine, then generates a spatial union of those buffers. The correction is the ratio of the actual turbine area to the area if the turbine buffers do not overlap.
+#'  
+#' ## Relationship to other spatial corrections
+#' 
+#' This correction is derived from different assumptions but fulfils a similar model correction as the \code{sqrt(n_turbines)} in previous CRM  \insertCite{Smales2013,Band2012}{collision}. 
+#' 
+#' Alternatively, a spatial probability surface can be generated by the analyst using kernel density estimate (KDE) or similar and used instead to estimate the proportion of flights at a given turbine estimation.  That approach requires more data and site coverage.
+#' 
+#'
+#' 
+#'    
+#'    
+#' @name cluster_correction    
+
+
+
+#' @rdname cluster_correction
+#'
+#' @param eff_detection_width numeric; the total effective detection width which is 
+#'     2 x effective detection radius or 2 x effective strip (half) width. 
+#' @param df_turbines data.frame; a data.frame with one row per turbine, 
+#'    containing the coordinates of each turbine as latitude, longitude pairs. 
+#'    Coordinate columns can be named (lon, long, longitude, x) and (lat, latitude, y)
+#' @param crs optional coordinate system epsg code for the coordinates.  Defaults to crs = 4326 (WGS 84 - lat/lon). 
+#' 
+#' 
+#' @return numeric; the minimum of 1 or a fractional scaling factor, calculated as the ratio of the total area covered by the turbines buffered by the EDA divided by the area if each turbine buffer were completely independent.
+#' 
+#' @examples
+#' 
+#' # generate inputs information
+#' 
+#' df_turbines <- data.frame(
+#'   turbine_id = c("T1", "T2", "T3"),
+#'   lat = c(-35.036287, -35.03900, -35.03200),
+#'   lon = c(145.947140, 145.94500, 145.94200)
+#' )
+#' 
+#' eff_detection_width <- 500 
+#' 
+#' cluster_correction_a(eff_detection_width = 500, df_turbines, crs = 4326)
+#' 
+#' @export
+cluster_correction_a <- function(eff_detection_width, df_turbines, crs = 4326){
+
+  stopifnot('eff_detection_width must be a single numeric value' = 
+              length(eff_detection_width) == 1)
+  check_num_bounds(eff_detection_width, min = 1, max = Inf)
+
+
+  sf_turbines <- df_to_projsf(df_turbines, crs) # internal helper function
+
+  ## total area with no overlap
+  area_separate <- nrow(df_turbines)* pi*(eff_detection_width*0.5)^2  #half width==radius
+
+  ## area buffered by EDA
+  sf_turb_buffer <- sf::st_make_valid(
+    sf::st_buffer(x = sf_turbines,
+                  dist = eff_detection_width*0.5,
+                  nQuadSegs = 360))
+
+  ## union
+  sf_turb_union <- sf::st_make_valid(sf::st_union(sf_turb_buffer))
+
+  ## Area
+  turb_area_buf <- sf::st_area(sf_turb_union) |> as.numeric()
+
+  ## return
+  turbine_prob <- min( turb_area_buf/area_separate, 1.0 ) 
+
+  return(turbine_prob)
+}
+
+
+#' @rdname cluster_correction
+#'
+#' @param avg_min_distance numeric; (optional) the mean distance between each turbine and its nearest neighbour. Alternatively you can input a data.frame with turbine coordinates
+#' @param eff_detection_width numeric; the effective detection (half) width or effective detection radius
+#' @param df_turbines data.frame; (optional) a data.frame with one row per turbine, 
+#'    containing the coordinates of each turbine as latitude, longitude pairs. 
+#'    Coordinate columns can be named (lon, long, longitude, x) and (lat, latitude, y)
+#' @param crs optional coordinate system epsg code for the coordinates.  Only used if df_turbines provided. Defaults to crs = 4326 (WGS 84 - lat/lon). 
+#'
+#'
+#' @return numeric; the minimum of 1 or a fractional scaling factor, calculated as the ratio of the average distance between turbine neighbours and the effective strip width
+#' 
+#' @examples
+#' 
+#' # generate input data
+#' 
+#' 
+#' df_turbines <- data.frame(
+#'   turbine_id = c("T1", "T2", "T3"),
+#'   lat = c(-35.036287, -35.03900, -35.03200),
+#'   lon = c(145.947140, 145.94500, 145.94200)
+#' )
+#'
+#' ##option 1 - provide the average distance between turbine neighbours
+#' 
+#' sf_turbines <- sf::st_as_sf(df_turbines, 
+#'   coords = c("lon", "lat"), 
+#'   crs = 4326)
+#'   
+#' idx <- sf::st_nearest_feature(sf_turbines)
+#' dist_avg <- sf::st_distance(sf_turbines, sf_turbines[idx,], by_element=TRUE) |>
+#'    mean() |>
+#'    as.numeric()
+#' 
+#' 
+#' corr1 <- cluster_correction_l(eff_detection_width = 800, 
+#'                               avg_min_distance = dist_avg)
+#'    
+#' ## option 2 - provide the data frame (might be slower)
+#'  
+#' corr2 <- cluster_correction_l(eff_detection_width = 800, 
+#'                               df_turbines = df_turbines,
+#'                               crs = 4326 # provide crs with turbines
+#'                               )
+#'  ## should be equal (or very very close)
+#'  corr1; corr2
+#' 
+#' @export
+cluster_correction_l <- function(eff_detection_width,
+                                 avg_min_distance = NULL,
+                                 df_turbines = NULL,
+                                 crs = 4326){
+
+  stopifnot('eff_detection_width must be a single numeric value' = 
+              length(eff_detection_width) == 1)
+
+
+  check_num_bounds(eff_detection_width, min = 1, max = Inf)
+
+
+  if( is.null(avg_min_distance ) & is.null(df_turbines)){
+    stop("Please enter one of avg_min_distance or df_turbines")
+  }
+
+  if( !is.null(avg_min_distance) & !is.null(df_turbines)){
+    stop("Please enter only one of avg_min_distance and df_turbines")
+  }
+
+
+  if(!is.null(avg_min_distance)){
+    stopifnot('avg_min_distance must be a single numeric value' = 
+                length(avg_min_distance) == 1)
+    check_num_bounds(avg_min_distance, min = 1, max = Inf)
+  }
+
+  if(!is.null(df_turbines)){
+
+    sf_turbines <- df_to_projsf(df_turbines, crs) # internal helper function
+
+    #index of nearest neighbour
+    idx <- sf::st_nearest_feature(sf_turbines)
+
+    #calculate average distance
+    avg_min_distance <- sf::st_distance(sf_turbines, 
+                                        sf_turbines[idx,], 
+                                        by_element=TRUE) |>
+      mean() |>
+      as.numeric()
+  }
+
+
+  return(min(avg_min_distance/eff_detection_width, 1.0)) 
+}
+
+
+# Convert df to sf
+# 
+# helper function shared by the cluster correction functions
+# 
+# @returns sf object in projected coords
+df_to_projsf <- function(df_turbines, crs){
+
+  ## checks
+  stopifnot(
+    "df_turbines must be a data.frame" =
+      is.data.frame(df_turbines)
+  )
+  stopifnot(
+    "df_turbines must have at least one row" =
+      nrow(df_turbines) >= 1
+  )
+
+  ## Grab coordinate data and check
+  lon_col <- names(df_turbines)[
+    tolower(names(df_turbines)) %in% c("lon",
+                                       "longitude",
+                                       "long",
+                                       "x",
+                                       "easting")
+  ]
+
+  lat_col <- names(df_turbines)[
+    tolower(names(df_turbines)) %in% c("lat",
+                                       "latitude",
+                                       "y",
+                                       "northing")
+  ]
+
+  stopifnot(
+    "df_turbines must contain one longitude column named lon, long, longitude, x, or easting" =
+      length(lon_col) == 1
+  )
+  stopifnot(
+    "df_turbines must contain one latitude column named lat, latitude, y, or northing" =
+      length(lat_col) == 1
+  )
+
+  sf_turbines <- sf::st_as_sf(df_turbines, 
+                              coords = c(lon_col, lat_col), 
+                              crs = crs)
+
+
+  # Convert to projected coord system for calc in metres
+  if(sf::st_is_longlat(sf_turbines)){
+
+    check_num_bounds(sf::st_coordinates(sf_turbines)[, "X"], -180, 180 )
+    check_num_bounds(sf::st_coordinates(sf_turbines)[, "Y"], -90, 90 )
+
+    projcrs <- lonlat_to_utm(sf::st_coordinates(sf_turbines))
+    sf_turbines <- sf::st_transform(sf_turbines, crs = projcrs)
+
+  }
+
+  check_num_bounds(sf::st_coordinates(sf_turbines)[, "X"], 150000, 1e6 )
+  check_num_bounds(abs(sf::st_coordinates(sf_turbines)[, "Y"]), 1.e6, 10.e6  )
+
+
+  return(sf_turbines)
+
+}

R/flux_scaling.R

@@ -11,7 +11,12 @@
 #' output flights / min).
 #' 
 #' @param obs_flux numeric; number of flights through one unit area of vertical space
-#'    per unit time as output by [obs_flux()] or similar
+#'    per unit time as output by [obs_flux()] or similar.
+#' @param spatial_correction numeric; factor to correct flux for spatial considerations. 
+#'    Either the cluster correction factor as output by [cluster_correction_a()]/ 
+#'    [cluster_correction_l()] or similar (assuming flat activity across the site) 
+#'    or a correction based on the spatial distribution of flights on the site. 
+#'    Defaults to 1.
 #' @param rotor_diameter numeric; rotor diameter. Must be in the equivalent units to 
 #'    the unit area of `obs_flux` (i.e., if the `obs_flux` is per m\eqn{^2}, `rotor_diameter` must be in m).
 #' @param hub_height numeric; hub height. Must be in the equivalent units to 
@@ -26,13 +31,17 @@
 #' @example examples/flux_example.R
 #'
 #' @export
-turbine_flights <- function(obs_flux, rotor_diameter, hub_height) {
+turbine_flights <- function(obs_flux, 
+                            spatial_correction = 1, 
+                            rotor_diameter, 
+                            hub_height) {
   check_num_bounds(obs_flux, min = 0)
+  check_num_bounds(spatial_correction, min = 0)
   check_num_bounds(rotor_diameter, min = 0)
   check_num_bounds(hub_height, min = 0)
 
   turbine_height <- hub_height + 0.5*rotor_diameter
-  res <- obs_flux * rotor_diameter * turbine_height
+  res <- obs_flux * spatial_correction * rotor_diameter * turbine_height
   return(res)
 }
 

R/flux_wrapper.R

@@ -63,6 +63,11 @@
 #' @param survey_type character; type of survey. Currently just point transects
 #'   but methods for line transects and digital areal surveys are in development.
 #' @inheritParams obs_flux
+#' @param spatial_correction numeric; factor to correct flux for spatial considerations. 
+#'    Either the cluster correction factor as output by [cluster_correction_a()]/ 
+#'    [cluster_correction_l()] or similar (assuming flat activity across the site) 
+#'    or a correction based on the spatial distribution of flights on the site. 
+#'    Defaults to 1 (no correction).
 #' @inheritParams turbine_flights
 #' @param time_units Time units of `encounter_rate`. Defaults to "min" (i.e. flights per minute).
 #' @param prop_day numeric; number between 0 and 1 representing the proportion
@@ -84,6 +89,7 @@ turbine_flights_year <- function(
     time_units = "min",
     eff_detection_width,
     eff_detection_height,
+    spatial_correction = 1,
     rotor_diameter,
     hub_height,
     prop_day,
@@ -106,6 +112,7 @@ turbine_flights_year <- function(
 
   flights_min <- turbine_flights(
     obs_flux = obs_flux,
+    spatial_correction = spatial_correction,
     rotor_diameter = rotor_diameter,
     hub_height = hub_height
   )

R/observed_flux.R

@@ -30,7 +30,8 @@
 #' @param encounter_rate numeric; number of flights observed per unit time of survey
 #'    as output by [encounter_rate()] or similar.
 #' @param eff_detection_width numeric; Allows you to manually specify the effective detection width,
-#'    which is usually 2 x effective detection radius. Must be in the same units as `eff_detection_height`.
+#'    which is usually 2 x effective detection radius or 2 x effective strip (half) width. 
+#'    Must be in the same units as `eff_detection_height`.
 #' @param eff_detection_height numeric; The effective detection height. Because the 
 #' density of flights typically decreases with increasing height, the assumption of uniform
 #' density required for traditional distance correction doesn't hold. So the observer's
@@ -84,4 +85,4 @@ obs_flux <- function(
   obs_flux <- encounter_rate / observer_vertical_area # flights/area/time
 
   return(obs_flux)
-}
+}