diff --git a/DESCRIPTION b/DESCRIPTION
index e442c15..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")),
@@ -22,6 +22,7 @@ Imports:
   Distance,
   methods,
   Rdpack,
+  sf,
   units
 Suggests: 
   tinytest,
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/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
new file mode 100644
index 0000000..ebaab9f
--- /dev/null
+++ b/R/cluster_correction.R
@@ -0,0 +1,253 @@
+#' Correct for turbine clustering
+#' 
+#' The construction of the CRM requires that the flight flux through the observed 
+#' window be scaled to the size of a turbine and then applied to each turbine. If 
+#' If the (effective) detection range of the observer overlaps multiple turbines
+#' the flux needs to be scaled appropriate over the turbines in the observation
+#' range.
+#' 
+#' @details
+#' If the turbine layout is more compact and/or your effective detection width is large, then the effective detection area contains more than one turbine, in which case the observed flux needs to be weighted to account for the effective turbines in the observed area. This accounts for the fact that the observed flux is an instantaneous measure and cannot interact with multiple turbines at once.
+#' 
+#' This package provides a 1D line correction [cluster_correction_l()] and 2D aerial [cluster_correction_a()].
+#' 
+#' ##  [cluster_correction_l()]
+#' 
+#' The 1D (line) correction calculates the ratio of the effective strip width (ESW) and the mean distance between each turbine and it's nearest neighbour. It's suitable for line transects or for interim CRM early in the project life-cycle where the turbine spacing can be estimated but the layout may not be available.
+#' 
+#' ## [cluster_correction_a()]
+#' 
+#' The 2D aerial correction places a buffer of radius equal to the effective detection radius (EDR) around each turbine, then generates a spatial union of those buffers. The correction is the ratio of the actual turbine area to the area if the turbine buffers do not overlap.
+#'  
+#' ## Relationship to other spatial corrections
+#' 
+#' This correction is derived from different assumptions but fulfils a similar model correction as the \code{sqrt(n_turbines)} in previous CRM  \insertCite{Smales2013,Band2012}{collision}. 
+#' 
+#' Alternatively, a spatial probability surface can be generated by the analyst using kernel density estimate (KDE) or similar and used instead to estimate the proportion of flights at a given turbine estimation.  That approach requires more data and site coverage.
+#' 
+#'
+#' 
+#'    
+#'    
+#' @name cluster_correction    
+
+
+
+#' @rdname cluster_correction
+#'
+#' @param eff_detection_width numeric; the total effective detection width which is 
+#'     2 x effective detection radius or 2 x effective strip (half) width. 
+#' @param df_turbines data.frame; a data.frame with one row per turbine, 
+#'    containing the coordinates of each turbine as latitude, longitude pairs. 
+#'    Coordinate columns can be named (lon, long, longitude, x) and (lat, latitude, y)
+#' @param crs optional coordinate system epsg code for the coordinates.  Defaults to crs = 4326 (WGS 84 - lat/lon). 
+#' 
+#' 
+#' @return numeric; the minimum of 1 or a fractional scaling factor, calculated as the ratio of the total area covered by the turbines buffered by the EDA divided by the area if each turbine buffer were completely independent.
+#' 
+#' @examples
+#' 
+#' # generate inputs information
+#' 
+#' df_turbines <- data.frame(
+#'   turbine_id = c("T1", "T2", "T3"),
+#'   lat = c(-35.036287, -35.03900, -35.03200),
+#'   lon = c(145.947140, 145.94500, 145.94200)
+#' )
+#' 
+#' eff_detection_width <- 500 
+#' 
+#' cluster_correction_a(eff_detection_width = 500, df_turbines, crs = 4326)
+#' 
+#' @export
+cluster_correction_a <- function(eff_detection_width, df_turbines, crs = 4326){
+  
+  stopifnot('eff_detection_width must be a single numeric value' = 
+              length(eff_detection_width) == 1)
+  check_num_bounds(eff_detection_width, min = 1, max = Inf)
+  
+  
+  sf_turbines <- df_to_projsf(df_turbines, crs) # internal helper function
+  
+  ## total area with no overlap
+  area_separate <- nrow(df_turbines)* pi*(eff_detection_width*0.5)^2  #half width==radius
+  
+  ## area buffered by EDA
+  sf_turb_buffer <- sf::st_make_valid(
+    sf::st_buffer(x = sf_turbines,
+                  dist = eff_detection_width*0.5,
+                  nQuadSegs = 360))
+  
+  ## union
+  sf_turb_union <- sf::st_make_valid(sf::st_union(sf_turb_buffer))
+  
+  ## Area
+  turb_area_buf <- sf::st_area(sf_turb_union) |> as.numeric()
+  
+  ## return
+  turbine_prob <- min( turb_area_buf/area_separate, 1.0 ) 
+  
+  return(turbine_prob)
+}
+
+
+#' @rdname cluster_correction
+#'
+#' @param avg_min_distance numeric; (optional) the mean distance between each turbine and its nearest neighbour. Alternatively you can input a data.frame with turbine coordinates
+#' @param eff_detection_width numeric; the effective detection (half) width or effective detection radius
+#' @param df_turbines data.frame; (optional) a data.frame with one row per turbine, 
+#'    containing the coordinates of each turbine as latitude, longitude pairs. 
+#'    Coordinate columns can be named (lon, long, longitude, x) and (lat, latitude, y)
+#' @param crs optional coordinate system epsg code for the coordinates.  Only used if df_turbines provided. Defaults to crs = 4326 (WGS 84 - lat/lon). 
+#'
+#'
+#' @return numeric; the minimum of 1 or a fractional scaling factor, calculated as the ratio of the average distance between turbine neighbours and the effective strip width
+#' 
+#' @examples
+#' 
+#' # generate input data
+#' 
+#' 
+#' df_turbines <- data.frame(
+#'   turbine_id = c("T1", "T2", "T3"),
+#'   lat = c(-35.036287, -35.03900, -35.03200),
+#'   lon = c(145.947140, 145.94500, 145.94200)
+#' )
+#'
+#' ##option 1 - provide the average distance between turbine neighbours
+#' 
+#' sf_turbines <- sf::st_as_sf(df_turbines, 
+#'   coords = c("lon", "lat"), 
+#'   crs = 4326)
+#'   
+#' idx <- sf::st_nearest_feature(sf_turbines)
+#' dist_avg <- sf::st_distance(sf_turbines, sf_turbines[idx,], by_element=TRUE) |>
+#'    mean() |>
+#'    as.numeric()
+#' 
+#' 
+#' corr1 <- cluster_correction_l(eff_detection_width = 800, 
+#'                               avg_min_distance = dist_avg)
+#'    
+#' ## option 2 - provide the data frame (might be slower)
+#'  
+#' corr2 <- cluster_correction_l(eff_detection_width = 800, 
+#'                               df_turbines = df_turbines,
+#'                               crs = 4326 # provide crs with turbines
+#'                               )
+#'  ## should be equal (or very very close)
+#'  corr1; corr2
+#' 
+#' @export
+cluster_correction_l <- function(eff_detection_width,
+                                 avg_min_distance = NULL,
+                                 df_turbines = NULL,
+                                 crs = 4326){
+
+  stopifnot('eff_detection_width must be a single numeric value' = 
+              length(eff_detection_width) == 1)
+
+
+  check_num_bounds(eff_detection_width, min = 1, max = Inf)
+  
+  
+  if( is.null(avg_min_distance ) & is.null(df_turbines)){
+    stop("Please enter one of avg_min_distance or df_turbines")
+  }
+
+  if( !is.null(avg_min_distance) & !is.null(df_turbines)){
+    stop("Please enter only one of avg_min_distance and df_turbines")
+  }
+    
+  
+  if(!is.null(avg_min_distance)){
+    stopifnot('avg_min_distance must be a single numeric value' = 
+                length(avg_min_distance) == 1)
+    check_num_bounds(avg_min_distance, min = 1, max = Inf)
+  }
+  
+  if(!is.null(df_turbines)){
+    
+    sf_turbines <- df_to_projsf(df_turbines, crs) # internal helper function
+    
+    #index of nearest neighbour
+    idx <- sf::st_nearest_feature(sf_turbines)
+    
+    #calculate average distance
+    avg_min_distance <- sf::st_distance(sf_turbines, 
+                                        sf_turbines[idx,], 
+                                        by_element=TRUE) |>
+      mean() |>
+      as.numeric()
+  }
+  
+    
+  return(min(avg_min_distance/eff_detection_width, 1.0)) 
+}
+
+
+# Convert df to sf
+# 
+# helper function shared by the cluster correction functions
+# 
+# @returns sf object in projected coords
+df_to_projsf <- function(df_turbines, crs){
+  
+  ## checks
+  stopifnot(
+    "df_turbines must be a data.frame" =
+      is.data.frame(df_turbines)
+  )
+  stopifnot(
+    "df_turbines must have at least one row" =
+      nrow(df_turbines) >= 1
+  )
+  
+  ## Grab coordinate data and check
+  lon_col <- names(df_turbines)[
+    tolower(names(df_turbines)) %in% c("lon",
+                                       "longitude",
+                                       "long",
+                                       "x",
+                                       "easting")
+  ]
+  
+  lat_col <- names(df_turbines)[
+    tolower(names(df_turbines)) %in% c("lat",
+                                       "latitude",
+                                       "y",
+                                       "northing")
+  ]
+  
+  stopifnot(
+    "df_turbines must contain one longitude column named lon, long, longitude, x, or easting" =
+      length(lon_col) == 1
+  )
+  stopifnot(
+    "df_turbines must contain one latitude column named lat, latitude, y, or northing" =
+      length(lat_col) == 1
+  )
+  
+  sf_turbines <- sf::st_as_sf(df_turbines, 
+                              coords = c(lon_col, lat_col), 
+                              crs = crs)
+
+  
+  # Convert to projected coord system for calc in metres
+  if(sf::st_is_longlat(sf_turbines)){
+    
+    check_num_bounds(sf::st_coordinates(sf_turbines)[, "X"], -180, 180 )
+    check_num_bounds(sf::st_coordinates(sf_turbines)[, "Y"], -90, 90 )
+    
+    projcrs <- lonlat_to_utm(sf::st_coordinates(sf_turbines))
+    sf_turbines <- sf::st_transform(sf_turbines, crs = projcrs)
+    
+  }
+    
+  check_num_bounds(sf::st_coordinates(sf_turbines)[, "X"], 150000, 1e6 )
+  check_num_bounds(abs(sf::st_coordinates(sf_turbines)[, "Y"]), 1.e6, 10.e6  )
+  
+  
+  return(sf_turbines)
+  
+}
\ No newline at end of file
diff --git a/R/flux_scaling.R b/R/flux_scaling.R
index 22abb25..f808cae 100644
--- a/R/flux_scaling.R
+++ b/R/flux_scaling.R
@@ -11,7 +11,12 @@
 #' output flights / min).
 #' 
 #' @param obs_flux numeric; number of flights through one unit area of vertical space
