"""Nested boundary-contour turbine layout."""
from __future__ import annotations
from functools import lru_cache
from typing import TYPE_CHECKING
import numpy as np
import shapely
from .._validation import _validate_positive
from .base import _grouped_by_piece
from .candidate import _CandidateLayout
from .fitters import SelectedContourFitter
if TYPE_CHECKING: # pragma: no cover
from ..boundaries.base import Boundary
from .fitters import SelectedCandidateFitter
__all__ = ["Contour"]
def _contour_lines(geometry):
"""Every exterior and hole ring of a polygon or multipolygon."""
polygons = geometry.geoms if geometry.geom_type == "MultiPolygon" else [geometry]
return [
ring for polygon in polygons for ring in (polygon.exterior, *polygon.interiors)
]
#: Relative to a boundary's larger side, the tolerance its deepest point is
#: found to.
_DEPTH_TOLERANCE = 1e-3
#: Fraction of the nominal spacing below which two candidates are too close.
_MIN_SEPARATION = 0.6
def _along(line, distances: np.ndarray) -> np.ndarray:
"""Points ``distances`` along ``line``, with a last axis for x and y.
The same points as ``shapely.line_interpolate_point`` to the last bit, as
GEOS finds them by summing segment lengths in order. That walked the
line from its start once per point, which on a smooth farm's insets of
hundreds of vertices took most of the time a Contour field took to build.
"""
coordinates = shapely.get_coordinates(line)
dx = coordinates[1:, 0] - coordinates[:-1, 0]
dy = coordinates[1:, 1] - coordinates[:-1, 1]
lengths = np.sqrt(dx * dx + dy * dy)
reached = np.cumsum(lengths)
before = np.concatenate([[0.0], reached[:-1]])
segment = np.minimum(
np.searchsorted(reached, distances, side="right"), len(lengths) - 1
)
# A distance past the end lands on the last segment, which may have no
# length, and is replaced by the end below.
with np.errstate(divide="ignore", invalid="ignore"):
fraction = (distances - before[segment]) / lengths[segment]
start = coordinates[segment]
stop = coordinates[segment + 1]
point = np.stack([dx[segment], dy[segment]], axis=-1) * fraction[..., None]
point += start
# GEOS returns a segment's ends for fractions outside it, and the end of
# the line past its length.
point = np.where((fraction <= 0.0)[..., None], start, point)
point = np.where((fraction >= 1.0)[..., None], stop, point)
return np.where((distances >= reached[-1])[..., None], coordinates[-1], point)
def _dropped(points: np.ndarray, center, separation: float) -> np.ndarray:
"""Which of one layer's ``points`` the field leaves out, as a mask.
``points`` has shape ``(phases, points, 2)``. Of each pair at most
``separation`` apart the later point goes, and every point that close to
``center``, unless it is ``None``. A chain of close points can lose one
more than a pass in order would drop, but such a pass made a Contour
layout half again as slow. Layers lie a scale apart, farther than the
separation, so each is thinned on its own, and one tree holds every
phase, set side by side.
"""
from scipy.spatial import cKDTree # ty: ignore[unresolved-import]
n_phases, count = points.shape[:2]
limit = separation * separation
dropped = np.zeros((n_phases, count), dtype=bool)
if center is not None:
gap = points - center
dropped |= gap[..., 0] * gap[..., 0] + gap[..., 1] * gap[..., 1] <= limit
if count < 2:
return dropped
width = float(np.ptp(points[..., 0])) + 4.0 * separation
side_by_side = points.reshape(-1, 2).copy()
side_by_side[:, 0] += np.repeat(width * np.arange(n_phases), count)
# The shift rounds, so pairs are found a little beyond the separation and
# measured again where they lie.
pairs = cKDTree(side_by_side).query_pairs(
separation * (1.0 + 1e-9), output_type="ndarray"
)
if len(pairs):
flat = points.reshape(-1, 2)
gap = flat[pairs[:, 0]] - flat[pairs[:, 1]]
close = gap[:, 0] * gap[:, 0] + gap[:, 1] * gap[:, 1] <= limit
dropped.reshape(-1)[pairs[close].max(axis=1)] = True
return dropped
@lru_cache(maxsize=1)
def _deepest(geometry):
"""A depth that no inset of ``geometry`` reaches.
It lies past the deepest point by more than the tolerance that point is
found to. A buffer rounds a corner with chords, which carry an inset past
the deepest point by about a thousandth of it, within that tolerance.
"""
minx, miny, maxx, maxy = geometry.bounds
tolerance = _DEPTH_TOLERANCE * max(maxx - minx, maxy - miny)
radius = shapely.maximum_inscribed_circle(geometry, tolerance).length
return radius + 2.0 * tolerance
@lru_cache(maxsize=1)
def _outermost(geometry, clearance: float):
"""``geometry`` inset by ``clearance``, the outermost contour at every scale.
The candidate search tries a hundred scales or more on one boundary, and
each would buffer it again.
"""
return geometry.buffer(-clearance)
[docs]
class Contour(_CandidateLayout):
"""Place turbines along nested inward offsets of the farm boundary.
``spacing_ratio`` is the target distance along each contour divided by the
distance between contour layers, and ``include_center`` adds a turbine at
an interior point. Layer count and turbine allocation follow from the
requested count and the available geometry. The innermost contour hugs
the deepest point of the farm, so it can pass a few meters from the
center or fold back on itself, and a contour turning an acute corner
cuts across it. A candidate closer than 0.6 of the smaller of the two
spacings to an earlier one, the center first and then the contours from
the outside in, is therefore left out.
"""
refine_candidates = False
def __init__(
self,
n_turbines: int | None = None,
*,
density: dict[str, float] | None = None,
spacing_ratio: float = 1.0,
include_center: bool = True,
fitter: SelectedCandidateFitter | None = None,
seed: int | np.random.Generator | None = None,
):
self.spacing_ratio = _validate_positive("spacing_ratio", spacing_ratio)
self.include_center = bool(include_center)
if fitter is None:
fitter = SelectedContourFitter()
super().__init__(
n_turbines,
density=density,
fitter=fitter,
seed=seed,
)
[docs]
def estimate_spacing(self, boundary: Boundary) -> float:
"""Return a starting contour-layer spacing in meters.
Assumes one square ``area / n_turbines`` cell per turbine. The estimate
only seeds the fitter's search and is not the final spacing.
"""
return float(np.sqrt(boundary.area / self._required_n))
def _candidate_origin(self, boundary):
return boundary.interior_point
def _required_candidates(self, boundary):
if not self.include_center:
return np.empty((0, 2))
return np.asarray(boundary.interior_point, dtype=float).reshape(1, 2)
@staticmethod
def _edge_clearance(boundary):
"""Inset of the outermost contour, so its points lie strictly inside."""
minx, miny, maxx, maxy = boundary.bounds
return max(float(np.hypot(maxx - minx, maxy - miny)) * 1e-7, 1e-9)
def _candidate_bounds(self, boundary, scale, phases, max_candidates):
# Where the second layer lies past the deepest point, the outermost
# contour is the only one, so its rings' lengths count the field.
clearance = self._edge_clearance(boundary)
deepest = _deepest(boundary.geometry)
outer = _outermost(boundary.geometry, clearance)
lengths = np.array(
[]
if clearance >= deepest or outer.is_empty
else [line.length for line in _contour_lines(outer)]
)
minx, miny, maxx, maxy = boundary.bounds
hull = boundary.geometry.convex_hull
separation = _MIN_SEPARATION * min(1.0, self.spacing_ratio)
def bound(scales):
scales = np.asarray(scales, dtype=float)
counts = np.round(lengths / (scales[:, None] * self.spacing_ratio))
totals = int(self.include_center) + np.maximum(1, counts.astype(int)).sum(
axis=1
)
# Nearer the edge, more layers may hold any number of points.
totals[clearance + scales < deepest] = np.iinfo(totals.dtype).max
# Thinning keeps no two points of a layer, or the center, within
# the separation, and layers lie a scale apart, farther still. So
# the bounding box holds at most one point per square whose
# diagonal is the separation, and Oler's inequality bounds how many
# the convex hull holds. The separation is a little smaller for
# rounding.
apart = separation * scales * (1.0 - 1e-9)
side = apart / np.sqrt(2.0)
boxed = (np.floor((maxx - minx) / side) + 1.0) * (
np.floor((maxy - miny) / side) + 1.0
)
oler = np.floor(
2.0 / np.sqrt(3.0) * hull.area / apart**2
+ hull.length / (2.0 * apart)
+ 1.0
)
# Capped below 2**63, so a tiny scale stays a large integer.
packed = np.minimum(np.minimum(boxed, oler), 2.0**62).astype(totals.dtype)
return np.repeat(np.minimum(totals, packed)[:, None], len(phases), axis=1)
return bound
def _candidate_fields(
self, boundary, scale, phases, max_candidates=None, min_candidates=None
):
if scale <= 0:
return [np.empty((0, 2)) for _ in phases]
clearance = self._edge_clearance(boundary)
minx, miny, maxx, maxy = boundary.bounds
maximum_layers = int(np.ceil(np.hypot(maxx - minx, maxy - miny) / scale)) + 1
if max_candidates is not None and maximum_layers > max_candidates:
return [None for _ in phases]
# Insets, their rings and point counts are shared by every phase.
contours = []
total = int(self.include_center)
deepest = _deepest(boundary.geometry)
inset = boundary.geometry
for layer_index in range(maximum_layers):
depth = clearance + layer_index * scale
# A buffer past the deepest point can leave a stray piece, which
# may lie outside the boundary, so no layer is laid there.
if depth >= deepest:
break
if layer_index == 0:
inset = _outermost(boundary.geometry, clearance)
elif layer_index == 1:
inset = boundary.geometry.buffer(-depth)
else:
# Insets by disks add, so a layer insets the one before it,
# whose smoother outline buffers several times faster than a
# cell cleared around many turbines.
inset = inset.buffer(-scale)
if inset.is_empty:
break
for contour_index, line in enumerate(_contour_lines(inset)):
length = line.length
count = max(1, int(round(length / (scale * self.spacing_ratio))))
total += count
if max_candidates is not None and total > max_candidates:
return [None for _ in phases]
contours.append((layer_index, contour_index, line, length, count))
# Insets are not clipped and thinning only drops points, so no phase
# holds more than ``total``.
if min_candidates is not None and total < min_candidates:
return [np.empty((0, 2)) for _ in phases]
if not contours:
# The edge clearance reaches the deepest point, so no layer fits
# and only a center remains.
center = np.asarray([boundary.interior_point], dtype=float)
return [center if self.include_center else np.empty((0, 2)) for _ in phases]
counts = np.asarray([count for *_, count in contours], dtype=int)
layers, rings, lengths = (
np.asarray([contour[column] for contour in contours], dtype=float)
for column in (0, 1, 3)
)
# Every phase at once, with the per-phase operations in the same order.
contour_phase, layer_phase = np.asarray(phases, dtype=float).reshape(-1, 2).T
offsets = (
contour_phase[:, None] + layers * layer_phase[:, None] + 0.5 * rings
) % 1.0
# The center, if any, comes first in every phase.
first = int(self.include_center)
points = np.empty((len(phases), first + int(counts.sum()), 2))
center = None
if self.include_center:
center = np.asarray(boundary.interior_point, dtype=float)
points[:, 0] = center
separation = _MIN_SEPARATION * scale * min(1.0, self.spacing_ratio)
dropped = np.zeros(points.shape[:2], dtype=bool)
stops = first + np.cumsum(counts)
starts = stops - counts
# The contours of a layer follow one another, from the outermost.
layer_ends = np.append(np.flatnonzero(np.diff(layers)) + 1, len(contours))
layer_starts = np.append(0, layer_ends[:-1])
# Layer by layer, innermost first, as most turbines the thinning drops
# lie there, so a scale it leaves short is mostly found before the
# rest are placed.
for begin, end in zip(layer_starts[::-1], layer_ends[::-1]):
for ring in range(begin, end):
steps = np.arange(counts[ring], dtype=float)
points[:, starts[ring] : stops[ring]] = _along(
contours[ring][2],
(steps + offsets[:, ring, None]) * lengths[ring] / counts[ring],
)
low, high = starts[begin], stops[end - 1]
dropped[:, low:high] = _dropped(points[:, low:high], center, separation)
if min_candidates is not None and np.all(
(~dropped).sum(axis=1) < min_candidates
):
return [np.empty((0, 2)) for _ in phases]
return [field[~mask] for field, mask in zip(points, dropped)]
def _nominal_scale(self, boundary, spacing, phase) -> float:
# Layers lie the scale apart and their points about the scale times
# spacing_ratio. Corners, necks and the center put a few points much
# closer, which would mislead a spacing measured from the field.
return spacing / min(1.0, self.spacing_ratio)
def _detailed_fit(self, boundary, positions, fit):
from ..boundaries.base import _farm_pieces
scale = float(fit["scale"])
clearance = self._edge_clearance(boundary)
distances = shapely.distance(
boundary.geometry.boundary, shapely.points(positions)
)
layer_indices = np.rint(np.maximum(0.0, distances - clearance) / scale).astype(
int
)
center_mask = np.zeros(len(positions), dtype=bool)
if self.include_center:
center = np.asarray(boundary.interior_point, dtype=float)
center_mask = np.all(np.isclose(positions, center), axis=1)
pieces = _farm_pieces(boundary)
# A farm filled over several parts lists its turbines part by part, so
# its layers are ordered and recorded within each part.
groups = (
[np.arange(len(positions))]
if len(pieces) == 1
else _grouped_by_piece(pieces, positions)
)
orders = []
layers = []
start = 0
for component_index, group in enumerate(groups):
group_layers = layer_indices[group]
group_center = center_mask[group]
if self.include_center:
non_center = group_layers[~group_center]
center_layer = int(non_center.max() + 1) if len(non_center) else 0
group_layers[group_center] = center_layer
order = np.argsort(group_layers, kind="stable")
group_layers = group_layers[order]
group_center = group_center[order]
group_distances = distances[group][order]
first = 0
for layer_index in np.unique(group_layers):
count = int(np.count_nonzero(group_layers == layer_index))
is_center = bool(np.all(group_center[first : first + count]))
layer = {
"layer_index": int(layer_index),
"kind": "center" if is_center else "contour",
"inset_distance": (
float(group_distances[first])
if is_center
else float(clearance + layer_index * scale)
),
"n_turbines": count,
"start": start,
"stop": start + count,
}
if len(groups) > 1:
layer = {"component_index": component_index, **layer}
layers.append(layer)
first += count
start += count
orders.append(group[order])
positions = np.asarray(positions)[np.concatenate(orders)]
return positions, {**fit, "spacing_ratio": self.spacing_ratio, "layers": layers}
