From 56599c2498845dd5f0015bb5d24c11ac938ad79b Mon Sep 17 00:00:00 2001 From: Peggy Yandell Date: Wed, 8 Apr 2026 06:56:01 +0000 Subject: [PATCH 01/10] started working on cluster correction functions --- R/cluster_correction.R | 72 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 R/cluster_correction.R diff --git a/R/cluster_correction.R b/R/cluster_correction.R new file mode 100644 index 0000000..8ce09c3 --- /dev/null +++ b/R/cluster_correction.R @@ -0,0 +1,72 @@ +#' Correct for turbine clusters +#' +#' TODO description +#' +#' @details +#' TODO +#' +#' @param df_turbines data.frame; a data.frame with one row per turbine, +#' containing the coordinates of each turbine in WGS 84 +#' @param eff_detection_width numeric; the effective detection width, +#' which is usually 2 x effective detection radius. +#' +#' @return numeric; the cluster correction factor, which is the average number +#' of turbines per EDA (effective detection area) +#' +#' @export +cluster_correction_a <- function(df_turbines, eff_detection_width){ + # df_turbines needs turbine ids and lat/lon + # lat column can be "lat", "latitude" or "y" + # lon column can be "lon", "long", "longitude" or "x" + # currently have it requiring WGS 84 but might be better to add a CRS input? + + EDA <- pi*(eff_detection_width/2)^2 + n_turbines <- nrow(df_turbines) + + lon_col <- names(df_turbines)[tolower(names(df_turbines)) %in% c("lon", + "longitude", + "long", + "x")] + + lat_col <- names(df_turbines)[tolower(names(df_turbines)) %in% c("lat", + "latitude", + "y")] + + sf_turb <- sf::st_as_sf(df_turbines, coords = c(lon_col, lat_col), crs = 4326) + + + sf_turb_buffer <- sf::st_make_valid( + sf::st_buffer(x = sf_turb, + dist = eff_detection_width/2)) + + sf_turb_union <- sf::st_union(sf_turb_buffer) + + turb_area <- units::drop_units(sf::st_area(sf_turb_union)) + + eff_turbines <- round(n_turbines*EDA/turb_area, 16) # number of turbines that flux is effectively spread across + + # might need to add a check to make sure it's >=1? + + return(eff_turbines) +} + + +#' Correct for turbine distance (offshore) +#' +#' TODO description +#' +#' @details +#' TODO +#' +#' @param avg_min_distance numeric; the average distance between each turbine and its nearest neighbour +#' @param eff_detection_width numeric; the effective detection width, +#' which is usually 2 x effective detection radius. +#' +#' @return numeric; the cluster correction factor, which is the average number +#' of turbines per ESW (effective strip width) +#' +#' @export +cluster_correction_l <- function(avg_min_distance, eff_detection_width){ + # I think this is literally it? + return(eff_detection_width/avg_min_distance) +} From 262ba03158c10739b36e8c99f981406c9acc83fb Mon Sep 17 00:00:00 2001 From: Peggy Yandell Date: Thu, 9 Apr 2026 02:19:56 +0000 Subject: [PATCH 02/10] rough outline of cluster correction functions and vignettes. Needs some polishing but is functional. --- NAMESPACE | 2 ++ R/cluster_correction.R | 38 ++++++++++++++++++------- R/flux_scaling.R | 13 +++++++-- R/flux_wrapper.R | 2 ++ R/observed_flux.R | 2 +- man/cluster_correction_a.Rd | 29 +++++++++++++++++++ man/cluster_correction_l.Rd | 24 ++++++++++++++++ man/turbine_flights.Rd | 11 +++++-- man/turbine_flights_year.Rd | 8 ++++++ vignettes/deterministic-example.Rmd | 12 ++++++-- vignettes/simple-simulation-example.Rmd | 17 +++++++++-- 11 files changed, 136 insertions(+), 22 deletions(-) create mode 100644 man/cluster_correction_a.Rd create mode 100644 man/cluster_correction_l.Rd diff --git a/NAMESPACE b/NAMESPACE index 433d95d..27bd949 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/cluster_correction.R b/R/cluster_correction.R index 8ce09c3..4cdfbf9 100644 --- a/R/cluster_correction.R +++ b/R/cluster_correction.R @@ -3,23 +3,33 @@ #' TODO description #' #' @details -#' TODO +#' If the turbine layout is more compact and/or your EDR is large, then the +#' effective detection area contains more than one turbine. +#' The observed flux needs to be weighted to account for the effective turbines +#' in the observed area. This accounts for the fact that an observed flight can +#' only interact with one turbine at any given time. TODO #' #' @param df_turbines data.frame; a data.frame with one row per turbine, #' containing the coordinates of each turbine in WGS 84 #' @param eff_detection_width numeric; the effective detection width, #' which is usually 2 x effective detection radius. #' -#' @return numeric; the cluster correction factor, which is the average number -#' of turbines per EDA (effective detection area) +#' @return numeric; the cluster correction factor - the average number +#' of turbines per EDA (effective detection area). #' #' @export cluster_correction_a <- function(df_turbines, eff_detection_width){ # df_turbines needs turbine ids and lat/lon - # lat column can be "lat", "latitude" or "y" - # lon column can be "lon", "long", "longitude" or "x" + # lat column can be "lat", "latitude" or "y" (any case) + # lon column can be "lon", "long", "longitude" or "x" (any case) # currently have it requiring WGS 84 but might be better to add a CRS input? + # TODO (low priority) I think this needs to be optimised better - it's kinda + # slow when doing it in a loop in the stochastic example + # to do it mathematically this is a good reference + # https://stackoverflow.com/questions/1667310/combined-area-of-overlapping-circles + # but I don't want to spend my day on that rn + EDA <- pi*(eff_detection_width/2)^2 n_turbines <- nrow(df_turbines) @@ -39,15 +49,18 @@ cluster_correction_a <- function(df_turbines, eff_detection_width){ sf::st_buffer(x = sf_turb, dist = eff_detection_width/2)) - sf_turb_union <- sf::st_union(sf_turb_buffer) + sf_turb_union <- sf::st_make_valid(sf::st_union(sf_turb_buffer)) turb_area <- units::drop_units(sf::st_area(sf_turb_union)) - eff_turbines <- round(n_turbines*EDA/turb_area, 16) # number of turbines that flux is effectively spread across + turbine_prob <- round(turb_area/(n_turbines*EDA), 16) # number of turbines that flux is effectively spread across + # rounding because sf sometimes has float issues - # might need to add a check to make sure it's >=1? + # may need to add a check to ensure it's always <= 1? I kind of can't think of + # a way that could happen but I think this function needs some testing to see + # if you can break it in that way - return(eff_turbines) + return(turbine_prob) } @@ -62,11 +75,14 @@ cluster_correction_a <- function(df_turbines, eff_detection_width){ #' @param eff_detection_width numeric; the effective detection width, #' which is usually 2 x effective detection radius. #' -#' @return numeric; the cluster correction factor, which is the average number -#' of turbines per ESW (effective strip width) +#' @return numeric; the (linear) cluster correction factor - the average number +#' of turbines per ESW (effective strip width). #' #' @export cluster_correction_l <- function(avg_min_distance, eff_detection_width){ # I think this is literally it? + # TODO: I would like to have the option to input df_turbines here and have it calculate the avg min distance + # but I think that requires better knowledge of R to do in a relatively speedy way + # ES suggested spatialdatatable::dtHaversine? return(eff_detection_width/avg_min_distance) } diff --git a/R/flux_scaling.R b/R/flux_scaling.R index 22abb25..3f6497c 100644 --- a/R/flux_scaling.R +++ b/R/flux_scaling.R @@ -11,7 +11,13 @@ #' 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. +#' TODO make better words #' @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 +32,14 @@ #' @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) } diff --git a/R/flux_wrapper.R b/R/flux_wrapper.R index 0fe08e0..1f6dc62 100644 --- a/R/flux_wrapper.R +++ b/R/flux_wrapper.R @@ -84,6 +84,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 +107,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 ) diff --git a/R/observed_flux.R b/R/observed_flux.R index 08cf3d5..7e36160 100644 --- a/R/observed_flux.R +++ b/R/observed_flux.R @@ -84,4 +84,4 @@ obs_flux <- function( obs_flux <- encounter_rate / observer_vertical_area # flights/area/time return(obs_flux) -} \ No newline at end of file +} diff --git a/man/cluster_correction_a.Rd b/man/cluster_correction_a.Rd new file mode 100644 index 0000000..3aa84b8 --- /dev/null +++ b/man/cluster_correction_a.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cluster_correction.R +\name{cluster_correction_a} +\alias{cluster_correction_a} +\title{Correct for turbine clusters} +\usage{ +cluster_correction_a(df_turbines, eff_detection_width) +} +\arguments{ +\item{df_turbines}{data.frame; a data.frame with one row per turbine, +containing the coordinates of each turbine in WGS 84} + +\item{eff_detection_width}{numeric; the effective detection width, +which is usually 2 x effective detection radius.} +} +\value{ +numeric; the cluster correction factor - the average number +of turbines per EDA (effective detection area). +} +\description{ +TODO description +} +\details{ +If the turbine layout is more compact and/or your EDR is large, then the +effective detection area contains more than one turbine. +The observed flux needs to be weighted to account for the effective turbines +in the observed area. This accounts for the fact that an observed flight can +only interact with one turbine at any given time. TODO +} diff --git a/man/cluster_correction_l.Rd b/man/cluster_correction_l.Rd new file mode 100644 index 0000000..a9c105c --- /dev/null +++ b/man/cluster_correction_l.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cluster_correction.R +\name{cluster_correction_l} +\alias{cluster_correction_l} +\title{Correct for turbine distance (offshore)} +\usage{ +cluster_correction_l(avg_min_distance, eff_detection_width) +} +\arguments{ +\item{avg_min_distance}{numeric; the average distance between each turbine and its nearest neighbour} + +\item{eff_detection_width}{numeric; the effective detection width, +which is usually 2 x effective detection radius.} +} +\value{ +numeric; the (linear) cluster correction factor - the average number +of turbines per ESW (effective strip width). +} +\description{ +TODO description +} +\details{ +TODO +} diff --git a/man/turbine_flights.Rd b/man/turbine_flights.Rd index ea89a56..d9eea90 100644 --- a/man/turbine_flights.Rd +++ b/man/turbine_flights.Rd @@ -4,11 +4,18 @@ \alias{turbine_flights} \title{Convert observed flight flux to number of flights through turbine plane (i.e. interactions).} \usage{ -turbine_flights(obs_flux, rotor_diameter, hub_height) +turbine_flights(obs_flux, spatial_correction = 1, rotor_diameter, hub_height) } \arguments{ \item{obs_flux}{numeric; number of flights through one unit area of vertical space -per unit time as output by \code{\link[=obs_flux]{obs_flux()}} or similar} +per unit time as output by \code{\link[=obs_flux]{obs_flux()}} or similar.} + +\item{spatial_correction}{numeric; factor to correct flux for spatial considerations. +Either the cluster correction factor as output by \code{\link[=cluster_correction_a]{cluster_correction_a()}}/ +\code{\link[=cluster_correction_l]{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. +TODO make better words} \item{rotor_diameter}{numeric; rotor diameter. Must be in the equivalent units to the unit area of \code{obs_flux} (i.e., if the \code{obs_flux} is per m\eqn{^2}, \code{rotor_diameter} must be in m).} diff --git a/man/turbine_flights_year.Rd b/man/turbine_flights_year.Rd index 028de03..8baa004 100644 --- a/man/turbine_flights_year.Rd +++ b/man/turbine_flights_year.Rd @@ -10,6 +10,7 @@ turbine_flights_year( time_units = "min", eff_detection_width, eff_detection_height, + spatial_correction = 1, rotor_diameter, hub_height, prop_day, @@ -39,6 +40,13 @@ the flux through the turbine. It is possible to use another method to estimate an effective detection height but we leave this to the analyst's judgement. Must be in the same units as \code{eff_detection_width}.} +\item{spatial_correction}{numeric; factor to correct flux for spatial considerations. +Either the cluster correction factor as output by \code{\link[=cluster_correction_a]{cluster_correction_a()}}/ +\code{\link[=cluster_correction_l]{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. +TODO make better words} + \item{rotor_diameter}{numeric; rotor diameter. Must be in the equivalent units to the unit area of \code{obs_flux} (i.e., if the \code{obs_flux} is per m\eqn{^2}, \code{rotor_diameter} must be in m).} diff --git a/vignettes/deterministic-example.Rmd b/vignettes/deterministic-example.Rmd index a05a72f..6f397ad 100644 --- a/vignettes/deterministic-example.Rmd +++ b/vignettes/deterministic-example.Rmd @@ -251,13 +251,15 @@ This is where we account for * The observer's detection model (`eff_detection_width` obtained using `edr_dist_model()`) * The size of the turbine (contained in the object made using `define_turbine()`) +* The turbine layout * The proportion of the day active (contained in the object made using `define_bird()`) * The proportion of the year active (contained in the object made using `define_bird()`) We first calculate the observed flux through a vertical airspace. This is in flights per unit time per unit area (assuming we are below max turbine height), in this example the units are minutes and metres squared, respectively. +Second we calculate the cluster correction factor (or the spatial flights pdf if you are doing a spatial model). -Second we apply this to the "turbine plane" (a rectangular "doorway" defined by the rotor diameter and maximum rotor swept height) using `turbine_flights()` to obtain the flights through the turbine plane per minute. +We then apply these to the "turbine plane" (a rectangular "doorway" defined by the rotor diameter and maximum rotor swept height) using `turbine_flights()` to obtain the flights through the turbine plane per minute. ```{r} @@ -268,10 +270,15 @@ obs_flux_min <- obs_flux( eff_detection_height = max_rsh ) +effective_flux <- cluster_correction_a( + df_turbines = df_turbines, + eff_detection_width = 2*edr_from_distmodel(ds_model)) + #flights through turbine / min turbine_flights_min <- turbine_flights( obs_flux = obs_flux_min, + spatial_correction = effective_flux, rotor_diameter = v90_single$rotor_diam, hub_height = v90_single$hh ) @@ -281,7 +288,6 @@ turbine_flights_min <- turbine_flights( The observed flight flux is `r obs_flux_min` flights per square meter per minute (below max turbine height). Scaling up to the turbine this corresponds to `r turbine_flights_min` flights through the turbine area per minute. - Here we correct it for daily and monthly variability. This can be done by hand but the package includes a helper function which can scale to a year. ```{r} @@ -296,7 +302,7 @@ turbine_flight_year <- flights_per_year( * The units for `turbine_flight_year` is flights / year (per turbine). * We expect `r turbine_flight_year` flights through the turbine area each year. -This is the expected number of interactions per turbine per year with no avoidance. +This is the expected number of interactions^["interaction" just refers to a flight through the turbine cylinder. Below the rotor swept height a bird could be a blade length from the tower and it would still count as an interaction for our purposes. This is then accounted for in the probability of collision.] per turbine per year with no avoidance. ### Step 3 - Probability of collision given interaction diff --git a/vignettes/simple-simulation-example.Rmd b/vignettes/simple-simulation-example.Rmd index 701dd6e..df51a2e 100644 --- a/vignettes/simple-simulation-example.Rmd +++ b/vignettes/simple-simulation-example.Rmd @@ -138,8 +138,8 @@ For each turbine we need a turbine ID, and a location. Location can be NA if you df_turbines <- data.frame( turbine_id = c("T01", "T02", "T03", "T04"), model = rep(c("Turbine 180", "Turbine 150"), each = 2), - lat = c(-32.512, -32.521, -32.523, -32.516), - lon = c(143.441, 143.442, 143.444, 143.447) + lat = c(-32.505, -32.521, -32.523, -32.516), + lon = c(143.441, 143.442, 143.425, 143.457) ) ``` @@ -372,6 +372,7 @@ Here we set up a loop to run the simulation. Note that each iteration of the loo * Step 0 - Set up simulation parameters * Step 1 - Sample inputs for bird, turbine, prop at height, prop below height, effective detection radius (EDR) and encounter rate. +* Step 1.5 (optional) - Calculate the `cluster_correction` or any spatial flight corrections. The cluster correction is only needed if the distance between turbines $\leq$ EDR. * 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 @@ -413,6 +414,13 @@ lst_results <- lapply(1:iterations, function(i) { prop_at_height2_i <- 1-prop_below_height2_i er2_i <- sample_input(flights_per_min2) + # Step 1.5 calculate the cluster correction + + df_turbines_results$cluster_corr <- cluster_correction_a( + df_turbines = df_turbines_results, + eff_detection_width = 2*edr_i + ) + # Step 2 calculate flux through turbine per year df_turbines_results[df_turbines_results$model == "Turbine 180", "flight_turbine_year"] <- turbine_flights_year( @@ -420,6 +428,8 @@ lst_results <- lapply(1:iterations, function(i) { encounter_rate = er1_i, eff_detection_width = 2*edr_i, eff_detection_height = max_rsh1, + spatial_correction = df_turbines_results[df_turbines_results$model == "Turbine 180", + "cluster_corr"], rotor_diameter = turbines1_i$rotor_diam, hub_height = turbines1_i$hh, prop_day = bird_i$prop_day, @@ -432,6 +442,8 @@ lst_results <- lapply(1:iterations, function(i) { encounter_rate = er2_i, eff_detection_width = 2*edr_i, eff_detection_height = max_rsh2, + spatial_correction = df_turbines_results[df_turbines_results$model == "Turbine 150", + "cluster_corr"], rotor_diameter = turbines2_i$rotor_diam, hub_height = turbines2_i$hh, prop_day = bird_i$prop_day, @@ -528,6 +540,7 @@ lst_results <- lapply(1:iterations, function(i) { }) ``` + 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. ```{r} From 12da43698d59d9a83fbeccb228c74d77db784998 Mon Sep 17 00:00:00 2001 From: sahinya17 Date: Tue, 14 Apr 2026 10:35:01 +1000 Subject: [PATCH 03/10] Added tests for cluster correction functions --- DESCRIPTION | 1 + NAMESPACE | 6 + R/cluster_correction.R | 100 ++++++++++- inst/tinytest/test_cluster_correction.R | 221 ++++++++++++++++++++++++ 4 files changed, 325 insertions(+), 3 deletions(-) create mode 100644 inst/tinytest/test_cluster_correction.R diff --git a/DESCRIPTION b/DESCRIPTION index e442c15..6721000 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,6 +22,7 @@ Imports: Distance, methods, Rdpack, + sf, units Suggests: tinytest, diff --git a/NAMESPACE b/NAMESPACE index 27bd949..7e5b3e1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,3 +26,9 @@ export(vertical_flux_plane_area) export(w_from_distmodel) import(Distance) importFrom(Rdpack,reprompt) +importFrom(sf,st_area) +importFrom(sf,st_as_sf) +importFrom(sf,st_buffer) +importFrom(sf,st_make_valid) +importFrom(sf,st_union) +importFrom(units,drop_units) diff --git a/R/cluster_correction.R b/R/cluster_correction.R index 4cdfbf9..9f78952 100644 --- a/R/cluster_correction.R +++ b/R/cluster_correction.R @@ -14,9 +14,11 @@ #' @param eff_detection_width numeric; the effective detection width, #' which is usually 2 x effective detection radius. #' -#' @return numeric; the cluster correction factor - the average number +#' @return numeric; the cluster correction factor - the average number #' of turbines per EDA (effective detection area). -#' +#' +#' @importFrom sf st_as_sf st_buffer st_make_valid st_union st_area +#' @importFrom units drop_units #' @export cluster_correction_a <- function(df_turbines, eff_detection_width){ # df_turbines needs turbine ids and lat/lon @@ -30,6 +32,15 @@ cluster_correction_a <- function(df_turbines, eff_detection_width){ # https://stackoverflow.com/questions/1667310/combined-area-of-overlapping-circles # but I don't want to spend my day on that rn + 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 + ) + EDA <- pi*(eff_detection_width/2)^2 n_turbines <- nrow(df_turbines) @@ -42,9 +53,65 @@ cluster_correction_a <- function(df_turbines, eff_detection_width){ "latitude", "y")] + stopifnot( + "df_turbines must contain a longitude column named lon, long, longitude, or x" = + length(lon_col) >= 1 + ) + stopifnot( + "df_turbines must contain a latitude column named lat, latitude, or y" = + length(lat_col) >= 1 + ) + + if (length(lon_col) > 1) { + warning( + "multiple possible longitude columns detected: ", + paste(lon_col, collapse = ", "), + ". Using the first match: ", lon_col[1] + ) + lon_col <- lon_col[1] + } + if (length(lat_col) > 1) { + warning( + "multiple possible latitude columns detected: ", + paste(lat_col, collapse = ", "), + ". Using the first match: ", lat_col[1] + ) + lat_col <- lat_col[1] + } + + stopifnot( + "NA values detected in longitude column" = + !any(is.na(df_turbines[[lon_col]])) + ) + stopifnot( + "NA values detected in latitude column" = + !any(is.na(df_turbines[[lat_col]])) + ) + + stopifnot( + "longitude values appear out of WGS 84 range (-180 to 180)" = + all(df_turbines[[lon_col]] >= -180 & df_turbines[[lon_col]] <= 180) + ) + stopifnot( + "latitude values appear out of WGS 84 range (-90 to 90)" = + all(df_turbines[[lat_col]] >= -90 & df_turbines[[lat_col]] <= 90) + ) + + stopifnot( + "eff_detection_width must be a single numeric value" = + is.numeric(eff_detection_width) && length(eff_detection_width) == 1 + ) + stopifnot( + "eff_detection_width must be greater than 0" = + eff_detection_width > 0 + ) + stopifnot( + "NA detected in eff_detection_width" = + !is.na(eff_detection_width) + ) + sf_turb <- sf::st_as_sf(df_turbines, coords = c(lon_col, lat_col), crs = 4326) - sf_turb_buffer <- sf::st_make_valid( sf::st_buffer(x = sf_turb, dist = eff_detection_width/2)) @@ -84,5 +151,32 @@ cluster_correction_l <- function(avg_min_distance, eff_detection_width){ # TODO: I would like to have the option to input df_turbines here and have it calculate the avg min distance # but I think that requires better knowledge of R to do in a relatively speedy way # ES suggested spatialdatatable::dtHaversine? + + stopifnot( + "avg_min_distance must be a single numeric value" = + is.numeric(avg_min_distance) && length(avg_min_distance) == 1 + ) + stopifnot( + "NA detected in avg_min_distance" = + !is.na(avg_min_distance) + ) + stopifnot( + "avg_min_distance must be greater than 0" = + avg_min_distance > 0 + ) + + stopifnot( + "eff_detection_width must be a single numeric value" = + is.numeric(eff_detection_width) && length(eff_detection_width) == 1 + ) + stopifnot( + "NA detected in eff_detection_width" = + !is.na(eff_detection_width) + ) + stopifnot( + "eff_detection_width must be greater than 0" = + eff_detection_width > 0 + ) + return(eff_detection_width/avg_min_distance) } diff --git a/inst/tinytest/test_cluster_correction.R b/inst/tinytest/test_cluster_correction.R new file mode 100644 index 0000000..6f95965 --- /dev/null +++ b/inst/tinytest/test_cluster_correction.R @@ -0,0 +1,221 @@ +# Tests for cluster_correction_a and cluster_correction_l + +# Helper: simple 2-turbine data.frame +two_turbines <- data.frame(lon = c(0, 0.01), lat = c(51, 51)) + +# cluster_correction_a---------------------- + +## Returns a single numeric value ------------------------------------------- +result_a <- cluster_correction_a(two_turbines, eff_detection_width = 500) +expect_true(is.numeric(result_a)) +expect_equal(length(result_a), 1L) +expect_true(is.finite(result_a)) + +## Single turbine: correction factor should be ~1 --------------------------- +# sf computes spherical area, which differs slightly from planar pi*r^2, +# so we allow a small tolerance rather than expecting exactly 1. +one_turbine <- data.frame(lon = 0, lat = 51) +expect_equal( + cluster_correction_a(one_turbine, eff_detection_width = 500), + 1, + tolerance = 0.05 +) + +## Very widely spaced turbines: factor close to 1 --------------------------- +# Turbines far apart so EDAs don't overlap; each turbine contributes one full EDA. +far_turbines <- data.frame(lon = c(0, 10), lat = c(51, 51)) +result_far <- cluster_correction_a(far_turbines, eff_detection_width = 500) +expect_equal(result_far, 1, tolerance = 0.05) + +## Overlapping EDAs: factor lower than non-overlapping case ----------------- +# Turbines very close together relative to detection width; their EDAs nearly +# fully overlap, so the correction factor should be much less than for +# non-overlapping turbines. +close_turbines <- data.frame(lon = c(0, 0.0001), lat = c(51, 51)) +result_close <- cluster_correction_a(close_turbines, eff_detection_width = 5000) +expect_true(result_close < result_far) + +## Result is positive and finite -------------------------------------------- +expect_true(result_close > 0) +expect_true(result_far > 0) + +## Accepts alternative column names ----------------------------------------- +df_alt_names <- data.frame(longitude = c(0, 0.01), latitude = c(51, 51)) +expect_equal( + cluster_correction_a(df_alt_names, 500), + result_a, + tolerance = 1e-6 +) + +df_xy <- data.frame(x = c(0, 0.01), y = c(51, 51)) +expect_equal(cluster_correction_a(df_xy, 500), result_a, tolerance = 1e-6) + +## Case-insensitive column matching ----------------------------------------- +df_upper <- data.frame(LON = c(0, 0.01), LAT = c(51, 51)) +expect_equal(cluster_correction_a(df_upper, 500), result_a, tolerance = 1e-6) + +## Warns when multiple possible lon/lat columns detected -------------------- +df_multi_lon <- data.frame(lon = c(0, 0.01), long = c(0, 0.01), lat = c(51, 51)) +expect_warning( + cluster_correction_a(df_multi_lon, 500), + 'multiple possible longitude columns detected' +) + +df_multi_lat <- data.frame( + lon = c(0, 0.01), + lat = c(51, 51), + latitude = c(51, 51) +) +expect_warning( + cluster_correction_a(df_multi_lat, 500), + 'multiple possible latitude columns detected' +) + +## Input validation: non-data.frame ----------------------------------------- +expect_error( + cluster_correction_a(list(lon = 0, lat = 51), 500), + 'df_turbines must be a data.frame' +) + +## Input validation: zero rows --------------------------------------------- +expect_error( + cluster_correction_a(data.frame(lon = numeric(0), lat = numeric(0)), 500), + 'df_turbines must have at least one row' +) + +## Input validation: missing lon column ------------------------------------ +expect_error( + cluster_correction_a(data.frame(easting = c(0, 1), lat = c(51, 51)), 500), + 'df_turbines must contain a longitude column' +) + +## Input validation: missing lat column ------------------------------------ +expect_error( + cluster_correction_a(data.frame(lon = c(0, 1), northing = c(51, 51)), 500), + 'df_turbines must contain a latitude column' +) + +## Input validation: NA coordinates ---------------------------------------- +expect_error( + cluster_correction_a(data.frame(lon = c(0, NA), lat = c(51, 51)), 500), + 'NA values detected in longitude column' +) +expect_error( + cluster_correction_a(data.frame(lon = c(0, 0.01), lat = c(51, NA)), 500), + 'NA values detected in latitude column' +) + +## Input validation: out-of-range coordinates ------------------------------ +expect_error( + cluster_correction_a(data.frame(lon = c(0, 200), lat = c(51, 51)), 500), + 'longitude values appear out of WGS 84 range' +) +expect_error( + cluster_correction_a(data.frame(lon = c(0, 0.01), lat = c(51, 100)), 500), + 'latitude values appear out of WGS 84 range' +) + +## Input validation: eff_detection_width type/length ----------------------- +# Note: a character value causes an arithmetic error before the stopifnot fires, +# because EDA = pi*(eff_detection_width/2)^2 is evaluated first. +expect_error(cluster_correction_a(two_turbines, eff_detection_width = '500')) +expect_error( + cluster_correction_a(two_turbines, eff_detection_width = c(500, 600)), + 'eff_detection_width must be a single numeric value' +) + +## Input validation: eff_detection_width bounds ---------------------------- +expect_error( + cluster_correction_a(two_turbines, eff_detection_width = 0), + 'eff_detection_width must be greater than 0' +) +expect_error( + cluster_correction_a(two_turbines, eff_detection_width = -100), + 'eff_detection_width must be greater than 0' +) + +## Input validation: NA eff_detection_width -------------------------------- +# is.numeric(NA) is FALSE, so the "single numeric value" check fires first. +expect_error( + cluster_correction_a(two_turbines, eff_detection_width = NA), + 'eff_detection_width must be a single numeric value' +) + +# cluster_correction_l---------------------------------------------------------- + +## Basic correctness: result equals eff_detection_width / avg_min_distance---- +expect_equal( + cluster_correction_l(avg_min_distance = 1000, eff_detection_width = 500), + 0.5 +) +expect_equal( + cluster_correction_l(avg_min_distance = 500, eff_detection_width = 500), + 1 +) +expect_equal( + cluster_correction_l(avg_min_distance = 200, eff_detection_width = 500), + 2.5 +) + +## Returns single numeric --------------------------------------------------- +res_l <- cluster_correction_l(1000, 500) +expect_true(is.numeric(res_l)) +expect_equal(length(res_l), 1L) + +## Input validation: avg_min_distance type/length --------------------------- +expect_error( + cluster_correction_l(avg_min_distance = '1000', eff_detection_width = 500), + 'avg_min_distance must be a single numeric value' +) +expect_error( + cluster_correction_l( + avg_min_distance = c(500, 1000), + eff_detection_width = 500 + ), + 'avg_min_distance must be a single numeric value' +) + +## Input validation: avg_min_distance bounds -------------------------------- +expect_error( + cluster_correction_l(avg_min_distance = 0, eff_detection_width = 500), + 'avg_min_distance must be greater than 0' +) +expect_error( + cluster_correction_l(avg_min_distance = -100, eff_detection_width = 500), + 'avg_min_distance must be greater than 0' +) + +## Input validation: NA avg_min_distance ------------------------------------ +expect_error( + cluster_correction_l(avg_min_distance = NA, eff_detection_width = 500), + 'avg_min_distance must be a single numeric value' +) + +## Input validation: eff_detection_width type/length ------------------------ +expect_error( + cluster_correction_l(avg_min_distance = 1000, eff_detection_width = '500'), + 'eff_detection_width must be a single numeric value' +) +expect_error( + cluster_correction_l( + avg_min_distance = 1000, + eff_detection_width = c(500, 600) + ), + 'eff_detection_width must be a single numeric value' +) + +## Input validation: eff_detection_width bounds ----------------------------- +expect_error( + cluster_correction_l(avg_min_distance = 1000, eff_detection_width = 0), + 'eff_detection_width must be greater than 0' +) +expect_error( + cluster_correction_l(avg_min_distance = 1000, eff_detection_width = -500), + 'eff_detection_width must be greater than 0' +) + +## Input validation: NA eff_detection_width --------------------------------- +expect_error( + cluster_correction_l(avg_min_distance = 1000, eff_detection_width = NA), + 'eff_detection_width must be a single numeric value' +) \ No newline at end of file From 31f4a7c172695e1f9d9ff1ad382b12a04bb005c5 Mon Sep 17 00:00:00 2001 From: sahinya17 Date: Tue, 14 Apr 2026 10:38:44 +1000 Subject: [PATCH 04/10] Added tests for cluster correction functions --- inst/tinytest/test_cluster_correction.R | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/inst/tinytest/test_cluster_correction.R b/inst/tinytest/test_cluster_correction.R index 6f95965..04d0583 100644 --- a/inst/tinytest/test_cluster_correction.R +++ b/inst/tinytest/test_cluster_correction.R @@ -1,6 +1,6 @@ # Tests for cluster_correction_a and cluster_correction_l -# Helper: simple 2-turbine data.frame +# Helper df two_turbines <- data.frame(lon = c(0, 0.01), lat = c(51, 51)) # cluster_correction_a---------------------- @@ -12,8 +12,6 @@ expect_equal(length(result_a), 1L) expect_true(is.finite(result_a)) ## Single turbine: correction factor should be ~1 --------------------------- -# sf computes spherical area, which differs slightly from planar pi*r^2, -# so we allow a small tolerance rather than expecting exactly 1. one_turbine <- data.frame(lon = 0, lat = 51) expect_equal( cluster_correction_a(one_turbine, eff_detection_width = 500), @@ -22,15 +20,11 @@ expect_equal( ) ## Very widely spaced turbines: factor close to 1 --------------------------- -# Turbines far apart so EDAs don't overlap; each turbine contributes one full EDA. far_turbines <- data.frame(lon = c(0, 10), lat = c(51, 51)) result_far <- cluster_correction_a(far_turbines, eff_detection_width = 500) expect_equal(result_far, 1, tolerance = 0.05) ## Overlapping EDAs: factor lower than non-overlapping case ----------------- -# Turbines very close together relative to detection width; their EDAs nearly -# fully overlap, so the correction factor should be much less than for -# non-overlapping turbines. close_turbines <- data.frame(lon = c(0, 0.0001), lat = c(51, 51)) result_close <- cluster_correction_a(close_turbines, eff_detection_width = 5000) expect_true(result_close < result_far) @@ -116,8 +110,6 @@ expect_error( ) ## Input validation: eff_detection_width type/length ----------------------- -# Note: a character value causes an arithmetic error before the stopifnot fires, -# because EDA = pi*(eff_detection_width/2)^2 is evaluated first. expect_error(cluster_correction_a(two_turbines, eff_detection_width = '500')) expect_error( cluster_correction_a(two_turbines, eff_detection_width = c(500, 600)), @@ -135,7 +127,6 @@ expect_error( ) ## Input validation: NA eff_detection_width -------------------------------- -# is.numeric(NA) is FALSE, so the "single numeric value" check fires first. expect_error( cluster_correction_a(two_turbines, eff_detection_width = NA), 'eff_detection_width must be a single numeric value' From 77ea8731840754394140cb66f64f28fea294a7fc Mon Sep 17 00:00:00 2001 From: techisdead Date: Wed, 15 Apr 2026 15:51:52 +1000 Subject: [PATCH 05/10] test failing test --- inst/tinytest/test_flux_scaling.R | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/inst/tinytest/test_flux_scaling.R b/inst/tinytest/test_flux_scaling.R index 4b37705..b1e158c 100644 --- a/inst/tinytest/test_flux_scaling.R +++ b/inst/tinytest/test_flux_scaling.R @@ -5,6 +5,7 @@ # obs_flux obs_flux <- 5 +spatial_correction <- 0.8 rotor_diameter <- 300 hub_height <- 1 ### Invalid Inputs----------------- @@ -52,13 +53,14 @@ expect_error(turbine_flights(obs_flux, rotor_diameter, -100) , "variable out of bounds") ## Zero obs_flux or diameter gives 0 turbine_flights---------- -expect_equal(turbine_flights(0, rotor_diameter, hub_height), 0) -expect_equal(turbine_flights(obs_flux, 0, hub_height), 0) -expect_equal(turbine_flights(0, 0, hub_height), 0) +expect_equal(turbine_flights(0, rotor_diameter, spatial_correction, hub_height), 0) +expect_equal(turbine_flights(obs_flux, 0, spatial_correction, hub_height), 0) +expect_equal(turbine_flights(0, 0, 0, hub_height), 0) ## Valid inputs give the correct result--------- -expect_equal(turbine_flights(obs_flux, rotor_diameter, hub_height) - , obs_flux * rotor_diameter * (hub_height + 0.5*rotor_diameter)) +expect_equal(turbine_flights(obs_flux, spatial_correction,rotor_diameter, hub_height) + , obs_flux * spatial_correction * + rotor_diameter * (hub_height + 0.5*rotor_diameter)) # flights per year-------------------------- ## Test if data inputs are checked appropriately---------------- From ce4911a67fcf87c4028192989ba6d7f216f15fcf Mon Sep 17 00:00:00 2001 From: techisdead Date: Fri, 17 Apr 2026 18:33:20 +1000 Subject: [PATCH 06/10] smoother cluster code in scratch - to be done in package --- R/cluster_correction.R | 32 +++++++++++-------- R/flux_scaling.R | 5 ++- R/scratch.R | 72 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 95 insertions(+), 14 deletions(-) diff --git a/R/cluster_correction.R b/R/cluster_correction.R index 9f78952..77fbb01 100644 --- a/R/cluster_correction.R +++ b/R/cluster_correction.R @@ -1,13 +1,17 @@ #' Correct for turbine clusters #' -#' TODO description +#' 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 EDR is large, then the -#' effective detection area contains more than one turbine. +#' 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. #' The observed flux needs to be weighted to account for the effective turbines -#' in the observed area. This accounts for the fact that an observed flight can -#' only interact with one turbine at any given time. TODO +#' in the observed area. This accounts for the fact that the observed flux is an instantaneous measure and cannot interact with multiple turbinesat once. +#' +#' This correction is derived from different assumptions but fulfils a similar model correction as the \code{sqrt(n_turbines)} in #' #' @param df_turbines data.frame; a data.frame with one row per turbine, #' containing the coordinates of each turbine in WGS 84 @@ -25,12 +29,10 @@ cluster_correction_a <- function(df_turbines, eff_detection_width){ # lat column can be "lat", "latitude" or "y" (any case) # lon column can be "lon", "long", "longitude" or "x" (any case) # currently have it requiring WGS 84 but might be better to add a CRS input? - - # TODO (low priority) I think this needs to be optimised better - it's kinda - # slow when doing it in a loop in the stochastic example - # to do it mathematically this is a good reference - # https://stackoverflow.com/questions/1667310/combined-area-of-overlapping-circles - # but I don't want to spend my day on that rn + # + # + # + # stopifnot( "df_turbines must be a data.frame" = @@ -47,11 +49,15 @@ cluster_correction_a <- function(df_turbines, eff_detection_width){ lon_col <- names(df_turbines)[tolower(names(df_turbines)) %in% c("lon", "longitude", "long", - "x")] + "x", + "easting", + "northing")] lat_col <- names(df_turbines)[tolower(names(df_turbines)) %in% c("lat", "latitude", - "y")] + "y", + "easting", + "northing")] stopifnot( "df_turbines must contain a longitude column named lon, long, longitude, or x" = diff --git a/R/flux_scaling.R b/R/flux_scaling.R index 3f6497c..42cb0ef 100644 --- a/R/flux_scaling.R +++ b/R/flux_scaling.R @@ -32,7 +32,10 @@ #' @example examples/flux_example.R #' #' @export -turbine_flights <- function(obs_flux, spatial_correction = 1, 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) diff --git a/R/scratch.R b/R/scratch.R index e69de29..d5d923b 100644 --- a/R/scratch.R +++ b/R/scratch.R @@ -0,0 +1,72 @@ +library(sf) + +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) +) + +crs = 4326 +eff_detection_width <- 500 +n_turbines <- nrow(df_turbines) + +area_separate <- nrow(df_turbines)* pi*(eff_detection_width/2)^2 + + +sf_turbines <- sf::st_as_sf(df_turbines, coords = c("lon", "lat"), crs = crs) + +sf::st_is_longlat(sf_turbines) +projcrs <- lonlat_to_utm(sf::st_coordinates(sf_turbines)) + +sf_turbines <- st_transform(sf_turbines, crs = projcrs) + +sf_turb_buffer <- sf::st_make_valid( + sf::st_buffer(x = sf_turbines, + dist = eff_detection_width/2, + nQuadSegs = 360)) + +plot(sf_turb_buffer) + +sf_turb_union <- sf::st_make_valid(sf::st_union(sf_turb_buffer)) +plot(sf_turb_union) + +turb_area <- units::drop_units(sf::st_area(sf_turb_union)) + +turb_area + +# == version 2 +sf_turb_intersect <- sf::st_make_valid(sf::st_intersection(sf_turb_buffer)) +plot(sf_turb_intersect) + +intersect_area <- units::drop_units(sf::st_area(sf_turb_intersect[sf_turb_intersect$n.overlaps > 1,])) + +turb_area2 <- area_separate - intersect_area +turb_area2 + + + +sf::st_distance(sf_turbines) + +str(sf_turbines) + + + + + +#' lon lat to UTM +#' +#' Equation from +#' +#' +lonlat_to_utm = function(lonlat) { + + utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1 + if (lonlat[2] > 0) { + utm + 32600 + } else{ + utm + 32700 + } + +} + + From 9832ca821c0c6aa6bcf7f3c27e539b0c176e91c4 Mon Sep 17 00:00:00 2001 From: techisdead Date: Mon, 20 Apr 2026 20:36:34 +1000 Subject: [PATCH 07/10] refactored the cluster correction. note the order of inputs have changed and I changed eff_width to be the half width to line up with its use in other functions. These tests pass but others fail - to be checked --- R/cluster_correction.R | 352 ++++++++++++++---------- R/scratch.R | 154 ++++++----- R/utils_spatial.R | 15 + inst/REFERENCES.bib | 9 + inst/tinytest/test_cluster_correction.R | 92 +++---- inst/tinytest/test_flux_wrapper.R | 2 +- 6 files changed, 363 insertions(+), 261 deletions(-) create mode 100644 R/utils_spatial.R diff --git a/R/cluster_correction.R b/R/cluster_correction.R index 77fbb01..d7e2c27 100644 --- a/R/cluster_correction.R +++ b/R/cluster_correction.R @@ -1,4 +1,4 @@ -#' Correct for turbine clusters +#' 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 @@ -7,182 +7,246 @@ #' 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. -#' 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 turbinesat once. +#' 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 correction is derived from different assumptions but fulfils a similar model correction as the \code{sqrt(n_turbines)} in +#' This package provides a 1D line correction [cluster_correction_l()] and 2D aerial [cluster_correction_a()]. #' -#' @param df_turbines data.frame; a data.frame with one row per turbine, -#' containing the coordinates of each turbine in WGS 84 -#' @param eff_detection_width numeric; the effective detection width, -#' which is usually 2 x effective detection radius. +#' ## [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. +#' +#' +#' +#' #' -#' @return numeric; the cluster correction factor - the average number -#' of turbines per EDA (effective detection area). +#' @name cluster_correction + + + +#' @rdname cluster_correction #' -#' @importFrom sf st_as_sf st_buffer st_make_valid st_union st_area -#' @importFrom units drop_units +#' @param eff_detection_width numeric; the effective detection (half) width or effective detection radius +#' @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(df_turbines, eff_detection_width){ - # df_turbines needs turbine ids and lat/lon - # lat column can be "lat", "latitude" or "y" (any case) - # lon column can be "lon", "long", "longitude" or "x" (any case) - # currently have it requiring WGS 84 but might be better to add a CRS input? - # - # - # - # - - 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 - ) +cluster_correction_a <- function(eff_detection_width, df_turbines, crs = 4326){ - EDA <- pi*(eff_detection_width/2)^2 - n_turbines <- nrow(df_turbines) + 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) - lon_col <- names(df_turbines)[tolower(names(df_turbines)) %in% c("lon", - "longitude", - "long", - "x", - "easting", - "northing")] - lat_col <- names(df_turbines)[tolower(names(df_turbines)) %in% c("lat", - "latitude", - "y", - "easting", - "northing")] + sf_turbines <- df_to_projsf(df_turbines, crs) # internal helper function - stopifnot( - "df_turbines must contain a longitude column named lon, long, longitude, or x" = - length(lon_col) >= 1 - ) - stopifnot( - "df_turbines must contain a latitude column named lat, latitude, or y" = - length(lat_col) >= 1 - ) - - if (length(lon_col) > 1) { - warning( - "multiple possible longitude columns detected: ", - paste(lon_col, collapse = ", "), - ". Using the first match: ", lon_col[1] - ) - lon_col <- lon_col[1] - } - if (length(lat_col) > 1) { - warning( - "multiple possible latitude columns detected: ", - paste(lat_col, collapse = ", "), - ". Using the first match: ", lat_col[1] - ) - lat_col <- lat_col[1] - } - - stopifnot( - "NA values detected in longitude column" = - !any(is.na(df_turbines[[lon_col]])) - ) - stopifnot( - "NA values detected in latitude column" = - !any(is.na(df_turbines[[lat_col]])) - ) - - stopifnot( - "longitude values appear out of WGS 84 range (-180 to 180)" = - all(df_turbines[[lon_col]] >= -180 & df_turbines[[lon_col]] <= 180) - ) - stopifnot( - "latitude values appear out of WGS 84 range (-90 to 90)" = - all(df_turbines[[lat_col]] >= -90 & df_turbines[[lat_col]] <= 90) - ) - - stopifnot( - "eff_detection_width must be a single numeric value" = - is.numeric(eff_detection_width) && length(eff_detection_width) == 1 - ) - stopifnot( - "eff_detection_width must be greater than 0" = - eff_detection_width > 0 - ) - stopifnot( - "NA detected in eff_detection_width" = - !is.na(eff_detection_width) - ) + ## total area with no overlap + area_separate <- nrow(df_turbines)* pi*(eff_detection_width)^2 #half width==radius - sf_turb <- sf::st_as_sf(df_turbines, coords = c(lon_col, lat_col), crs = 4326) - + ## area buffered by EDA sf_turb_buffer <- sf::st_make_valid( - sf::st_buffer(x = sf_turb, - dist = eff_detection_width/2)) + sf::st_buffer(x = sf_turbines, + dist = eff_detection_width, + nQuadSegs = 360)) + ## union sf_turb_union <- sf::st_make_valid(sf::st_union(sf_turb_buffer)) - turb_area <- units::drop_units(sf::st_area(sf_turb_union)) - - turbine_prob <- round(turb_area/(n_turbines*EDA), 16) # number of turbines that flux is effectively spread across - # rounding because sf sometimes has float issues + ## Area + turb_area_buf <- sf::st_area(sf_turb_union) |> as.numeric() - # may need to add a check to ensure it's always <= 1? I kind of can't think of - # a way that could happen but I think this function needs some testing to see - # if you can break it in that way + ## return + turbine_prob <- min( turb_area_buf/area_separate, 1.0 ) return(turbine_prob) } -#' Correct for turbine distance (offshore) +#' @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 #' -#' TODO description +#' @examples #' -#' @details -#' TODO +#' # generate input data #' -#' @param avg_min_distance numeric; the average distance between each turbine and its nearest neighbour -#' @param eff_detection_width numeric; the effective detection width, -#' which is usually 2 x effective detection radius. -#' -#' @return numeric; the (linear) cluster correction factor - the average number -#' of turbines per ESW (effective strip width). +#' +#' 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_col, lat_col), +#' 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(avg_min_distance, eff_detection_width){ - # I think this is literally it? - # TODO: I would like to have the option to input df_turbines here and have it calculate the avg min distance - # but I think that requires better knowledge of R to do in a relatively speedy way - # ES suggested spatialdatatable::dtHaversine? +cluster_correction_l <- function(eff_detection_width = NULL, + 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 eff_detection_width or avg_min_distance") + } + + if( !is.null(avg_min_distance) & !is.null(df_turbines)){ + stop("Please enter only one of eff_detection_width and avg_min_distance") + } + + + 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( - "avg_min_distance must be a single numeric value" = - is.numeric(avg_min_distance) && length(avg_min_distance) == 1 - ) - stopifnot( - "NA detected in avg_min_distance" = - !is.na(avg_min_distance) + "df_turbines must be a data.frame" = + is.data.frame(df_turbines) ) stopifnot( - "avg_min_distance must be greater than 0" = - avg_min_distance > 0 + "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( - "eff_detection_width must be a single numeric value" = - is.numeric(eff_detection_width) && length(eff_detection_width) == 1 - ) - stopifnot( - "NA detected in eff_detection_width" = - !is.na(eff_detection_width) + "df_turbines must contain one longitude column named lon, long, longitude, x, or easting" = + length(lon_col) == 1 ) stopifnot( - "eff_detection_width must be greater than 0" = - eff_detection_width > 0 + "df_turbines must contain one latitude column named lat, latitude, y, or northing" = + length(lat_col) == 1 ) - return(eff_detection_width/avg_min_distance) -} + 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) + +} \ No newline at end of file diff --git a/R/scratch.R b/R/scratch.R index d5d923b..5b37de1 100644 --- a/R/scratch.R +++ b/R/scratch.R @@ -1,72 +1,92 @@ -library(sf) - -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) -) - -crs = 4326 -eff_detection_width <- 500 -n_turbines <- nrow(df_turbines) - -area_separate <- nrow(df_turbines)* pi*(eff_detection_width/2)^2 - - -sf_turbines <- sf::st_as_sf(df_turbines, coords = c("lon", "lat"), crs = crs) - -sf::st_is_longlat(sf_turbines) -projcrs <- lonlat_to_utm(sf::st_coordinates(sf_turbines)) - -sf_turbines <- st_transform(sf_turbines, crs = projcrs) - -sf_turb_buffer <- sf::st_make_valid( - sf::st_buffer(x = sf_turbines, - dist = eff_detection_width/2, - nQuadSegs = 360)) - -plot(sf_turb_buffer) - -sf_turb_union <- sf::st_make_valid(sf::st_union(sf_turb_buffer)) -plot(sf_turb_union) - -turb_area <- units::drop_units(sf::st_area(sf_turb_union)) - -turb_area - -# == version 2 -sf_turb_intersect <- sf::st_make_valid(sf::st_intersection(sf_turb_buffer)) -plot(sf_turb_intersect) - -intersect_area <- units::drop_units(sf::st_area(sf_turb_intersect[sf_turb_intersect$n.overlaps > 1,])) - -turb_area2 <- area_separate - intersect_area -turb_area2 - - - -sf::st_distance(sf_turbines) - -str(sf_turbines) - - - - - -#' lon lat to UTM +#' library(sf) +#' +#' 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) +#' ) +#' +#' crs = 4326 +#' eff_detection_width <- 500 +#' n_turbines <- nrow(df_turbines) +#' +#' area_separate <- nrow(df_turbines)* pi*(eff_detection_width)^2 +#' #' -#' Equation from +#' sf_turbines <- sf::st_as_sf(df_turbines, coords = c("lon", "lat"), crs = crs) #' +#' sf::st_is_longlat(sf_turbines) +#' projcrs <- lonlat_to_utm(sf::st_coordinates(sf_turbines)) #' -lonlat_to_utm = function(lonlat) { - - utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1 - if (lonlat[2] > 0) { - utm + 32600 - } else{ - utm + 32700 - } - -} +#' sf_turbines <- st_transform(sf_turbines, crs = projcrs) +#' +#' sf_turb_buffer <- sf::st_make_valid( +#' sf::st_buffer(x = sf_turbines, +#' dist = eff_detection_width/2, +#' nQuadSegs = 360)) +#' +#' plot(sf_turb_buffer) +#' +#' sf_turb_union <- sf::st_make_valid(sf::st_union(sf_turb_buffer)) +#' plot(sf_turb_union) +#' +#' turb_area <- units::drop_units(sf::st_area(sf_turb_union)) +#' +#' turb_area +#' +#' # == version 2 +#' sf_turb_intersect <- sf::st_make_valid(sf::st_intersection(sf_turb_buffer)) +#' plot(sf_turb_intersect) +#' +#' intersect_area <- units::drop_units(sf::st_area(sf_turb_intersect[sf_turb_intersect$n.overlaps > 1,])) +#' +#' turb_area2 <- area_separate - intersect_area +#' turb_area2 +#' +#' +#' +#' sf::st_distance(sf_turbines) +#' +#' str(sf_turbines) +#' +#' +#' +#' +#' +#' #' lon lat to UTM +#' #' +#' #' Equation from +#' #' +#' #' +#' lonlat_to_utm = function(lonlat) { +#' +#' utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1 +#' if (lonlat[2] > 0) { +#' utm + 32600 +#' } else{ +#' utm + 32700 +#' } +#' +#' } +#' +#' +#' +#' 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) +#' ) +#' +#' sf_turbines <- sf::st_as_sf(df_turbines, coords = c("lon", "lat"), crs = 4326) +#' +#' all(sf::st_is_valid(sf_turbines)) +#' +#' stopifnot("dataset contains invalid coordinates. Check values and crs", +#' all(sf::st_is_valid(sf_turbines)) +#' ) +#' +#' +#' + diff --git a/R/utils_spatial.R b/R/utils_spatial.R new file mode 100644 index 0000000..9d03150 --- /dev/null +++ b/R/utils_spatial.R @@ -0,0 +1,15 @@ +#' lon lat to UTM +#' +#' Equation from https://r.geocompx.org/reproj-geo-data#which-crs +#' +#' +lonlat_to_utm = function(lonlat) { + + utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1 + if (lonlat[2] > 0) { + utm + 32600 + } else{ + utm + 32700 + } + +} \ No newline at end of file diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 0ef07f7..25b8d38 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -49,3 +49,12 @@ @article{Wilson1927 doi = {10.1080/01621459.1927.10502953}, URL = {https://www.jstor.org/stable/2276774} } + +@techreport{Band2012, + title = {Using a Collision Risk Model to Assess Bird Collision Risks for Offshore Windfarms}, + author = {Band, Bill}, + year = 2012, + institution = {prepared for The Crown Estate as part of the Strategic Ornithological Support Services programme, project SOSS02.}, + url = {https://www.bto.org/sites/default/files/u28/downloads/Projects/Final_Report_SOSS02_Band1ModelGuidance.pdf}, + langid = {english} +} diff --git a/inst/tinytest/test_cluster_correction.R b/inst/tinytest/test_cluster_correction.R index 04d0583..5cef3fc 100644 --- a/inst/tinytest/test_cluster_correction.R +++ b/inst/tinytest/test_cluster_correction.R @@ -6,27 +6,28 @@ two_turbines <- data.frame(lon = c(0, 0.01), lat = c(51, 51)) # cluster_correction_a---------------------- ## Returns a single numeric value ------------------------------------------- -result_a <- cluster_correction_a(two_turbines, eff_detection_width = 500) +result_a <- cluster_correction_a(df_turbines = two_turbines, eff_detection_width = 500) expect_true(is.numeric(result_a)) expect_equal(length(result_a), 1L) expect_true(is.finite(result_a)) ## Single turbine: correction factor should be ~1 --------------------------- one_turbine <- data.frame(lon = 0, lat = 51) + expect_equal( - cluster_correction_a(one_turbine, eff_detection_width = 500), + cluster_correction_a(eff_detection_width = 500, one_turbine), 1, tolerance = 0.05 ) ## Very widely spaced turbines: factor close to 1 --------------------------- -far_turbines <- data.frame(lon = c(0, 10), lat = c(51, 51)) -result_far <- cluster_correction_a(far_turbines, eff_detection_width = 500) +far_turbines <- data.frame(lon = c(1, 10), lat = c(51, 51)) +result_far <- cluster_correction_a(df_turbines = far_turbines, eff_detection_width = 500) expect_equal(result_far, 1, tolerance = 0.05) ## Overlapping EDAs: factor lower than non-overlapping case ----------------- close_turbines <- data.frame(lon = c(0, 0.0001), lat = c(51, 51)) -result_close <- cluster_correction_a(close_turbines, eff_detection_width = 5000) +result_close <- cluster_correction_a(df_turbines = close_turbines, eff_detection_width = 5000) expect_true(result_close < result_far) ## Result is positive and finite -------------------------------------------- @@ -36,23 +37,23 @@ expect_true(result_far > 0) ## Accepts alternative column names ----------------------------------------- df_alt_names <- data.frame(longitude = c(0, 0.01), latitude = c(51, 51)) expect_equal( - cluster_correction_a(df_alt_names, 500), + cluster_correction_a(500, df_turbines = df_alt_names), result_a, tolerance = 1e-6 ) df_xy <- data.frame(x = c(0, 0.01), y = c(51, 51)) -expect_equal(cluster_correction_a(df_xy, 500), result_a, tolerance = 1e-6) +expect_equal(cluster_correction_a(500, df_xy), result_a, tolerance = 1e-6) ## Case-insensitive column matching ----------------------------------------- df_upper <- data.frame(LON = c(0, 0.01), LAT = c(51, 51)) -expect_equal(cluster_correction_a(df_upper, 500), result_a, tolerance = 1e-6) +expect_equal(cluster_correction_a(500, df_upper), result_a, tolerance = 1e-6) ## Warns when multiple possible lon/lat columns detected -------------------- df_multi_lon <- data.frame(lon = c(0, 0.01), long = c(0, 0.01), lat = c(51, 51)) -expect_warning( - cluster_correction_a(df_multi_lon, 500), - 'multiple possible longitude columns detected' +expect_error( + cluster_correction_a(500, df_multi_lon), + 'df_turbines must contain one longitude column named lon, long, longitude, x, or easting' ) df_multi_lat <- data.frame( @@ -60,75 +61,68 @@ df_multi_lat <- data.frame( lat = c(51, 51), latitude = c(51, 51) ) -expect_warning( - cluster_correction_a(df_multi_lat, 500), - 'multiple possible latitude columns detected' +expect_error( + cluster_correction_a(500, df_multi_lat) ) ## Input validation: non-data.frame ----------------------------------------- expect_error( - cluster_correction_a(list(lon = 0, lat = 51), 500), + cluster_correction_a(500, list(lon = 0, lat = 51)), 'df_turbines must be a data.frame' ) ## Input validation: zero rows --------------------------------------------- expect_error( - cluster_correction_a(data.frame(lon = numeric(0), lat = numeric(0)), 500), + cluster_correction_a(500, data.frame(lon = numeric(0), lat = numeric(0))), 'df_turbines must have at least one row' ) ## Input validation: missing lon column ------------------------------------ expect_error( - cluster_correction_a(data.frame(easting = c(0, 1), lat = c(51, 51)), 500), - 'df_turbines must contain a longitude column' + cluster_correction_a(500, data.frame(longit = c(0, 1), lat = c(51, 51))), + 'df_turbines must contain one longitude column named lon, long, longitude, x, or easting' ) ## Input validation: missing lat column ------------------------------------ expect_error( - cluster_correction_a(data.frame(lon = c(0, 1), northing = c(51, 51)), 500), - 'df_turbines must contain a latitude column' + cluster_correction_a(500, data.frame(lon = c(0, 1), lati = c(51, 51))) ) ## Input validation: NA coordinates ---------------------------------------- expect_error( - cluster_correction_a(data.frame(lon = c(0, NA), lat = c(51, 51)), 500), - 'NA values detected in longitude column' + cluster_correction_a(500, data.frame(lon = c(0, NA), lat = c(51, 51))), + 'missing values in coordinates not allowed' ) expect_error( - cluster_correction_a(data.frame(lon = c(0, 0.01), lat = c(51, NA)), 500), - 'NA values detected in latitude column' + cluster_correction_a(500, data.frame(lon = c(0, 0.01), lat = c(51, NA))), + 'missing values in coordinates not allowed' ) ## Input validation: out-of-range coordinates ------------------------------ expect_error( - cluster_correction_a(data.frame(lon = c(0, 200), lat = c(51, 51)), 500), - 'longitude values appear out of WGS 84 range' -) -expect_error( - cluster_correction_a(data.frame(lon = c(0, 0.01), lat = c(51, 100)), 500), - 'latitude values appear out of WGS 84 range' + cluster_correction_a(500, data.frame(lon = c(0, -10), lat = c(51, 51))), ) ## Input validation: eff_detection_width type/length ----------------------- -expect_error(cluster_correction_a(two_turbines, eff_detection_width = '500')) +expect_error(cluster_correction_a(df_turbines = two_turbines, eff_detection_width = '500')) expect_error( - cluster_correction_a(two_turbines, eff_detection_width = c(500, 600)), + cluster_correction_a(df_turbines = two_turbines, eff_detection_width = c(500, 600)), 'eff_detection_width must be a single numeric value' ) ## Input validation: eff_detection_width bounds ---------------------------- expect_error( - cluster_correction_a(two_turbines, eff_detection_width = 0), - 'eff_detection_width must be greater than 0' + cluster_correction_a(df_turbines = two_turbines, eff_detection_width = 0), + 'variable out of bounds' ) expect_error( - cluster_correction_a(two_turbines, eff_detection_width = -100), - 'eff_detection_width must be greater than 0' + cluster_correction_a(df_turbines = two_turbines, eff_detection_width = -100), + 'variable out of bounds' ) ## Input validation: NA eff_detection_width -------------------------------- expect_error( - cluster_correction_a(two_turbines, eff_detection_width = NA), + cluster_correction_a(df_turbines = two_turbines, eff_detection_width = NA), 'eff_detection_width must be a single numeric value' ) @@ -137,7 +131,7 @@ expect_error( ## Basic correctness: result equals eff_detection_width / avg_min_distance---- expect_equal( cluster_correction_l(avg_min_distance = 1000, eff_detection_width = 500), - 0.5 + 1 ) expect_equal( cluster_correction_l(avg_min_distance = 500, eff_detection_width = 500), @@ -145,18 +139,18 @@ expect_equal( ) expect_equal( cluster_correction_l(avg_min_distance = 200, eff_detection_width = 500), - 2.5 + 0.4 ) ## Returns single numeric --------------------------------------------------- -res_l <- cluster_correction_l(1000, 500) +res_l <- cluster_correction_l(500,1000) expect_true(is.numeric(res_l)) expect_equal(length(res_l), 1L) ## Input validation: avg_min_distance type/length --------------------------- expect_error( cluster_correction_l(avg_min_distance = '1000', eff_detection_width = 500), - 'avg_min_distance must be a single numeric value' + 'Numeric input expected' ) expect_error( cluster_correction_l( @@ -169,23 +163,23 @@ expect_error( ## Input validation: avg_min_distance bounds -------------------------------- expect_error( cluster_correction_l(avg_min_distance = 0, eff_detection_width = 500), - 'avg_min_distance must be greater than 0' + 'variable out of bounds' ) expect_error( cluster_correction_l(avg_min_distance = -100, eff_detection_width = 500), - 'avg_min_distance must be greater than 0' + 'variable out of bounds' ) ## Input validation: NA avg_min_distance ------------------------------------ expect_error( cluster_correction_l(avg_min_distance = NA, eff_detection_width = 500), - 'avg_min_distance must be a single numeric value' + 'Numeric input expected' ) ## Input validation: eff_detection_width type/length ------------------------ expect_error( cluster_correction_l(avg_min_distance = 1000, eff_detection_width = '500'), - 'eff_detection_width must be a single numeric value' + 'Numeric input expected' ) expect_error( cluster_correction_l( @@ -198,15 +192,15 @@ expect_error( ## Input validation: eff_detection_width bounds ----------------------------- expect_error( cluster_correction_l(avg_min_distance = 1000, eff_detection_width = 0), - 'eff_detection_width must be greater than 0' + 'variable out of bounds' ) expect_error( cluster_correction_l(avg_min_distance = 1000, eff_detection_width = -500), - 'eff_detection_width must be greater than 0' + 'variable out of bounds' ) ## Input validation: NA eff_detection_width --------------------------------- expect_error( cluster_correction_l(avg_min_distance = 1000, eff_detection_width = NA), - 'eff_detection_width must be a single numeric value' -) \ No newline at end of file + 'Numeric input expected' +) diff --git a/inst/tinytest/test_flux_wrapper.R b/inst/tinytest/test_flux_wrapper.R index bf7ec6c..ce4c2cc 100644 --- a/inst/tinytest/test_flux_wrapper.R +++ b/inst/tinytest/test_flux_wrapper.R @@ -25,7 +25,7 @@ expect_error( rotor_diameter, hub_height, prop_day, - prop_year + prop_year, ), "Only point transects are currently implemented. Line transects and digital areal surveys are in development." From d466e0b34c328d95ccc9ed78be12336d16392eb8 Mon Sep 17 00:00:00 2001 From: techisdead Date: Tue, 21 Apr 2026 15:33:53 +1000 Subject: [PATCH 08/10] reviewed by PY --- NAMESPACE | 6 -- R/cluster_correction.R | 11 ++-- R/observed_flux.R | 3 +- man/cluster_correction.Rd | 110 ++++++++++++++++++++++++++++++++++++ man/cluster_correction_a.Rd | 29 ---------- man/cluster_correction_l.Rd | 24 -------- man/lonlat_to_utm.Rd | 11 ++++ 7 files changed, 129 insertions(+), 65 deletions(-) create mode 100644 man/cluster_correction.Rd delete mode 100644 man/cluster_correction_a.Rd delete mode 100644 man/cluster_correction_l.Rd create mode 100644 man/lonlat_to_utm.Rd diff --git a/NAMESPACE b/NAMESPACE index 7e5b3e1..27bd949 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,9 +26,3 @@ export(vertical_flux_plane_area) export(w_from_distmodel) import(Distance) importFrom(Rdpack,reprompt) -importFrom(sf,st_area) -importFrom(sf,st_as_sf) -importFrom(sf,st_buffer) -importFrom(sf,st_make_valid) -importFrom(sf,st_union) -importFrom(units,drop_units) diff --git a/R/cluster_correction.R b/R/cluster_correction.R index d7e2c27..c5d3432 100644 --- a/R/cluster_correction.R +++ b/R/cluster_correction.R @@ -35,7 +35,8 @@ #' @rdname cluster_correction #' -#' @param eff_detection_width numeric; the effective detection (half) width or effective detection radius +#' @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) @@ -69,12 +70,12 @@ cluster_correction_a <- function(eff_detection_width, df_turbines, crs = 4326){ 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)^2 #half width==radius + 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, + dist = eff_detection_width*0.5, nQuadSegs = 360)) ## union @@ -138,7 +139,7 @@ cluster_correction_a <- function(eff_detection_width, df_turbines, crs = 4326){ #' corr1; corr2 #' #' @export -cluster_correction_l <- function(eff_detection_width = NULL, +cluster_correction_l <- function(eff_detection_width, avg_min_distance = NULL, df_turbines = NULL, crs = 4326){ @@ -244,7 +245,7 @@ df_to_projsf <- function(df_turbines, crs){ } 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 ) + check_num_bounds(abs(sf::st_coordinates(sf_turbines)[, "Y"]), 1.e6, 10.e6 ) return(sf_turbines) diff --git a/R/observed_flux.R b/R/observed_flux.R index 7e36160..cfced0d 100644 --- a/R/observed_flux.R +++ b/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 diff --git a/man/cluster_correction.Rd b/man/cluster_correction.Rd new file mode 100644 index 0000000..38f66e2 --- /dev/null +++ b/man/cluster_correction.Rd @@ -0,0 +1,110 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cluster_correction.R +\name{cluster_correction} +\alias{cluster_correction} +\alias{cluster_correction_a} +\alias{cluster_correction_l} +\title{Correct for turbine clustering} +\usage{ +cluster_correction_a(eff_detection_width, df_turbines, crs = 4326) + +cluster_correction_l( + eff_detection_width = NULL, + avg_min_distance = NULL, + df_turbines = NULL, + crs = 4326 +) +} +\arguments{ +\item{eff_detection_width}{numeric; the effective detection (half) width or effective detection radius} + +\item{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)} + +\item{crs}{optional coordinate system epsg code for the coordinates. Only used if df_turbines provided. Defaults to crs = 4326 (WGS 84 - lat/lon).} + +\item{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} +} +\value{ +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. + +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 +} +\description{ +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 \code{\link[=cluster_correction_l]{cluster_correction_l()}} and 2D aerial \code{\link[=cluster_correction_a]{cluster_correction_a()}}. +\subsection{\code{\link[=cluster_correction_l]{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. +} + +\subsection{\code{\link[=cluster_correction_a]{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. +} + +\subsection{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. +} +} +\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) + + +# 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_col, lat_col), + 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 + +} diff --git a/man/cluster_correction_a.Rd b/man/cluster_correction_a.Rd deleted file mode 100644 index 3aa84b8..0000000 --- a/man/cluster_correction_a.Rd +++ /dev/null @@ -1,29 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/cluster_correction.R -\name{cluster_correction_a} -\alias{cluster_correction_a} -\title{Correct for turbine clusters} -\usage{ -cluster_correction_a(df_turbines, eff_detection_width) -} -\arguments{ -\item{df_turbines}{data.frame; a data.frame with one row per turbine, -containing the coordinates of each turbine in WGS 84} - -\item{eff_detection_width}{numeric; the effective detection width, -which is usually 2 x effective detection radius.} -} -\value{ -numeric; the cluster correction factor - the average number -of turbines per EDA (effective detection area). -} -\description{ -TODO description -} -\details{ -If the turbine layout is more compact and/or your EDR is large, then the -effective detection area contains more than one turbine. -The observed flux needs to be weighted to account for the effective turbines -in the observed area. This accounts for the fact that an observed flight can -only interact with one turbine at any given time. TODO -} diff --git a/man/cluster_correction_l.Rd b/man/cluster_correction_l.Rd deleted file mode 100644 index a9c105c..0000000 --- a/man/cluster_correction_l.Rd +++ /dev/null @@ -1,24 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/cluster_correction.R -\name{cluster_correction_l} -\alias{cluster_correction_l} -\title{Correct for turbine distance (offshore)} -\usage{ -cluster_correction_l(avg_min_distance, eff_detection_width) -} -\arguments{ -\item{avg_min_distance}{numeric; the average distance between each turbine and its nearest neighbour} - -\item{eff_detection_width}{numeric; the effective detection width, -which is usually 2 x effective detection radius.} -} -\value{ -numeric; the (linear) cluster correction factor - the average number -of turbines per ESW (effective strip width). -} -\description{ -TODO description -} -\details{ -TODO -} diff --git a/man/lonlat_to_utm.Rd b/man/lonlat_to_utm.Rd new file mode 100644 index 0000000..02f78fc --- /dev/null +++ b/man/lonlat_to_utm.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils_spatial.R +\name{lonlat_to_utm} +\alias{lonlat_to_utm} +\title{lon lat to UTM} +\usage{ +lonlat_to_utm(lonlat) +} +\description{ +Equation from https://r.geocompx.org/reproj-geo-data#which-crs +} From ce83780f0905dc3fdcb21f5c4234ef9178aedc99 Mon Sep 17 00:00:00 2001 From: techisdead Date: Wed, 22 Apr 2026 15:27:30 +1000 Subject: [PATCH 09/10] bumped version, checks and tests pass. Redocumented. Ready to merge --- DESCRIPTION | 2 +- NEWS.md | 3 ++- R/cluster_correction.R | 2 +- R/flux_scaling.R | 1 - R/flux_wrapper.R | 5 +++++ R/utils_spatial.R | 1 + inst/tinytest/test_cluster_correction.R | 3 ++- inst/tinytest/test_flux_wrapper.R | 1 + man/cluster_correction.Rd | 4 ++-- man/lonlat_to_utm.Rd | 3 +++ man/obs_flux.Rd | 3 ++- man/turbine_flights.Rd | 3 +-- man/turbine_flights_year.Rd | 6 +++--- vignettes/deterministic-example.Rmd | 9 +++++---- vignettes/simple-simulation-example.Rmd | 4 ++-- 15 files changed, 31 insertions(+), 19 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6721000..90bbe40 100644 --- a/DESCRIPTION +++ b/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")), diff --git a/NEWS.md b/NEWS.md index 5b7af7d..a68519a 100644 --- a/NEWS.md +++ b/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. diff --git a/R/cluster_correction.R b/R/cluster_correction.R index c5d3432..5421696 100644 --- a/R/cluster_correction.R +++ b/R/cluster_correction.R @@ -117,7 +117,7 @@ cluster_correction_a <- function(eff_detection_width, df_turbines, crs = 4326){ #' ##option 1 - provide the average distance between turbine neighbours #' #' sf_turbines <- sf::st_as_sf(df_turbines, -#' coords = c(lon_col, lat_col), +#' coords = c("lon", "lat"), #' crs = 4326) #' #' idx <- sf::st_nearest_feature(sf_turbines) diff --git a/R/flux_scaling.R b/R/flux_scaling.R index 42cb0ef..f808cae 100644 --- a/R/flux_scaling.R +++ b/R/flux_scaling.R @@ -17,7 +17,6 @@ #' [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. -#' TODO make better words #' @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 diff --git a/R/flux_wrapper.R b/R/flux_wrapper.R index 1f6dc62..46dec33 100644 --- a/R/flux_wrapper.R +++ b/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 diff --git a/R/utils_spatial.R b/R/utils_spatial.R index 9d03150..4a8d409 100644 --- a/R/utils_spatial.R +++ b/R/utils_spatial.R @@ -2,6 +2,7 @@ #' #' Equation from https://r.geocompx.org/reproj-geo-data#which-crs #' +#' @param lonlat a lon, lat coordinate pair #' lonlat_to_utm = function(lonlat) { diff --git a/inst/tinytest/test_cluster_correction.R b/inst/tinytest/test_cluster_correction.R index 5cef3fc..ea833ba 100644 --- a/inst/tinytest/test_cluster_correction.R +++ b/inst/tinytest/test_cluster_correction.R @@ -123,7 +123,7 @@ expect_error( ## Input validation: NA eff_detection_width -------------------------------- expect_error( cluster_correction_a(df_turbines = two_turbines, eff_detection_width = NA), - 'eff_detection_width must be a single numeric value' + 'Numeric input expected' ) # cluster_correction_l---------------------------------------------------------- @@ -204,3 +204,4 @@ expect_error( cluster_correction_l(avg_min_distance = 1000, eff_detection_width = NA), 'Numeric input expected' ) + diff --git a/inst/tinytest/test_flux_wrapper.R b/inst/tinytest/test_flux_wrapper.R index ce4c2cc..bdfe576 100644 --- a/inst/tinytest/test_flux_wrapper.R +++ b/inst/tinytest/test_flux_wrapper.R @@ -143,6 +143,7 @@ expect_equal(turbine_flights_year( time_units = "min", eff_detection_width, eff_detection_height, + 1, rotor_diameter, hub_height, prop_day, diff --git a/man/cluster_correction.Rd b/man/cluster_correction.Rd index 38f66e2..d6a6eba 100644 --- a/man/cluster_correction.Rd +++ b/man/cluster_correction.Rd @@ -9,7 +9,7 @@ cluster_correction_a(eff_detection_width, df_turbines, crs = 4326) cluster_correction_l( - eff_detection_width = NULL, + eff_detection_width, avg_min_distance = NULL, df_turbines = NULL, crs = 4326 @@ -86,7 +86,7 @@ df_turbines <- data.frame( ##option 1 - provide the average distance between turbine neighbours sf_turbines <- sf::st_as_sf(df_turbines, - coords = c(lon_col, lat_col), + coords = c("lon", "lat"), crs = 4326) idx <- sf::st_nearest_feature(sf_turbines) diff --git a/man/lonlat_to_utm.Rd b/man/lonlat_to_utm.Rd index 02f78fc..dc93bfd 100644 --- a/man/lonlat_to_utm.Rd +++ b/man/lonlat_to_utm.Rd @@ -6,6 +6,9 @@ \usage{ lonlat_to_utm(lonlat) } +\arguments{ +\item{lonlat}{a lon, lat coordinate pair} +} \description{ Equation from https://r.geocompx.org/reproj-geo-data#which-crs } diff --git a/man/obs_flux.Rd b/man/obs_flux.Rd index 9512b35..efe3218 100644 --- a/man/obs_flux.Rd +++ b/man/obs_flux.Rd @@ -11,7 +11,8 @@ obs_flux(encounter_rate, eff_detection_width, eff_detection_height) as output by \code{\link[=encounter_rate]{encounter_rate()}} or similar.} \item{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 \code{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 \code{eff_detection_height}.} \item{eff_detection_height}{numeric; The effective detection height. Because the density of flights typically decreases with increasing height, the assumption of uniform diff --git a/man/turbine_flights.Rd b/man/turbine_flights.Rd index d9eea90..12d82c4 100644 --- a/man/turbine_flights.Rd +++ b/man/turbine_flights.Rd @@ -14,8 +14,7 @@ per unit time as output by \code{\link[=obs_flux]{obs_flux()}} or similar.} Either the cluster correction factor as output by \code{\link[=cluster_correction_a]{cluster_correction_a()}}/ \code{\link[=cluster_correction_l]{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. -TODO make better words} +Defaults to 1.} \item{rotor_diameter}{numeric; rotor diameter. Must be in the equivalent units to the unit area of \code{obs_flux} (i.e., if the \code{obs_flux} is per m\eqn{^2}, \code{rotor_diameter} must be in m).} diff --git a/man/turbine_flights_year.Rd b/man/turbine_flights_year.Rd index 8baa004..d8e478e 100644 --- a/man/turbine_flights_year.Rd +++ b/man/turbine_flights_year.Rd @@ -27,7 +27,8 @@ as output by \code{\link[=encounter_rate]{encounter_rate()}} or similar.} \item{time_units}{Time units of \code{encounter_rate}. Defaults to "min" (i.e. flights per minute).} \item{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 \code{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 \code{eff_detection_height}.} \item{eff_detection_height}{numeric; The effective detection height. Because the density of flights typically decreases with increasing height, the assumption of uniform @@ -44,8 +45,7 @@ Must be in the same units as \code{eff_detection_width}.} Either the cluster correction factor as output by \code{\link[=cluster_correction_a]{cluster_correction_a()}}/ \code{\link[=cluster_correction_l]{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. -TODO make better words} +Defaults to 1 (no correction).} \item{rotor_diameter}{numeric; rotor diameter. Must be in the equivalent units to the unit area of \code{obs_flux} (i.e., if the \code{obs_flux} is per m\eqn{^2}, \code{rotor_diameter} must be in m).} diff --git a/vignettes/deterministic-example.Rmd b/vignettes/deterministic-example.Rmd index 6f397ad..e5c9e10 100644 --- a/vignettes/deterministic-example.Rmd +++ b/vignettes/deterministic-example.Rmd @@ -228,8 +228,9 @@ dt_survey <- copy(df_survey) setDT(dt_obs) setDT(dt_survey) -dt_obs_survey <- dt_obs[height <= max_rsh, .("size" = sum(size)), survey_id][dt_survey, - on = "survey_id"] +dt_obs_survey <- dt_obs[height <= max_rsh, .("size" = sum(size)), survey_id + ][dt_survey, + on = "survey_id"] # need sum(size) because we want one row per survey (not observation) ``` @@ -271,8 +272,8 @@ obs_flux_min <- obs_flux( ) effective_flux <- cluster_correction_a( - df_turbines = df_turbines, - eff_detection_width = 2*edr_from_distmodel(ds_model)) + eff_detection_width = 2*edr_from_distmodel(ds_model), + df_turbines = df_turbines) #flights through turbine / min diff --git a/vignettes/simple-simulation-example.Rmd b/vignettes/simple-simulation-example.Rmd index df51a2e..793960e 100644 --- a/vignettes/simple-simulation-example.Rmd +++ b/vignettes/simple-simulation-example.Rmd @@ -417,8 +417,8 @@ lst_results <- lapply(1:iterations, function(i) { # Step 1.5 calculate the cluster correction df_turbines_results$cluster_corr <- cluster_correction_a( - df_turbines = df_turbines_results, - eff_detection_width = 2*edr_i + eff_detection_width = 2*edr_i, + df_turbines = df_turbines_results ) # Step 2 calculate flux through turbine per year From f5fbad12a90eb0d3a95e9d5df158248db344c5bb Mon Sep 17 00:00:00 2001 From: sahinya17 Date: Thu, 23 Apr 2026 10:53:44 +1000 Subject: [PATCH 10/10] Fixed test case in cluster_correction_l function --- R/cluster_correction.R | 6 +++--- inst/tinytest/test_cluster_correction.R | 21 ++++++++++++++++++++- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/R/cluster_correction.R b/R/cluster_correction.R index 5421696..ebaab9f 100644 --- a/R/cluster_correction.R +++ b/R/cluster_correction.R @@ -152,11 +152,11 @@ cluster_correction_l <- function(eff_detection_width, if( is.null(avg_min_distance ) & is.null(df_turbines)){ - stop("Please enter one of eff_detection_width or avg_min_distance") + 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 eff_detection_width and avg_min_distance") + stop("Please enter only one of avg_min_distance and df_turbines") } diff --git a/inst/tinytest/test_cluster_correction.R b/inst/tinytest/test_cluster_correction.R index ea833ba..fccfa34 100644 --- a/inst/tinytest/test_cluster_correction.R +++ b/inst/tinytest/test_cluster_correction.R @@ -128,7 +128,7 @@ expect_error( # cluster_correction_l---------------------------------------------------------- -## Basic correctness: result equals eff_detection_width / avg_min_distance---- +## Basic correctness: result equals avg_min_distance / eff_detection_width---- expect_equal( cluster_correction_l(avg_min_distance = 1000, eff_detection_width = 500), 1 @@ -205,3 +205,22 @@ expect_error( 'Numeric input expected' ) +## Input validation: mutual exclusivity of avg_min_distance and df_turbines -- +expect_error( + cluster_correction_l(eff_detection_width = 500), + 'Please enter one of avg_min_distance or df_turbines' +) + +two_turbines_l <- data.frame( + lon = c(145.947140, 145.94500), + lat = c(-38.491772, -38.49100) +) +expect_error( + cluster_correction_l( + eff_detection_width = 500, + avg_min_distance = 1000, + df_turbines = two_turbines_l + ), + 'Please enter only one of avg_min_distance and df_turbines' +) +