-#'    per unit time as output by [obs_flux()] or similar
+#'    per unit time as output by [obs_flux()] or similar.
+#' @param spatial_correction numeric; factor to correct flux for spatial considerations. 
+#'    Either the cluster correction factor as output by [cluster_correction_a()]/ 
+#'    [cluster_correction_l()] or similar (assuming flat activity across the site) 
+#'    or a correction based on the spatial distribution of flights on the site. 
+#'    Defaults to 1.
 #' @param rotor_diameter numeric; rotor diameter. Must be in the equivalent units to 
 #'    the unit area of `obs_flux` (i.e., if the `obs_flux` is per m\eqn{^2}, `rotor_diameter` must be in m).
 #' @param hub_height numeric; hub height. Must be in the equivalent units to 
@@ -26,13 +31,17 @@
 #' @example examples/flux_example.R
 #'
 #' @export
-turbine_flights <- function(obs_flux, rotor_diameter, hub_height) {
+turbine_flights <- function(obs_flux, 
+                            spatial_correction = 1, 
+                            rotor_diameter, 
+                            hub_height) {
   check_num_bounds(obs_flux, min = 0)
+  check_num_bounds(spatial_correction, min = 0)
   check_num_bounds(rotor_diameter, min = 0)
   check_num_bounds(hub_height, min = 0)
   
   turbine_height <- hub_height + 0.5*rotor_diameter
-  res <- obs_flux * rotor_diameter * turbine_height
+  res <- obs_flux * spatial_correction * rotor_diameter * turbine_height
   return(res)
 }
 
diff --git a/R/flux_wrapper.R b/R/flux_wrapper.R
index 0fe08e0..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
@@ -84,6 +89,7 @@ turbine_flights_year <- function(
     time_units = "min",
     eff_detection_width,
     eff_detection_height,
+    spatial_correction = 1,
     rotor_diameter,
     hub_height,
     prop_day,
@@ -106,6 +112,7 @@ turbine_flights_year <- function(
 
   flights_min <- turbine_flights(
     obs_flux = obs_flux,
+    spatial_correction = spatial_correction,
     rotor_diameter = rotor_diameter,
     hub_height = hub_height
   )
diff --git a/R/observed_flux.R b/R/observed_flux.R
index 08cf3d5..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
@@ -84,4 +85,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/R/scratch.R b/R/scratch.R
index e69de29..5b37de1 100644
--- a/R/scratch.R
+++ b/R/scratch.R
@@ -0,0 +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
+#' 
+#' 
+#' 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
+#'   }
+#'   
+#' }
+#' 
+#' 
+#' 
+#' 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..4a8d409
--- /dev/null
+++ b/R/utils_spatial.R
@@ -0,0 +1,16 @@
+#' lon lat to UTM
+#' 
+#' Equation from https://r.geocompx.org/reproj-geo-data#which-crs
+#' 
+#' @param lonlat a lon, lat coordinate pair
+#' 
+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
new file mode 100644
index 0000000..fccfa34
--- /dev/null
+++ b/inst/tinytest/test_cluster_correction.R
@@ -0,0 +1,226 @@
+# Tests for cluster_correction_a and cluster_correction_l
+
+# Helper df
+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(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(eff_detection_width = 500, one_turbine),
+  1,
+  tolerance = 0.05
+)
+
+## Very widely spaced turbines: factor close to 1 ---------------------------
+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(df_turbines = 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(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(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(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_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(
+  lon = c(0, 0.01),
+  lat = c(51, 51),
+  latitude = c(51, 51)
+)
+expect_error(
+  cluster_correction_a(500, df_multi_lat)
+)
+
+## Input validation: non-data.frame -----------------------------------------
+expect_error(
+  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(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(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(500, data.frame(lon = c(0, 1), lati = c(51, 51)))
+)
+
+## Input validation: NA coordinates ----------------------------------------
+expect_error(
+  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(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(500, data.frame(lon = c(0, -10), lat = c(51, 51))),
+)
+
+## Input validation: eff_detection_width type/length -----------------------
+expect_error(cluster_correction_a(df_turbines = two_turbines, eff_detection_width = '500'))
+expect_error(
+  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(df_turbines = two_turbines, eff_detection_width = 0),
+  'variable out of bounds'
+)
+expect_error(
+  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(df_turbines = two_turbines, eff_detection_width = NA),
+  'Numeric input expected'
+)
+
+# cluster_correction_l----------------------------------------------------------
+
+## Basic correctness: result equals avg_min_distance / eff_detection_width----
+expect_equal(
+  cluster_correction_l(avg_min_distance = 1000, eff_detection_width = 500),
+  1
+)
+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),
+  0.4
+)
+
+## Returns single numeric ---------------------------------------------------
+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),
+  'Numeric input expected'
+)
+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),
+  'variable out of bounds'
+)
+expect_error(
+  cluster_correction_l(avg_min_distance = -100, eff_detection_width = 500),
+  'variable out of bounds'
+)
+
+## Input validation: NA avg_min_distance ------------------------------------
+expect_error(
+  cluster_correction_l(avg_min_distance = NA, eff_detection_width = 500),
+  'Numeric input expected'
+)
+
+## Input validation: eff_detection_width type/length ------------------------
+expect_error(
+  cluster_correction_l(avg_min_distance = 1000, eff_detection_width = '500'),
+  'Numeric input expected'
+)
+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),
+  'variable out of bounds'
+)
+expect_error(
+  cluster_correction_l(avg_min_distance = 1000, eff_detection_width = -500),
+  'variable out of bounds'
+)
+
+## Input validation: NA eff_detection_width ---------------------------------
+expect_error(
+  cluster_correction_l(avg_min_distance = 1000, eff_detection_width = NA),
+  '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'
+)
+
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----------------
diff --git a/inst/tinytest/test_flux_wrapper.R b/inst/tinytest/test_flux_wrapper.R
index bf7ec6c..bdfe576 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."
@@ -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
new file mode 100644
index 0000000..d6a6eba
--- /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,
+  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", "lat"), 
+  crs = 4326)
+  
+idx <- sf::st_nearest_feature(sf_turbines)
+dist_avg <- sf::st_distance(sf_turbines, sf_turbines[idx,], by_element=TRUE) |>
+   mean() |>
+   as.numeric()
+
+
+corr1 <- cluster_correction_l(eff_detection_width = 800, 
+                              avg_min_distance = dist_avg)
+   
+## option 2 - provide the data frame (might be slower)
+ 
+corr2 <- cluster_correction_l(eff_detection_width = 800, 
+                              df_turbines = df_turbines,
+                              crs = 4326 # provide crs with turbines
+                              )
+ ## should be equal (or very very close)
+ corr1; corr2
+
+}
diff --git a/man/lonlat_to_utm.Rd b/man/lonlat_to_utm.Rd
new file mode 100644
index 0000000..dc93bfd
--- /dev/null
+++ b/man/lonlat_to_utm.Rd
@@ -0,0 +1,14 @@
+% 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)
+}
+\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 ea89a56..12d82c4 100644
--- a/man/turbine_flights.Rd
+++ b/man/turbine_flights.Rd
@@ -4,11 +4,17 @@
 \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.}
 
 \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..d8e478e 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,
@@ -26,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
@@ -39,6 +41,12 @@ 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 (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 a05a72f..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)
 ```
 
@@ -251,13 +252,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 +271,15 @@ obs_flux_min <- obs_flux(
   eff_detection_height = max_rsh
 )
 
+effective_flux <- cluster_correction_a(
+  eff_detection_width = 2*edr_from_distmodel(ds_model),
+  df_turbines = df_turbines)
+
 #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 +289,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 +303,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..793960e 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(
+    eff_detection_width = 2*edr_i,
+    df_turbines = df_turbines_results
+  )
+  
   # 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}