"""Spatial metrics for generated wind-farm layouts.
Each metric accepts a :class:`~oowl.layouts.base.Farm` or an ``(N, 2)`` array of
turbine positions, ``N`` zero included, and a metric undefined for so few
turbines is ``nan``. Area-based metrics use ``Farm.boundary`` when available and
fall back to the convex hull of the turbines. Pass an explicit ``boundary``, a
:class:`~oowl.boundaries.base.Boundary` or a Shapely polygon, to override, or
``boundary=False`` to force the hull fallback.
"""
from __future__ import annotations
from collections.abc import Iterable, Sequence
import numpy as np
import shapely
from shapely import contains_xy
from shapely.geometry.base import BaseGeometry
from ._validation import _validate_nonnegative, _validate_positive, _validate_real_array
#: Most hexagons :func:`hexagon_centers` lays out over a polygon's box.
_MAX_HEXAGONS = 10_000_000
__all__ = [
"nearest_neighbor_distances",
"n_turbines",
"nearest_neighbor_mean",
"nearest_neighbor_variation",
"hull_area",
"region_area",
"density",
"clark_evans",
"aspect_ratio",
"bounding_box_diagonal",
"ripley_k",
"ripley_l",
"quadrat_window",
"hexagon_vertices",
"hexagon_centers",
"quadrat_counts",
"quadrat_p_value",
"SCALAR_METRICS",
"DEFAULT_SUMMARY_METRICS",
"summary",
]
def _as_positions(obj: object) -> np.ndarray:
"""Validated ``(N, 2)`` finite float positions from a Farm or array.
A layout without turbines is valid, and a metric undefined for it is ``nan``.
"""
positions = obj.positions if hasattr(obj, "positions") else obj
positions = _validate_real_array("positions", positions)
if positions.size == 0:
positions = positions.reshape(0, 2)
if positions.ndim != 2 or positions.shape[1] != 2:
raise ValueError("positions must have shape (N, 2)")
if not np.all(np.isfinite(positions)):
raise ValueError("positions must contain only finite values")
return positions
def _as_radii(radii: Sequence[float] | np.ndarray) -> np.ndarray:
"""Radii in meters as a 1-D finite float array, shape ``(len(radii),)``."""
values = _validate_real_array("radii", radii)
if values.ndim != 1:
raise ValueError("radii must have shape (M,)")
if not np.all(np.isfinite(values)):
raise ValueError("radii must contain only finite values")
return values
def _resolve_boundary(obj, boundary):
"""Effective boundary, explicit argument first, else ``Farm.boundary``.
``boundary=False`` returns ``None`` to force the convex-hull fallback.
"""
from .boundaries.base import Boundary
if boundary is False:
return None
if boundary is not None:
if not isinstance(boundary, (Boundary, BaseGeometry)):
raise TypeError(
"boundary must be a Boundary, a Shapely geometry, None or False, "
f"not {type(boundary).__name__}"
)
return boundary
return getattr(obj, "boundary", None)
def _boundary_geometry(boundary):
"""Shapely geometry for a Boundary, or the geometry itself, or ``None``."""
if boundary is None:
return None
return getattr(boundary, "geometry", boundary)
def _region_geometry(obj, boundary) -> BaseGeometry:
"""Study-region geometry, the resolved boundary or else the position hull."""
geometry = _boundary_geometry(_resolve_boundary(obj, boundary))
if geometry is None:
geometry = shapely.multipoints(_as_positions(obj)).convex_hull
return geometry
[docs]
def n_turbines(obj: object) -> int:
"""Number of turbines in the layout, from a Farm or an ``(N, 2)`` array."""
return int(len(_as_positions(obj)))
[docs]
def nearest_neighbor_distances(obj: object) -> np.ndarray:
"""Nearest-neighbor distance in meters per turbine, shape ``(N,)``.
``nan`` for a single turbine, which has no neighbor.
"""
layout = _as_positions(obj)
if len(layout) < 2:
return np.full(len(layout), np.nan)
from .layouts.base import _nearest_neighbor_distances
return _nearest_neighbor_distances(layout)
[docs]
def nearest_neighbor_mean(obj: object) -> float:
"""Mean nearest-neighbor distance in meters between turbine positions.
``nan`` for fewer than two turbines.
"""
nnd = nearest_neighbor_distances(obj)
return float(nnd.mean()) if len(nnd) else float("nan")
[docs]
def nearest_neighbor_variation(obj: object) -> float:
"""Coefficient of variation of the nearest-neighbor distances.
Near 0 for a regular layout, larger for irregular or clustered ones. ``nan``
when the mean distance is zero or undefined.
"""
nnd = nearest_neighbor_distances(obj)
mean = float(nnd.mean()) if len(nnd) else float("nan")
if not np.isfinite(mean) or mean == 0.0:
return float("nan")
return float(nnd.std() / mean)
[docs]
def hull_area(obj: object) -> float:
"""Convex-hull area in square meters of the positions, ignoring any boundary.
``nan`` when the hull encloses no area.
"""
layout = _as_positions(obj)
if len(layout) < 3:
return float("nan")
from scipy.spatial import ConvexHull, QhullError
try:
return float(ConvexHull(layout).volume)
except QhullError:
return float("nan")
[docs]
def region_area(obj: object, boundary: object = None) -> float:
"""Study-region area in square meters, from the boundary or else the hull."""
resolved = _resolve_boundary(obj, boundary)
if resolved is not None:
area = getattr(resolved, "area", None)
if area is None:
area = _boundary_geometry(resolved).area
return float(area)
return hull_area(obj)
[docs]
def density(obj: object, boundary: object = None) -> float:
"""Turbines per square meter over the study region."""
a = region_area(obj, boundary)
return len(_as_positions(obj)) / a if a > 0 else float("nan")
[docs]
def clark_evans(obj: object, boundary: object = None) -> float:
"""Clark-Evans ratio R of the mean nearest-neighbor distance to its CSR value.
Under complete spatial randomness (CSR) the expected mean is
``0.5 * sqrt(A / N)`` with ``A`` the study-region area. ``R > 1`` indicates
regularity, ``R < 1`` clustering and ``R ~ 1`` CSR. ``nan`` when the area is
not positive or there are no turbines.
See Clark and Evans (1954), Ecology 35(4), https://doi.org/10.2307/1931034.
"""
a = region_area(obj, boundary)
layout = _as_positions(obj)
if not a > 0 or not len(layout):
return float("nan")
expected = 0.5 * np.sqrt(a / len(layout))
return float(nearest_neighbor_mean(layout) / expected)
[docs]
def aspect_ratio(obj: object) -> float:
"""Square root of the largest over the smallest position covariance eigenvalue.
This is the ratio of the spreads along the principal axes, 1 for an
isotropic layout and large for an elongated one. ``nan`` for a single
turbine or a degenerate spread, such as a straight row.
"""
layout = _as_positions(obj)
if len(layout) < 2:
return float("nan")
cov = np.cov(layout.T)
eig = np.linalg.eigvalsh(cov)
# Rounding in the covariance sums leaves collinear positions a smallest
# eigenvalue of order N eps times the largest, which is not a real spread.
if eig.min() <= len(layout) * np.finfo(float).eps * eig.max():
return float("nan")
return float(np.sqrt(eig.max() / eig.min()))
[docs]
def bounding_box_diagonal(obj: object) -> float:
"""Diagonal length in meters of the axis-aligned bounding box of the positions.
``nan`` for a layout without turbines.
"""
layout = _as_positions(obj)
if not len(layout):
return float("nan")
dx = layout[:, 0].max() - layout[:, 0].min()
dy = layout[:, 1].max() - layout[:, 1].min()
return float(np.sqrt(dx**2 + dy**2))
[docs]
def ripley_k(
obj: object,
radii: Sequence[float] | np.ndarray,
boundary: object = None,
) -> np.ndarray:
"""Ripley's K(r) in square meters at the given radii in meters.
``K(r) = A / N^2 * sum_{i != j} 1[d_ij <= r]`` with ``A`` the study-region
area. No edge correction is applied, so with the bounding box of the
positions as ``boundary`` it matches ``pointpats.distance_statistics.k``.
``radii`` must be one-dimensional. Returns shape ``(len(radii),)``, all
``nan`` for fewer than two turbines.
See Ripley (1977), Journal of the Royal Statistical Society B 39(2),
https://doi.org/10.1111/j.2517-6161.1977.tb01615.x.
"""
layout = _as_positions(obj)
radii = _as_radii(radii)
n = len(layout)
a = region_area(obj, boundary)
if n < 2 or not a > 0:
return np.full(radii.shape, np.nan)
from scipy.spatial.distance import pdist
pair_distances = np.sort(pdist(layout))
counts = np.searchsorted(pair_distances, radii, side="right")
return a * (2.0 * counts) / (n**2)
[docs]
def ripley_l(
obj: object,
radii: Sequence[float] | np.ndarray,
boundary: object = None,
) -> np.ndarray:
"""Ripley's L(r) = sqrt(K(r) / pi) in meters, shape ``(len(radii),)``.
``L(r) - r > 0`` indicates clustering and ``< 0`` regularity.
See https://doi.org/10.1111/j.2517-6161.1977.tb01615.x.
"""
return np.sqrt(ripley_k(obj, radii, boundary) / np.pi)
[docs]
def quadrat_window(
obj: object,
boundary: object = None,
hexagon_fraction: float = 0.15,
buffer_factor: float = 1.0,
) -> tuple[BaseGeometry, float]:
"""Study-region polygon and hexagon side length for the quadrat test.
The hexagon side length is ``hexagon_fraction * sqrt(region area)`` in
meters, so the grid depends on the region and not on turbine count. The
region is buffered outward by ``buffer_factor`` side lengths so the lattice
covers every turbine. ``hexagon_fraction`` must be finite and > 0, and
``buffer_factor`` finite and >= 0. Returns ``(polygon, side_length)``. A
region without area, such as the hull of collinear turbines, gives a ``nan``
or zero side length and is not buffered.
"""
hexagon_fraction = _validate_positive("hexagon_fraction", hexagon_fraction)
buffer_factor = _validate_nonnegative("buffer_factor", buffer_factor)
geometry = _region_geometry(obj, boundary)
side = hexagon_fraction * np.sqrt(region_area(obj, boundary))
poly = (
geometry.buffer(buffer_factor * side)
if buffer_factor and side > 0
else geometry
)
return poly, side
[docs]
def hexagon_vertices(
center_x: float, center_y: float, side_length: float
) -> np.ndarray:
"""Closed vertex ring ``(7, 2)`` of a flat-top hexagon centered at the point."""
cx, cy, side = center_x, center_y, side_length
sh = side * np.cos(np.pi / 6.0) # apothem (center-to-edge)
return np.array(
[
[cx + side, cy],
[cx + side / 2.0, cy + sh],
[cx - side / 2.0, cy + sh],
[cx - side, cy],
[cx - side / 2.0, cy - sh],
[cx + side / 2.0, cy - sh],
[cx + side, cy],
]
)
[docs]
def hexagon_centers(polygon: BaseGeometry, side_length: float) -> np.ndarray:
"""Centers ``(M, 2)`` of a flat-top hexagon lattice inside ``polygon``.
The cells tessellate, so nearest-center assignment equals hexagon membership.
Raises ``TypeError`` for a ``polygon`` that is no Shapely geometry, and
``ValueError`` for a ``side_length`` that is not finite and positive or
that would lay out more than ten million hexagons over the polygon's box.
"""
if not isinstance(polygon, BaseGeometry):
raise TypeError(
f"polygon must be a Shapely geometry, not {type(polygon).__name__}"
)
poly, side = polygon, _validate_positive("side_length", side_length)
minx, miny, maxx, maxy = poly.bounds
cells = (maxx - minx + side) * (maxy - miny + side) / (1.5 * np.sqrt(3.0) * side**2)
if cells > _MAX_HEXAGONS:
raise ValueError(
f"side_length {side:.4g} would lay out about {cells:.3g} hexagons "
"over the polygon, so choose a larger one"
)
sh = side * np.cos(np.pi / 6.0)
col_x = np.arange(minx, maxx + side, 1.5 * side)
centers = []
for col, cx in enumerate(col_x):
offset = sh if col % 2 == 1 else 0.0
cy = np.arange(miny + offset, maxy + sh, 2.0 * sh)
centers.append(np.column_stack([np.full(cy.shape, cx), cy]))
centers = np.vstack(centers)
inside = contains_xy(poly, centers[:, 0], centers[:, 1])
return centers[inside]
[docs]
def quadrat_counts(obj: object, centers: np.ndarray) -> np.ndarray:
"""Point count per hexagon, shape ``(len(centers),)``, by nearest cell center."""
layout = _as_positions(obj)
from scipy.spatial.distance import cdist
nearest = cdist(layout, centers).argmin(axis=1)
return np.bincount(nearest, minlength=len(centers)).astype(float)
def _cell_areas(centers: np.ndarray, side: float, region: BaseGeometry) -> np.ndarray:
"""Area of each flat-top hexagon cell inside ``region``, shape ``(len(centers),)``.
Cells that ``region`` holds properly keep their full area without an overlay.
"""
sh = side * np.cos(np.pi / 6.0)
cx = centers[:, 0, None]
cy = centers[:, 1, None]
dx = np.array([side, side / 2.0, -side / 2.0, -side, -side / 2.0, side / 2.0, side])
dy = np.array([0.0, sh, sh, 0.0, -sh, -sh, 0.0])
cells = shapely.polygons(np.stack([cx + dx, cy + dy], axis=-1))
shapely.prepare(region)
inside = shapely.contains_properly(region, cells)
area = shapely.area(cells)
outside = ~inside
area[outside] = shapely.area(shapely.intersection(cells[outside], region))
return area
[docs]
def quadrat_p_value(
obj: object,
boundary: object = None,
hexagon_fraction: float = 0.15,
buffer_factor: float = 1.0,
) -> float:
"""Chi-square CSR p-value from hexagonal quadrat counts.
Cells come from :func:`quadrat_window` and :func:`hexagon_centers`.
Each cell's expected count is proportional to its area inside the study
region, and cells with no such area are dropped.
A low value rejects complete spatial randomness, as under clustering. A high
value only fails to reject it and is not evidence of regularity. ``nan``
when fewer than two usable cells remain or no turbine lies in them.
"""
poly, side = quadrat_window(obj, boundary, hexagon_fraction, buffer_factor)
if not side > 0:
return float("nan")
centers = hexagon_centers(poly, side)
if len(centers) < 2:
return float("nan")
region = _region_geometry(obj, boundary)
cell_area = _cell_areas(centers, side, region)
usable = cell_area > 0.0
if usable.sum() < 2:
return float("nan")
counts = quadrat_counts(obj, centers)[usable]
# Turbines outside the study region leave no observation to test.
if not counts.sum() > 0:
return float("nan")
cell_area = cell_area[usable]
expected = counts.sum() * cell_area / cell_area.sum()
from scipy.stats import chisquare
return float(chisquare(counts, f_exp=expected).pvalue)
SCALAR_METRICS = {
"n_turbines": n_turbines,
"nearest_neighbor_mean": nearest_neighbor_mean,
"nearest_neighbor_variation": nearest_neighbor_variation,
"density": density,
"region_area": region_area,
"clark_evans": clark_evans,
"aspect_ratio": aspect_ratio,
"bounding_box_diagonal": bounding_box_diagonal,
"quadrat_p_value": quadrat_p_value,
"hull_area": hull_area,
}
"""Registry mapping public scalar metric names to their implementations.
Each entry takes a single Farm or ``(N, 2)`` array, and boundary-aware metrics
resolve the Farm's boundary automatically. Ripley's K and L take ``radii``, so
they are not scalar and are not included.
"""
DEFAULT_SUMMARY_METRICS = (
"n_turbines",
"nearest_neighbor_mean",
"nearest_neighbor_variation",
"density",
"region_area",
)
"""Metric names included by default in :func:`summary`."""
def _metric_names(metrics: Iterable[str] | None) -> tuple[str, ...]:
"""The requested metric names, :data:`DEFAULT_SUMMARY_METRICS` for ``None``."""
if isinstance(metrics, (str, bytes)):
raise TypeError(
"metrics must be a collection of metric names, not a single "
f"{type(metrics).__name__}"
)
names = tuple(DEFAULT_SUMMARY_METRICS) if metrics is None else tuple(metrics)
for name in names:
if not isinstance(name, str):
raise TypeError(f"metric names must be strings, not {type(name).__name__}")
return names
[docs]
def summary(obj: object, metrics: Iterable[str] | None = None) -> dict[str, float]:
"""Ordered ``name -> value`` dict of scalar metrics for a farm or array.
``metrics`` names entries of :data:`SCALAR_METRICS`, defaulting to
:data:`DEFAULT_SUMMARY_METRICS`. A string instead of a collection of names,
or a name that is not a string, raises ``TypeError``. Unknown names raise
``ValueError``.
"""
names = _metric_names(metrics)
result = {}
for name in names:
if name not in SCALAR_METRICS:
raise ValueError(
f"unknown metric {name!r}; available: {sorted(SCALAR_METRICS)}"
)
result[name] = SCALAR_METRICS[name](obj)
return result
