"""Boundary-grid layout parameterization from Stanley and Ning (2019).
Reference: https://doi.org/10.5194/wes-4-663-2019
"""
from __future__ import annotations
from functools import lru_cache
from typing import TYPE_CHECKING
import numpy as np
import shapely
from .._validation import _validate_unit_interval
from .base import (
_MAX_FILL_TURBINES,
LayoutAlgorithm,
_apart_from,
_closest_distance,
_NoFitError,
_spread_selection,
_too_many_error,
)
from .perimeter import _contained_perimeter_points, _spaced_perimeter_points
if TYPE_CHECKING: # pragma: no cover
from ..boundaries.base import Boundary
__all__ = ["BoundaryGrid"]
#: Row spacing over column spacing, ``dy = 4 dx`` from Stanley and Ning (2019).
_GRID_SPACING_RATIO = 4.0
#: Row offset angle in degrees from Stanley and Ning (2019).
_ROW_OFFSET_ANGLE = 20.0
#: Perimeter starts and grid rotations tried, keeping the best spaced. Eight
#: rotations stratified over a half turn spaced 80 farms as well as sixteen
#: drawn over the full circle, in a third of the time, and filled 2 % fewer
#: turbines at a spacing.
_N_REALIZATIONS = 8
#: Relative widening of the filled grid's column spacing, so rounding in the
#: lattice coordinates never brings a pair below the spacing.
_LATTICE_MARGIN = 1e-9
#: Row offset per unit of column spacing.
_ROW_OFFSET = float(np.tan(np.deg2rad(_ROW_OFFSET_ANGLE)))
def _realizations(rng):
"""Random perimeter starts beside grid rotations stratified over 180 degrees.
The lattice is symmetric under a half turn about its center, so rotations
drawn over the full circle tried half their lattices twice. One rotation
is drawn in each of :data:`_N_REALIZATIONS` equal sectors of a half turn.
"""
width = 180.0 / _N_REALIZATIONS
for index in range(_N_REALIZATIONS):
yield float(rng.random()), float((index + rng.random()) * width)
@lru_cache(maxsize=64)
def _rotation_transpose(rotation: float) -> np.ndarray:
"""The transposed rotation matrix of ``rotation`` degrees, read-only."""
angle = np.deg2rad(rotation)
matrix = np.array(
[
[np.cos(angle), -np.sin(angle)],
[np.sin(angle), np.cos(angle)],
]
).T
matrix.flags.writeable = False
return matrix
def _grid_basis(x_spacing, rotation):
y_spacing = _GRID_SPACING_RATIO * x_spacing
row_offset = _ROW_OFFSET * x_spacing
basis = np.array([[x_spacing, 0.0], [row_offset, y_spacing]])
return basis @ _rotation_transpose(rotation)
def _grid_points(geometry, center, x_spacing, rotation):
if x_spacing <= 0 or geometry.is_empty:
return np.empty((0, 2))
indices = _lattice_indices(geometry, center, x_spacing, rotation)
points = indices @ _grid_basis(x_spacing, rotation) + center
minx, miny, maxx, maxy = geometry.bounds
# The rotated index range covers far more than the bounding box, and no
# point outside the box can touch the boundary.
in_box = (
(points[:, 0] >= minx)
& (points[:, 0] <= maxx)
& (points[:, 1] >= miny)
& (points[:, 1] <= maxy)
)
points = points[in_box]
mask = shapely.intersects_xy(geometry, points[:, 0], points[:, 1])
return points[mask]
def _lattice_indices(geometry, center, x_spacing, rotation):
"""Lattice indices whose points at ``x_spacing`` or more may lie in the box.
Larger spacings only narrow the range, so the indices cover them too.
"""
basis = _grid_basis(x_spacing, rotation)
minx, miny, maxx, maxy = geometry.bounds
corners = np.array([[minx, miny], [minx, maxy], [maxx, miny], [maxx, maxy]])
coefficients = (corners - center) @ np.linalg.inv(basis)
lower = np.floor(coefficients.min(axis=0)).astype(int) - 1
upper = np.ceil(coefficients.max(axis=0)).astype(int) + 1
first = np.arange(lower[0], upper[0] + 1)
second = np.arange(lower[1], upper[1] + 1)
return np.column_stack([np.tile(first, len(second)), np.repeat(second, len(first))])
class _GridCounter:
"""The number of grid points in ``geometry`` at spacings inside a bracket.
Each point moves along a line from the center as the spacing grows, so
over a bracket it sweeps a short segment. A segment farther from the edge
than rounding could carry it stays on one side, and only the points whose
segments come near the edge are tested at each spacing. The count equals
that of :func:`_grid_points`, so a search over it decides alike.
"""
def __init__(self, geometry, center, rotation):
self.geometry = geometry
self.center = center
self.rotation = rotation
self.edge = geometry.boundary
shapely.prepare(self.edge)
minx, miny, maxx, maxy = geometry.bounds
self.box = (minx, miny, maxx, maxy)
size = max(abs(minx), abs(miny), abs(maxx), abs(maxy), *np.abs(center))
# Far beyond the rounding of a coordinate, far below a spacing.
self.margin = 1e-9 * max(size, 1.0)
def _points(self, indices, x_spacing):
return indices @ _grid_basis(x_spacing, self.rotation) + self.center
def classify(self, lower, upper):
"""Sort the points over ``lower`` to ``upper`` into sure and near ones."""
indices = _lattice_indices(self.geometry, self.center, lower, self.rotation)
start = self._points(indices, lower)
stop = self._points(indices, upper)
minx, miny, maxx, maxy = self.box
# A segment wholly past one side of the box lies outside.
beyond = (
(np.maximum(start[:, 0], stop[:, 0]) < minx - self.margin)
| (np.minimum(start[:, 0], stop[:, 0]) > maxx + self.margin)
| (np.maximum(start[:, 1], stop[:, 1]) < miny - self.margin)
| (np.minimum(start[:, 1], stop[:, 1]) > maxy + self.margin)
)
keep = np.flatnonzero(~beyond)
segments = shapely.linestrings(np.stack((start[keep], stop[keep]), axis=1))
clear = ~shapely.dwithin(self.edge, segments, self.margin)
sure = keep[clear]
inside = shapely.intersects_xy(self.geometry, start[sure, 0], start[sure, 1])
self.sure_inside = int(np.count_nonzero(inside))
self.near = indices[keep[~clear]]
def count(self, x_spacing):
"""Grid points in the geometry at ``x_spacing``, inside the bracket."""
points = self._points(self.near, x_spacing)
near = shapely.intersects_xy(self.geometry, points[:, 0], points[:, 1])
return self.sure_inside + int(np.count_nonzero(near))
def _fit_grid(geometry, count, rotation):
if count == 0:
return np.empty((0, 2)), float("inf"), 0
center = _grid_center(geometry)
if count == 1:
# The lattice is centered on a point inside, which it keeps at every
# spacing, so the search for the widest grid holding one point grew
# its spacing 80 times and bisected past 1e14 m. Beyond the bounding
# box's diagonal the center is the only point left.
min_x, min_y, max_x, max_y = geometry.bounds
return center[None, :], float(np.hypot(max_x - min_x, max_y - min_y)), 1
initial = max(np.sqrt(geometry.area / (count * _GRID_SPACING_RATIO)), 1e-12)
lower = initial
for _ in range(200):
if len(_grid_points(geometry, center, lower, rotation)) >= count:
break
lower /= 1.5
else:
return np.empty((0, 2)), 0.0, 0
upper = lower * 1.5
for _ in range(80):
if len(_grid_points(geometry, center, upper, rotation)) < count:
break
lower = upper
upper *= 1.5
# One classification over the whole bracket leaves a few dozen points near
# the edge, so each step tests those alone.
counter = _GridCounter(geometry, center, rotation)
counter.classify(lower, upper)
for _ in range(50):
if upper - lower <= 1e-6 * max(lower, 1.0):
break
midpoint = 0.5 * (lower + upper)
if counter.count(midpoint) >= count:
lower = midpoint
else:
upper = midpoint
candidates = _grid_points(geometry, center, lower, rotation)
return candidates, float(lower), int(len(candidates))
def _grid_center(geometry):
"""The centroid of ``geometry``, or an interior point when it falls outside."""
centroid = geometry.centroid
if not geometry.covers(centroid):
centroid = geometry.representative_point()
return np.array(centroid.coords[0])
def _inner_geometry(geometry, boundary_count: int, grid_count: int):
"""``geometry`` inset so an inner grid of ``grid_count`` keeps off the ring.
The inset is the spacing of ``boundary_count`` turbines around the
exterior or the row spacing of the grid over the whole area, whichever is
less, so a pair across ring and grid lies about as far apart as the closer
of the two.
"""
ring_spacing = geometry.exterior.length / boundary_count
row_spacing = np.sqrt(geometry.area / (grid_count * _GRID_SPACING_RATIO))
return geometry.buffer(-min(ring_spacing, row_spacing))
def _grid_option(geometry, boundary_points, grid_count: int, rotation, rng):
"""The ring and a grid of ``grid_count`` fitted in ``geometry``, or ``None``.
Returns the positions, the grid's spacing and candidate count, and the
closest pair, or ``None`` when the grid does not fit.
"""
candidates, x_spacing, candidate_count = _fit_grid(geometry, grid_count, rotation)
if len(candidates) < grid_count:
return None
grid_points = _select_interior(candidates, grid_count, boundary_points, rng)
positions = np.vstack([boundary_points, grid_points])
return {
"positions": positions,
"x_spacing": x_spacing,
"candidate_count": candidate_count,
"minimum_distance": _closest_distance(positions),
}
def _crowds_ring(option, geometry, boundary_count: int, grid_count: int) -> bool:
"""Whether a pair of ``option`` lies under half its ring or row spacing.
The inner grid is fitted to the whole boundary, so a grid turbine could
sit on the edge, 105 m from a ring turbine where both spaced 250 m and
more. Such a layout is compared with a grid inset from the ring, which a
few turbines in a small cell would crowd instead.
"""
if not boundary_count or not grid_count:
return False
ring_spacing = geometry.exterior.length / boundary_count
return option["minimum_distance"] < 0.5 * min(ring_spacing, option["x_spacing"])
def _ring_cap(fraction, grid_count):
"""Most ring turbines that hold at most ``fraction`` of the farm, ``None`` if any."""
if fraction >= 1.0:
return None
# The exact-count rule, ring count at most floor(fraction * total), solved
# for the ring, then stepped down where rounding overshot it.
count = int(np.floor(fraction * grid_count / (1.0 - fraction))) + 1
while count > int(np.floor(fraction * (count + grid_count))):
count -= 1
return count
def _ringed_grid(geometry, candidates, start, spacing, fraction):
"""The spaced ring around ``geometry`` and the ``candidates`` clear of it.
The ring holds at most ``fraction`` of the turbines. Grid turbines keep
the spacing from each ring turbine rather than from the whole edge, so
those between ring turbines stay.
"""
ring = _spaced_perimeter_points(
geometry, start, spacing, _ring_cap(fraction, len(candidates))
)
grid = _apart_from(candidates, ring, spacing)
# Grid turbines dropped near the ring can leave it above its share.
while (most := _ring_cap(fraction, len(grid))) is not None and len(ring) > most:
ring = _spaced_perimeter_points(geometry, start, spacing, most)
grid = _apart_from(candidates, ring, spacing)
return ring, grid
def _select_interior(candidates, count, boundary_points, rng):
if count == 0:
return np.empty((0, 2))
if len(candidates) == count:
return candidates
return _spread_selection(candidates, count, boundary_points, rng)
[docs]
class BoundaryGrid(LayoutAlgorithm):
"""Automatic boundary-grid layout based on Stanley and Ning (2019).
Up to ``boundary_fraction`` of the turbines are equally spaced around the
exterior boundary. The remainder occupy a centered, rotated inner grid.
Of the perimeter starts and grid rotations tried, the one with the widest
closest pair is kept, and when a grid turbine crowds the ring it is
refitted inside the boundary inset by the ring's spacing or the grid's
row spacing, whichever is less, keeping the wider of the two. The paper's
initialization values ``dy = 4 dx`` and a 20 degree row offset are fixed.
The seed chooses perimeter start and grid rotation.
With ``minimum_spacing`` the ring and the grid are both filled at that
spacing, the ring holding at most ``boundary_fraction`` of the turbines.
"""
def __init__(
self,
n_turbines: int | None = None,
*,
density: dict[str, float] | None = None,
boundary_fraction: float = 0.45,
seed: int | np.random.Generator | None = None,
):
self.boundary_fraction = _validate_unit_interval(
"boundary_fraction", boundary_fraction
)
super().__init__(
n_turbines,
density=density,
seed=seed,
)
def _generate_single(self, boundary, rng, seed):
geometry = boundary.geometry
exterior = geometry.exterior
n_turbines = self._required_n
maximum_boundary_count = int(np.floor(self.boundary_fraction * n_turbines))
best = None
for start, rotation in _realizations(rng):
for boundary_count in range(maximum_boundary_count, -1, -1):
boundary_points = _contained_perimeter_points(
geometry, boundary_count, start
)
grid_count = n_turbines - boundary_count
option = _grid_option(
geometry, boundary_points, grid_count, rotation, rng
)
if option is None:
continue
option.update(
boundary_count=boundary_count,
grid_count=grid_count,
start=start,
rotation=rotation,
)
if (
best is None
or option["minimum_distance"] > best["minimum_distance"]
):
best = option
break
if best is None:
raise ValueError("could not construct a feasible boundary-grid layout")
if _crowds_ring(best, geometry, best["boundary_count"], best["grid_count"]):
# Refit the chosen realization with the grid inset from its ring.
inner = _inner_geometry(
geometry, best["boundary_count"], best["grid_count"]
)
cleared = (
None
if inner.is_empty
else _grid_option(
inner,
_contained_perimeter_points(
geometry, best["boundary_count"], best["start"]
),
best["grid_count"],
best["rotation"],
rng,
)
)
if (
cleared is not None
and cleared["minimum_distance"] > best["minimum_distance"]
):
best = {**best, **cleared}
return self._single_result(
best["positions"],
boundary,
seed,
self._fit_record(best, best["start"] * exterior.length),
)
def _fit_record(self, best: dict, boundary_start: float | None) -> dict:
"""The fit record of the ``best`` realization."""
# Without inner-grid turbines there is no grid, so no grid spacing.
x_spacing = best["x_spacing"] if best["grid_count"] else None
realized_minimum = best["minimum_distance"]
if not np.isfinite(realized_minimum):
realized_minimum = None
return {
"method": "boundary_grid",
"boundary_fraction": self.boundary_fraction,
"boundary_count": best["boundary_count"],
"grid_count": best["grid_count"],
"boundary_start": boundary_start,
"boundary_start_fraction": best["start"],
"x_spacing": x_spacing,
"y_spacing": None if x_spacing is None else _GRID_SPACING_RATIO * x_spacing,
"row_offset": (
None
if x_spacing is None
else float(np.tan(np.deg2rad(_ROW_OFFSET_ANGLE)) * x_spacing)
),
"row_offset_angle": _ROW_OFFSET_ANGLE,
"rotation": best["rotation"],
"candidate_count": best["candidate_count"],
"realized_minimum_spacing": realized_minimum,
}
_fills_directly = True
def _fill_positions(
self, boundary: Boundary, spacing: float, rng: np.random.Generator
) -> tuple[np.ndarray, dict]:
"""Fill each part with a spaced ring and a spaced grid inside it.
The ring holds the most turbines no closer than the spacing, but at
most ``boundary_fraction`` of its part, spread evenly when it holds
fewer. The grid's column spacing is the spacing, the lattice's
shortest vector since its rows lie four times as far apart, and it
keeps only turbines that far from the ring. A part whose grid holds
more without a ring goes without one. Of the perimeter starts and
rotations the exact count tries, the one holding the most turbines
wins, ties going to the larger closest pair. ``boundary_start`` is
recorded for a single part only, since each part starts its ring at
the same fraction of its own exterior.
"""
from ..boundaries.base import _farm_pieces
fraction = self.boundary_fraction
x_spacing = spacing * (1.0 + _LATTICE_MARGIN)
# A cell in several parts rings each, as the exterior of each is its edge.
pieces = [
polygon
for piece in _farm_pieces(boundary)
for polygon in getattr(piece, "geoms", (piece,))
]
# The grid is drawn over each part's bounding box before it is clipped.
cell = _GRID_SPACING_RATIO * x_spacing**2
most = 0
for piece in pieces:
minx, miny, maxx, maxy = piece.bounds
most += int(piece.length // spacing + (maxx - minx) * (maxy - miny) / cell)
if most > _MAX_FILL_TURBINES:
raise _too_many_error(spacing, most, "ring and grid")
centers = [_grid_center(piece) for piece in pieces]
# Any realization beats this key, so ``best`` is always replaced.
best_key: tuple[int, float] = (-1, 0.0)
best: dict = {}
positions = np.empty((0, 2))
for start, rotation in _realizations(rng):
placed = np.empty((0, 2))
boundary_count = grid_count = candidate_count = 0
for index, piece in enumerate(pieces):
candidates = _grid_points(piece, centers[index], x_spacing, rotation)
ring, inner = np.empty((0, 2)), candidates
if fraction > 0:
ring_option = _ringed_grid(
piece, candidates, start, spacing, fraction
)
# The ring is the method, so it wins ties.
if sum(map(len, ring_option)) >= len(candidates):
ring, inner = ring_option
# Parts may lie closer than the spacing, so later ones give way.
ring = _apart_from(ring, placed, spacing)
inner = _apart_from(inner, placed, spacing)
placed = np.vstack([placed, ring, inner])
boundary_count += len(ring)
grid_count += len(inner)
candidate_count += len(candidates)
closest = _closest_distance(placed)
if (len(placed), closest) > best_key:
best_key, positions = (len(placed), closest), placed
best = {
"boundary_count": boundary_count,
"grid_count": grid_count,
"start": start,
"rotation": rotation,
"x_spacing": x_spacing,
"candidate_count": candidate_count,
"minimum_distance": closest,
}
if len(positions) < 2:
# The base error would claim no layout fits, where a ring alone may.
raise _NoFitError(
f"no boundary grid of two or more turbines with at most "
f"{fraction:g} of them on the boundary fits "
f"{float(boundary.area) / 1e6:.3g} km2 with {spacing:.3g} m spacing"
)
boundary_start = best["start"] * pieces[0].exterior.length
return positions, self._fit_record(
best, boundary_start if len(pieces) == 1 else None
)
