"""Base class for boundary geometries.
A :class:`Boundary` wraps a Shapely ``Polygon`` or ``MultiPolygon`` in meters.
Exclusions, setbacks and partitions are
:class:`~oowl.boundaries.modifiers.BoundaryModifier` objects chained by a
:class:`~oowl.boundaries.composites.Composite`.
"""
from __future__ import annotations
from abc import ABC, abstractmethod
from copy import copy
from functools import cached_property
from typing import TYPE_CHECKING, Any
import numpy as np
import shapely
from shapely import affinity
from shapely.geometry import MultiPolygon, Polygon
from shapely.geometry.base import BaseGeometry
from .._random import resolve_seed
from .._repr import _ArgumentRepr, _describe
from .._validation import _validate_positive, _validate_real_array
if TYPE_CHECKING: # pragma: no cover
from collections.abc import Callable
import matplotlib.axes
from numpy.typing import ArrayLike
from .composites import Composite
from .modifiers import BoundaryModifier
__all__ = ["Boundary"]
Point2D = tuple[float, float]
def _resolve_center(center: Point2D | None) -> tuple[Point2D, bool]:
"""Return a usable center and whether the caller supplied it explicitly."""
if center is None:
return (0.0, 0.0), False
values = _validate_real_array("center", center)
if values.shape != (2,) or not np.all(np.isfinite(values)):
raise ValueError("center must contain two finite coordinates or be None")
return (float(values[0]), float(values[1])), True
def _centroid(geometry: BaseGeometry) -> Point2D:
"""Return the centroid of ``geometry`` as ``(x, y)``."""
centroid = geometry.centroid
return (centroid.x, centroid.y)
def _oriented(
geometry: BaseGeometry, rotation: float | None, rng: np.random.Generator
) -> BaseGeometry:
"""Rotate ``geometry`` about its centroid by ``rotation`` degrees.
A ``rotation`` of ``None`` draws a uniform heading from ``rng``.
"""
angle = rng.uniform(0.0, 360.0) if rotation is None else rotation
return affinity.rotate(geometry, angle, origin="centroid", use_radians=False)
def _translated_copies(geometry: BaseGeometry, offsets) -> np.ndarray:
"""Return ``geometry`` translated by each ``(dx, dy)`` row of ``offsets``."""
offsets = np.asarray(offsets, dtype=float).reshape(-1, 2)
shifts = np.repeat(offsets, int(shapely.get_num_coordinates(geometry)), axis=0)
def shift(coordinates):
# The full affine product of ``affinity.translate``, so each copy is
# bit-identical to translating ``geometry`` one offset at a time.
x, y = coordinates.T
return np.stack(
[
1.0 * x + 0.0 * y + shifts[:, 0],
0.0 * x + 1.0 * y + shifts[:, 1],
]
).T
return shapely.transform(np.full(len(offsets), geometry, dtype=object), shift)
def _project_inside(
geometry: BaseGeometry, interior: BaseGeometry, points
) -> np.ndarray:
"""Move points outside ``geometry`` to their nearest point in ``interior``."""
projected = np.array(points, dtype=float).reshape(-1, 2)
# A point is covered exactly when it intersects, and testing coordinates
# builds no point geometries.
outside = ~shapely.intersects_xy(geometry, projected[:, 0], projected[:, 1])
if np.any(outside):
lines = shapely.shortest_line(interior, shapely.points(projected[outside]))
projected[outside] = shapely.get_coordinates(lines).reshape(-1, 2, 2)[:, 0]
return projected
def _edge_clearance(geometry: BaseGeometry) -> float:
"""Inset in meters that keeps a point projected inside after rounding.
A tiny fraction of the extent, but at least 128 coordinate units in the
last place. Shapely's inward buffer misplaces the edge below about 112 of
them, which matters for a small farm far from the origin.
"""
minx, miny, maxx, maxy = geometry.bounds
magnitude = max(abs(minx), abs(miny), abs(maxx), abs(maxy))
return max(
float(np.hypot(maxx - minx, maxy - miny)) * 1e-10,
1e-12,
128.0 * float(np.spacing(magnitude)),
)
def _inset_interior(geometry: BaseGeometry) -> BaseGeometry:
"""Return ``geometry`` shrunk by a rounding-sized clearance to project into.
A point projected onto the edge itself can test outside after rounding.
"""
interior = geometry.buffer(-_edge_clearance(geometry))
return geometry if interior.is_empty else interior
def _validate_target_area(area: float) -> float:
"""Return ``area`` as a finite positive float in square meters."""
return _validate_positive("area", area)
def _optional_area(area: float | None) -> float | None:
"""``area`` validated as for :func:`_validate_target_area`, or ``None``."""
return None if area is None else _validate_target_area(area)
def _check_geometry(geometry: BaseGeometry) -> None:
"""Reject an empty, zero-area or invalid boundary geometry."""
if geometry.is_empty:
raise ValueError("boundary geometry must have positive area")
# The shoelace area of a symmetric bow-tie cancels to zero, so an invalid
# shape of zero area names its defect too.
flat = geometry.area <= 0
if not geometry.is_valid:
need = "have positive area and be valid" if flat else "be valid"
raise ValueError(
f"boundary geometry must {need}, not {shapely.is_valid_reason(geometry)}"
)
if flat:
raise ValueError("boundary geometry must have positive area")
[docs]
class Boundary(_ArgumentRepr, ABC):
"""Abstract wind-farm boundary backed by a Shapely geometry in meters.
Subclasses build an unscaled shape about the origin, which is scaled about
its centroid to ``target_area`` and then moved so its centroid lies at
``center``, where either is set. Placing the shape last keeps a site far
from the origin, as in UTM coordinates, as accurate as one at it.
"""
# ``center`` shows the value passed, since an omitted center resolves to
# the origin but leaves the boundary free to be placed.
_repr_attributes = {"area": "target_area", "center": "_passed_center"}
#: Whether a caller supplied an explicit seed.
_pinned_seed = False
#: Whether resizing must re-apply something given in meters, see
#: :class:`Composite`. A plain shape scales as a whole.
_in_meters = False
#: Whether the shape only describes a form, built at unit size until an
#: ``area`` sizes it. A geometry given in meters is already sized.
_unit_shape = True
def _resolve_own_seed(self, seed) -> int:
"""Resolve ``seed``, remembering whether the caller pinned it."""
self._pinned_seed = seed is not None
return resolve_seed(seed)
@abstractmethod
def _build_geometry(self) -> BaseGeometry:
"""Construct and return the Shapely geometry for this boundary."""
def _scale_to_target_area(self, geometry: BaseGeometry) -> BaseGeometry:
target_area = getattr(self, "target_area", None)
if target_area is None:
if not self._unit_shape:
return geometry
# Not yet sized, so built at unit size to be sized later.
target_area = 1.0
factor = float(np.sqrt(target_area / geometry.area))
return affinity.scale(
geometry, xfact=factor, yfact=factor, origin=_centroid(geometry)
)
@property
def _sized(self) -> bool:
"""Whether something gave this boundary a size in meters.
A shape built without ``area`` is unit size, to be scaled by its user
or sized by a sampler, and a farm is not laid out on it.
"""
return getattr(self, "target_area", None) is not None or not self._unit_shape
def _translate_to_center(self, geometry: BaseGeometry) -> BaseGeometry:
target = getattr(self, "center", None)
if target is None:
return geometry
cx, cy = _centroid(geometry)
return affinity.translate(
geometry,
xoff=float(target[0]) - cx,
yoff=float(target[1]) - cy,
)
@cached_property
def _unscaled(self) -> Polygon | MultiPolygon:
"""The checked shape before scaling and moving, kept by resized copies."""
geometry = self._build_geometry()
if not isinstance(geometry, (Polygon, MultiPolygon)):
raise TypeError("boundary geometry must be a Polygon or MultiPolygon")
_check_geometry(geometry)
return geometry
@cached_property
def geometry(self) -> Polygon | MultiPolygon:
"""The Shapely ``Polygon`` or ``MultiPolygon`` of this boundary, in meters.
Raises ``TypeError`` for another geometry type and ``ValueError`` for an
empty, zero-area or invalid one.
"""
geometry = self._scale_to_target_area(self._unscaled)
geometry = self._translate_to_center(geometry)
# Scaling and moving can round a nearly degenerate geometry invalid.
_check_geometry(geometry)
return geometry
@property
def area(self) -> float:
"""Enclosed area in square meters."""
return self.geometry.area
@cached_property
def bounds(self) -> tuple[float, float, float, float]:
"""Axis-aligned bounding box as ``(min_x, min_y, max_x, max_y)`` in meters."""
return self.geometry.bounds
@cached_property
def centroid(self) -> Point2D:
"""Geometric center of mass as ``(x, y)`` in meters.
It may lie outside a concave or multipart boundary. See
:attr:`interior_point`.
"""
return _centroid(self.geometry)
@cached_property
def interior_point(self) -> Point2D:
"""A point inside the boundary, as ``(x, y)`` in meters.
The centroid if contained, otherwise Shapely's representative point.
"""
centroid = self.centroid
if bool(self._contains(np.asarray(centroid, dtype=float).reshape(1, 2))[0]):
return centroid
point = self.geometry.representative_point()
return (point.x, point.y)
@cached_property
def _enclosing_halfplanes(self) -> tuple[np.ndarray, np.ndarray]:
"""Outward unit normals ``(4, 2)`` and limits ``(4,)`` of a box around it.
The box is the minimum rotated rectangle, widened by a margin that
covers its rounding, so a point with ``normals @ p > limits`` for any
row lies outside the boundary and no point of its edge is cut off.
"""
corners = np.asarray(
shapely.minimum_rotated_rectangle(self.geometry).exterior.coords,
dtype=float,
)[:4]
edges = np.roll(corners, -1, axis=0) - corners
normals = np.column_stack([edges[:, 1], -edges[:, 0]])
normals /= np.linalg.norm(normals, axis=1, keepdims=True)
inward = np.einsum("ij,ij->i", normals, corners - corners.mean(axis=0)) < 0
normals[inward] *= -1.0
minx, miny, maxx, maxy = self.bounds
margin = 1e-9 * float(np.hypot(maxx - minx, maxy - miny)) + 1e-6
return normals, np.einsum("ij,ij->i", normals, corners) + margin
@cached_property
def _edges(self) -> tuple[np.ndarray, np.ndarray]:
"""Start and end points ``(E, 2)`` of every edge of every ring."""
rings = shapely.get_rings(shapely.get_parts(self.geometry))
coordinates = [shapely.get_coordinates(ring) for ring in rings]
return (
np.concatenate([ring[:-1] for ring in coordinates]),
np.concatenate([ring[1:] for ring in coordinates]),
)
@property
def parts(self) -> tuple["Boundary", ...]:
"""The boundaries a farm is built from, this one alone without a partition."""
return (self,)
[docs]
def contains(self, points: ArrayLike) -> np.ndarray:
"""Test which points lie in the boundary, counting the edge as inside.
Parameters
----------
points:
Array-like of shape ``(N, 2)``, or one point of shape ``(2,)``.
Returns
-------
numpy.ndarray
Boolean array of length ``N``.
Raises
------
TypeError
For points that are not numbers.
ValueError
For another shape or a point that is not finite.
"""
from .._validation import _validate_real_array
pts = _validate_real_array("points", points)
if pts.size == 0:
pts = pts.reshape(0, 2)
elif pts.shape == (2,):
pts = pts.reshape(1, 2)
if pts.ndim != 2 or pts.shape[1] != 2:
raise ValueError(f"points must have shape (N, 2) or (2,), not {pts.shape}")
if not np.all(np.isfinite(pts)):
raise ValueError("points must contain only finite values")
return self._contains(pts)
def _contains(self, points: np.ndarray) -> np.ndarray:
""":meth:`contains` for finite float points of shape ``(N, 2)``."""
return shapely.intersects_xy(self.geometry, points[:, 0], points[:, 1])
@property
def _is_random(self) -> bool:
"""Whether anything in this boundary is drawn rather than given."""
return hasattr(self, "seed")
def _uncached_copy(self) -> "Boundary":
"""Return a shallow copy whose cached geometry is rebuilt on first use."""
boundary = copy(self)
for name in (
"_unscaled",
"geometry",
"bounds",
"centroid",
"interior_point",
"_enclosing_halfplanes",
"_edges",
):
boundary.__dict__.pop(name, None)
return boundary
def _reseeded(self, seed: int) -> "Boundary":
"""Return a copy whose random geometry is rebuilt from ``seed``."""
boundary = self._uncached_copy()
if hasattr(boundary, "seed"):
setattr(boundary, "seed", int(seed))
return boundary
def _at_base_area(self, base_area: float) -> "Boundary":
"""This boundary sized for a spacing search, see ``Composite._at_base_area``.
A boundary whose shape scales as a whole has no base to size apart, so
this is a resize to ``base_area``.
"""
return self._resized(base_area)
def _resized(self, target_area: float) -> "Boundary":
"""Return a copy of the same shape, same seed, scaled to ``target_area``.
No builder reads the target area, so the shape is built once here and
shared by every resized copy. A composite whose modifiers are given in
meters resizes its base instead, see :class:`Composite`.
"""
target_area = _validate_target_area(target_area)
unscaled = self._unscaled
boundary = self._uncached_copy()
setattr(boundary, "target_area", target_area)
boundary.__dict__["_unscaled"] = unscaled
return boundary
def _resample(self, seed: int) -> "Boundary":
"""Return an unpositioned copy with independently resolved randomness."""
boundary = self._reseeded(seed)
setattr(boundary, "_positioned", False)
return boundary
[docs]
def modify(self, *modifiers: "BoundaryModifier") -> "Composite":
"""Return a :class:`Composite` applying ``modifiers`` in order.
This boundary is left unchanged.
"""
from .composites import Composite
return Composite(self, modifiers)
[docs]
def plot(
self,
ax: "matplotlib.axes.Axes | None" = None,
**kwargs: Any,
) -> "matplotlib.axes.Axes":
"""Plot the boundary and return its Matplotlib axes."""
from ..plotting import _plot_boundary
return _plot_boundary(self.geometry, ax=ax, **kwargs)
class _GeometryBoundary(Boundary):
"""Internal boundary view over an existing polygon geometry.
``cells`` are the partition cells of the boundary it stands in for, so a
farm placed in them keeps its parts in the same order.
"""
def __init__(
self,
geometry: Polygon | MultiPolygon,
*,
positioned: bool = True,
cells: tuple[BaseGeometry, ...] | None = None,
):
self._geometry = geometry
self._positioned = bool(positioned)
self._cells = cells
def _build_geometry(self) -> BaseGeometry:
return self._geometry
@property
def parts(self) -> tuple[Boundary, ...]:
"""The cells as boundaries, or this view alone without them."""
if self._cells is None:
return (self,)
cells = self._cells
target_area = getattr(self, "target_area", None)
if target_area is not None:
# The cells take the similarity transform the shape does.
factor = float(np.sqrt(target_area / self._geometry.area))
origin = _centroid(self._geometry)
cells = tuple(
affinity.scale(cell, xfact=factor, yfact=factor, origin=origin)
for cell in cells
)
return tuple(_GeometryBoundary(cell) for cell in cells)
def __repr__(self) -> str:
# Printed as the public call that wraps a geometry the same way.
area = getattr(self, "target_area", None)
return (
f"Polygon.from_geometry({_describe(self._geometry)}, "
f"area={_describe(area)}, center=None)"
)
def metadata(self) -> dict[str, Any]:
"""JSON-compatible record naming the public class of the same boundary."""
return {"type": "Polygon"}
# A geometry wrapped as it is has its size in meters already.
_unit_shape = False
def _cells(
boundary: Boundary,
transform: Callable[[BaseGeometry], BaseGeometry] | None = None,
) -> tuple[BaseGeometry, ...] | None:
"""Cell geometries of a boundary with several parts, through ``transform``.
``None`` for a boundary of one part, whose farm fills its polygons.
"""
parts = boundary.parts
if len(parts) < 2:
return None
return tuple(
part.geometry if transform is None else transform(part.geometry)
for part in parts
)
def _translation(offset) -> Callable[[BaseGeometry], BaseGeometry]:
"""A transform moving a geometry by ``offset``, an ``(x, y)`` pair in meters."""
xoff, yoff = (float(value) for value in offset)
return lambda geometry: affinity.translate(geometry, xoff=xoff, yoff=yoff)
def _farm_pieces(boundary: Boundary) -> tuple[BaseGeometry, ...]:
"""The geometries a farm on ``boundary`` fills one by one, in order.
The partition cells in ``parts`` order when there are several, so the
gap between cells does not change how a farm is placed, else the
polygons of the shape.
"""
cells = _cells(boundary)
if cells is not None:
return cells
geometry = boundary.geometry
return tuple(getattr(geometry, "geoms", (geometry,)))
