"""Offshore wind-farm boundaries sampled from a fitted model.
A support-function PCA generates a convex envelope, then a second PCA generates
inward edge paths on a beta-binomially sampled subset of its sides. The packaged
model was fitted to quality-controlled, operating EMODnet wind-farm polygons.
"""
from __future__ import annotations
import json
from dataclasses import dataclass
from functools import cache
from pathlib import Path
import numpy as np
from shapely.geometry import LineString, Polygon
from shapely.geometry.base import BaseGeometry
from shapely.geometry.polygon import orient
from .._random import as_rng
from .._validation import _validate_finite, _validate_integer
from .base import (
Boundary,
Point2D,
_optional_area,
_oriented,
_resolve_center,
)
from .fourier import _simplify_to_v
__all__ = ["OffshorePolygon"]
_DEFAULT_MODEL_PATH = Path(__file__).parents[1] / "data" / "offshore_polygon_model.npz"
@dataclass
class _OffshoreModel:
"""Fitted support-envelope and edge-path model of offshore lease outlines."""
support_mean: np.ndarray
support_components: np.ndarray
support_score_std: np.ndarray
support_score_location: float
support_score_scale: float
hull_vertex_log_mean: float
hull_vertex_log_std: float
hull_vertex_min: int
hull_vertex_max: int
path_mean: np.ndarray
path_components: np.ndarray
path_feature_mean: np.ndarray
path_feature_cholesky: np.ndarray
path_landmarks: int
path_count_min: int
path_count_max: int
path_posterior_degrees: int
path_posterior_scale: float
occupancy_alpha: float
occupancy_beta: float
metadata: dict
@classmethod
def load(cls, path: str | Path) -> "_OffshoreModel":
"""Load a packaged model from the ``.npz`` archive at ``path``.
Reading is pickle-free and gated on ``schema_version``. Any version
other than 1 raises ``ValueError``.
"""
with np.load(Path(path).with_suffix(".npz"), allow_pickle=False) as data:
schema_version = int(data["schema_version"])
if schema_version != 1:
raise ValueError(
f"unsupported offshore polygon model schema {schema_version}"
)
return cls(
support_mean=data["support_mean"],
support_components=data["support_components"],
support_score_std=data["support_score_std"],
support_score_location=float(data["support_score_location"]),
support_score_scale=float(data["support_score_scale"]),
hull_vertex_log_mean=float(data["hull_vertex_log_mean"]),
hull_vertex_log_std=float(data["hull_vertex_log_std"]),
hull_vertex_min=int(data["hull_vertex_min"]),
hull_vertex_max=int(data["hull_vertex_max"]),
path_mean=data["path_mean"],
path_components=data["path_components"],
path_feature_mean=data["path_feature_mean"],
path_feature_cholesky=data["path_feature_cholesky"],
path_landmarks=int(data["path_landmarks"]),
path_count_min=int(data["path_count_min"]),
path_count_max=int(data["path_count_max"]),
path_posterior_degrees=int(data["path_posterior_degrees"]),
path_posterior_scale=float(data["path_posterior_scale"]),
occupancy_alpha=float(data["occupancy_alpha"]),
occupancy_beta=float(data["occupancy_beta"]),
metadata=json.loads(str(data["metadata_json"].item())),
)
@cache
def _load_default_model() -> _OffshoreModel:
return _OffshoreModel.load(_DEFAULT_MODEL_PATH)
def _clip_halfplane(points, normal_x, normal_y, offset):
"""Clip a convex ring of ``(x, y)`` tuples to ``x @ normal <= offset``.
The ring is a few dozen points, so Python floats beat NumPy calls, and
the plain products round as NumPy's do.
"""
result = []
previous_x, previous_y = points[-1]
previous_value = previous_x * normal_x + previous_y * normal_y - offset
for x, y in points:
value = x * normal_x + y * normal_y - offset
previous_inside = previous_value <= 1e-12
inside = value <= 1e-12
if inside != previous_inside:
fraction = previous_value / (previous_value - value)
result.append(
(
previous_x + fraction * (x - previous_x),
previous_y + fraction * (y - previous_y),
)
)
if inside:
result.append((x, y))
previous_x, previous_y, previous_value = x, y, value
return result
def _polygon_from_support(support):
"""Intersect the half-planes of equally spaced support values, or ``None``."""
radius = 4.0 * float(np.max(np.abs(support)))
points = [
(-radius, -radius),
(radius, -radius),
(radius, radius),
(-radius, radius),
]
angles = np.arange(len(support)) * (2.0 * np.pi / len(support))
normals = np.column_stack([np.cos(angles), np.sin(angles)]).tolist()
for (normal_x, normal_y), offset in zip(normals, support.tolist(), strict=True):
points = _clip_halfplane(points, normal_x, normal_y, offset)
if len(points) < 3:
return None
polygon = Polygon(points)
if not polygon.is_valid or polygon.area <= np.finfo(float).eps:
return None
return orient(polygon, sign=1.0)
def _exact_vertices(coordinates, count):
"""Drop the least significant vertices of a ring until ``count`` remain.
Significance is the cross product of the two adjacent edges, twice the
area the removal loses. Returns ``None`` when no removal keeps the
polygon valid.
"""
coordinates = np.asarray(coordinates).copy()
while len(coordinates) > count:
previous = np.roll(coordinates, 1, axis=0)
following = np.roll(coordinates, -1, axis=0)
importance = np.abs(
(coordinates[:, 0] - previous[:, 0]) * (following[:, 1] - coordinates[:, 1])
- (coordinates[:, 1] - previous[:, 1])
* (following[:, 0] - coordinates[:, 0])
)
for index in np.argsort(importance):
candidate = np.delete(coordinates, index, axis=0)
polygon = Polygon(candidate)
if polygon.is_valid and polygon.area > 0:
coordinates = candidate
break
else:
return None
return coordinates
def _inscribed_share(corners: int) -> float:
"""Least share of a convex shape's area its largest ``corners``-gon covers.
The regular polygon in a circle attains it (Sas 1939), 0.41 for a
triangle and 0.64 for a quadrilateral.
"""
return corners / (2.0 * np.pi) * np.sin(2.0 * np.pi / corners)
def _reduce_open_path(coordinates, intermediate_count):
"""Drop the least significant interior vertices of an open path."""
coordinates = np.asarray(coordinates).copy()
target = intermediate_count + 2
while len(coordinates) > target:
previous = coordinates[:-2]
current = coordinates[1:-1]
following = coordinates[2:]
importance = np.abs(
(current[:, 0] - previous[:, 0]) * (following[:, 1] - current[:, 1])
- (current[:, 1] - previous[:, 1]) * (following[:, 0] - current[:, 0])
)
coordinates = np.delete(coordinates, int(np.argmin(importance)) + 1, axis=0)
return coordinates
def _sample_envelope(model, rng):
"""Draw a convex envelope with a sampled corner count, or ``None``."""
standardized = model.support_score_location + (
model.support_score_scale * rng.standard_normal(len(model.support_score_std))
)
scores = model.support_score_std * standardized
support = model.support_mean + scores @ model.support_components
envelope = _polygon_from_support(support)
if envelope is None:
return None
target = int(
np.rint(
np.exp(
model.hull_vertex_log_mean + rng.normal(scale=model.hull_vertex_log_std)
)
+ 2.0
)
)
if not model.hull_vertex_min <= target <= model.hull_vertex_max:
return None
coordinates = np.asarray(envelope.exterior.coords)[:-1]
if len(coordinates) > target:
reduced = _simplify_to_v(coordinates, target)
if reduced is not None:
reduced = _exact_vertices(np.asarray(reduced.exterior.coords)[:-1], target)
# Douglas-Peucker keeps the corners farthest from a chord, and at three
# corners kept as little as 0.01 % of the envelope, a sliver 18,000
# times longer than wide. Every convex shape holds a polygon of
# ``target`` corners with at least ``_inscribed_share`` of its area,
# so a smaller one is redone by dropping the corners that cost least.
if reduced is None or (
Polygon(reduced).area < _inscribed_share(target) * envelope.area
):
reduced = _exact_vertices(coordinates, target)
if reduced is None:
return None
coordinates = reduced
while len(coordinates) < target:
edges = np.roll(coordinates, -1, axis=0) - coordinates
index = int(np.argmax(np.linalg.norm(edges, axis=1)))
midpoint = 0.5 * (
coordinates[index] + coordinates[(index + 1) % len(coordinates)]
)
coordinates = np.insert(coordinates, index + 1, midpoint, axis=0)
polygon = Polygon(coordinates)
return orient(polygon, sign=1.0) if polygon.is_valid and polygon.area > 0 else None
def _sample_path(model, rng):
"""Draw an inward edge path in local coordinates from (0, 0) to (1, 0).
Returns ``None`` for a draw outside the fitted corner counts, or one that
is degenerate or self-intersecting.
"""
standardized = (
model.path_posterior_scale
* rng.standard_normal(len(model.path_feature_mean))
/ np.sqrt(
rng.chisquare(model.path_posterior_degrees) / model.path_posterior_degrees
)
)
features = model.path_feature_mean + model.path_feature_cholesky @ standardized
# A complex vector-matrix product through BLAS spins every BLAS thread for a
# few microseconds of work, so einsum forms it directly.
shape = model.path_mean + np.einsum("i,ij->j", features[:-2], model.path_components)
indentation_area = float(np.exp(features[-2]))
count = int(np.rint(np.exp(features[-1])))
if not model.path_count_min <= count <= model.path_count_max:
return None
x_coordinates = shape.real
profile = np.maximum(shape.imag, 0.0)
trial = np.column_stack([np.r_[0.0, x_coordinates, 1.0], np.r_[0.0, profile, 0.0]])
normalizer = float(Polygon(trial).area)
if normalizer <= np.finfo(float).eps:
return None
coordinates = np.column_stack(
[
np.r_[0.0, x_coordinates, 1.0],
np.r_[0.0, indentation_area * profile / normalizer, 0.0],
]
)
if not LineString(coordinates).is_simple:
return None
return _reduce_open_path(coordinates, count)
[docs]
class OffshorePolygon(Boundary):
"""Boundary sampled from the morphology of operating offshore wind farms.
Sampling is parametric and copies no reference polygon. The shape is scaled
to ``area`` and moved to ``center``. A heading is drawn per sample unless
``rotation`` pins it in degrees. Raises ``RuntimeError`` after
``max_attempts`` invalid draws.
"""
def __init__(
self,
*,
area: float | None = None,
center: Point2D | None = None,
rotation: float | None = None,
seed: int | np.random.Generator | None = None,
max_attempts: int = 256,
) -> None:
self.target_area = _optional_area(area)
self.center, self._positioned = _resolve_center(center)
self.rotation = (
None if rotation is None else _validate_finite("rotation", rotation)
)
self.seed = self._resolve_own_seed(seed)
self.max_attempts = _validate_integer("max_attempts", max_attempts)
def _build_geometry(self) -> BaseGeometry:
model = _load_default_model()
rng = as_rng(self.seed)
propensity = rng.beta(model.occupancy_alpha, model.occupancy_beta)
for _ in range(self.max_attempts):
envelope = _sample_envelope(model, rng)
if envelope is None:
continue
coordinates = np.asarray(envelope.exterior.coords)[:-1]
edges = np.roll(coordinates, -1, axis=0) - coordinates
lengths = np.linalg.norm(edges, axis=1)
irregular_count = int(rng.binomial(len(coordinates), propensity))
selected = set(
map(
int,
rng.choice(
len(coordinates),
size=irregular_count,
replace=False,
p=lengths / lengths.sum(),
),
)
)
result = []
failed = False
for index, start in enumerate(coordinates):
result.append(start)
if index not in selected:
continue
local = _sample_path(model, rng)
if local is None:
failed = True
break
edge = edges[index]
length = lengths[index]
direction = edge / length
inward = np.asarray([-direction[1], direction[0]])
result.extend(
start
+ local[1:-1, 0, None] * edge
+ local[1:-1, 1, None] * length * inward
)
if failed:
continue
polygon = Polygon(result)
if polygon.is_valid and polygon.area > np.finfo(float).eps:
return orient(_oriented(polygon, self.rotation, rng), sign=1.0)
raise RuntimeError(
f"could not generate a valid offshore polygon in {self.max_attempts} attempts"
)
