"""Boundary partition modifiers used by farms and clusters."""
from __future__ import annotations
from abc import ABC, abstractmethod
from collections.abc import Sequence
from copy import copy
import numpy as np
import shapely
from shapely import affinity
from shapely.errors import GEOSException
from shapely.geometry import LineString, MultiPoint, MultiPolygon, Point, Polygon
from shapely.geometry.base import BaseGeometry
from shapely.ops import split, unary_union
from .._random import as_rng, derive_seed, resolve_seed, spawn_seeds
from .._validation import _validate_integer, _validate_nonnegative, _validate_real_array
from ..strings import Straight, StringGuide
from .base import _translated_copies
from .modifiers import BoundaryModifier
__all__ = ["Partition"]
#: Largest accepted error of a fitted power-cell area, as a fraction of the parent.
_AREA_FRACTION_TOLERANCE = 0.01
#: Site draws tried when a power-diagram fit is locally singular.
_MAX_SITE_FIT_ATTEMPTS = 4
#: Area-fraction error at which the power-weight fit stops, far below the
#: accepted error so the cell areas follow their weights exactly.
_POWER_FIT_TOLERANCE = 1e-10
#: Newton steps and halvings of each step the power-weight fit tries.
_POWER_FIT_ITERATIONS = 100
_POWER_LINE_SEARCH_STEPS = 40
#: Relative change of the power-weight dual within its rounding, from areas
#: and second moments summed over the cells.
_DUAL_ROUNDING = 1e-12
#: Length, as a share of the region's span, that the power-weight fit adds
#: to every cell pair's shared edge in its curvature. An edge across a gap
#: or hole, or a cell emptied by a step, has no length inside the region,
#: so the areas stay flat as it moves, and this still points the step on.
#: Near the answer it changes a Newton step by about this share. Over 900
#: diagrams 0.001 fitted as fast as 0.0001 and 25 % faster than 0.01.
_PAIR_EDGE_SHARE = 0.001
#: Relative component area below which clipping artifacts are discarded.
_SPLINTER_AREA = 1e-12
#: Segments per quarter circle of the half gap grown around each cell. The
#: grown cell circumscribes its arcs, which widens the gap by a share of
#: ``1 / cos(pi / 128) - 1``, 3e-4.
_GAP_QUAD_SEGMENTS = 32
def _scaled_details(details: dict, factor: float, origin) -> dict:
"""Return a partition record scaled by ``factor`` about ``origin``.
A boundary scaled as a whole scales its cells, so the sites and channel
lines recorded beside them scale alike, and power weights, squared
distances, by the square of ``factor``.
"""
origin = np.asarray(origin, dtype=float)
def scaled(points):
return (origin + factor * (np.asarray(points, dtype=float) - origin)).tolist()
record = dict(details)
if "sites" in record:
record["sites"] = scaled(record["sites"])
if "channels" in record:
record["channels"] = [scaled(line) for line in record["channels"]]
if "power_weights_m2" in record:
weights = np.asarray(record["power_weights_m2"], dtype=float)
record["power_weights_m2"] = (factor**2 * weights).tolist()
return record
# Grids for cutting gaps, relative to the span of the partitioned boundary.
# Floating-point overlay fails where a vertex lies within rounding error of a
# cutter edge. Snap rounding on a grid coarser than the one power cells lie on
# separates them, while rounding onto that same grid fails often, so these
# stay well above ``_partition_precision``. A coarser grid is tried if needed.
_CUT_GRIDS = (1e-9, 1e-8, 1e-7)
def _span(geometry) -> float:
"""Return the bounding-box diagonal of ``geometry``, at least 1 m."""
min_x, min_y, max_x, max_y = geometry.bounds
return max(float(np.hypot(max_x - min_x, max_y - min_y)), 1.0)
def _partition_precision(geometry) -> float:
"""Return a scale-aware grid that makes shared power-cell edges identical."""
return _span(geometry) * 1e-12
def _polygon_components(geometry) -> list[Polygon]:
"""Return every genuine positive-area polygon component."""
if isinstance(geometry, Polygon):
parts = [geometry]
elif isinstance(geometry, MultiPolygon):
parts = list(geometry.geoms)
else:
parts = [
part for part in shapely.get_parts(geometry) if isinstance(part, Polygon)
]
return [part for part in parts if not part.is_empty and part.area > 0.0]
def _polygons(geometry) -> list[Polygon]:
"""Return polygon components of an overlay result, dropping splinters.
Components below ``_SPLINTER_AREA`` of the largest are overlay artifacts.
"""
parts = _polygon_components(geometry)
if len(parts) < 2:
return parts
largest = max(part.area for part in parts)
return [part for part in parts if part.area >= _SPLINTER_AREA * largest]
#: Share of one farm's area below which a piece of the boundary joins its
#: nearest neighbor rather than taking a farm of its own.
_MIN_PIECE_SHARE = 0.5
def _farm_components(geometry, n_farms: int) -> list[Polygon | MultiPolygon]:
"""The polygon components of ``geometry``, merged until each holds a farm.
A setback or exclusion zone can split a boundary into more pieces than
farms, or pinch off a piece too small for one. The smallest group joins
the one nearest it while there are more groups than farms or it holds
under :data:`_MIN_PIECE_SHARE` of a farm's area, so a farm takes whole
neighboring pieces as one cell. A sliver pinched off had taken a farm
of its own, and the cut its neighbor needed vanished.
"""
groups = _polygon_components(geometry)
total = sum(group.area for group in groups)
while len(groups) > 1:
areas = [group.area for group in groups]
smallest = int(np.argmin(areas))
enough = areas[smallest] * n_farms >= _MIN_PIECE_SHARE * total
if len(groups) <= n_farms and enough:
break
piece = groups.pop(smallest)
distances = shapely.distance(piece, np.asarray(groups, dtype=object))
nearest = int(np.argmin(distances))
groups[nearest] = MultiPolygon(
[*_polygon_components(groups[nearest]), *_polygon_components(piece)]
)
return groups
def _allocate_components(total: int, geometries: list) -> np.ndarray:
"""Split ``total`` cells across ``geometries`` by area, at least one each.
``geometries`` holds at most ``total`` entries, as :func:`_farm_components`
gives. Largest-remainder rounding, with ties broken by order.
"""
areas = np.asarray([geometry.area for geometry in geometries], dtype=float)
quotas = total * areas / areas.sum()
counts = np.floor(quotas).astype(int)
remainder = total - int(counts.sum())
if remainder:
# Round fractional quotas for deterministic ties.
fractions = np.round(quotas - counts, 9)
order = np.argsort(-fractions, kind="stable")
counts[order[:remainder]] += 1
# Ensure each positive component receives one cell.
for missing in np.flatnonzero(counts == 0):
donors = np.flatnonzero(counts > 1)
surplus = np.round(counts[donors] - quotas[donors], 9)
donor = donors[int(np.argmax(surplus))]
counts[donor] -= 1
counts[missing] = 1
return counts
def _sample_points(geometry: Polygon, count: int, rng) -> np.ndarray:
"""Draw ``count`` uniform points strictly inside ``geometry`` by rejection."""
min_x, min_y, max_x, max_y = geometry.bounds
points: list[tuple[float, float]] = []
batch = max(32, 4 * count)
attempts = 0
while len(points) < count:
candidates = rng.uniform((min_x, min_y), (max_x, max_y), size=(batch, 2))
inside = shapely.contains_xy(geometry, candidates[:, 0], candidates[:, 1])
points.extend(map(tuple, candidates[inside]))
attempts += batch
if attempts > 100_000 and len(points) < count:
raise ValueError("could not sample Voronoi sites inside the boundary")
return np.asarray(points[:count], dtype=float)
def _sites_by_component(components, sites: np.ndarray) -> list[np.ndarray]:
"""``sites`` grouped by the component nearest each, every one given one.
Raises ``ValueError`` when a component is left without a site.
"""
distances = shapely.distance(
np.asarray(components, dtype=object)[:, None],
shapely.points(sites)[None, :],
)
nearest = np.argmin(distances, axis=0)
grouped = [sites[nearest == index] for index in range(len(components))]
if any(not len(group) for group in grouped):
raise ValueError("the carried sites leave a part of the boundary without one")
return grouped
def _validate_sites(sites, n_farms: int) -> np.ndarray | None:
"""Return ``sites`` as ``n_farms`` distinct finite points, or ``None``."""
if sites is None:
return None
values = _validate_real_array("sites", sites)
if values.shape != (n_farms, 2) or not np.all(np.isfinite(values)):
raise ValueError("sites must contain one finite (x, y) point per farm")
if len(np.unique(values, axis=0)) != len(values):
raise ValueError("sites must be distinct")
return values
def _voronoi_cells(
geometry: Polygon, sites: np.ndarray
) -> list[Polygon | MultiPolygon]:
"""Return the Voronoi cell of each site clipped to ``geometry``, in site order."""
if len(sites) == 1:
return [geometry]
envelope = geometry.envelope.buffer(_span(geometry))
diagram = shapely.voronoi_polygons(MultiPoint(sites), extend_to=envelope)
candidates = list(shapely.get_parts(diagram))
cells: list[Polygon | MultiPolygon] = []
for site in sites:
point = Point(site)
matches = [cell for cell in candidates if cell.covers(point)]
if len(matches) != 1:
raise ValueError("Voronoi construction did not produce one cell per site")
clipped = matches[0].intersection(geometry)
polygons = _polygons(clipped)
if not polygons:
raise ValueError("Voronoi construction produced an empty cell")
cell = unary_union(polygons)
if not cell.is_valid or cell.area <= 0:
raise ValueError("Voronoi construction produced an invalid cell")
cells.append(cell)
return cells
def _gap_cut(cells, extent, gap: float) -> list:
"""Each of ``cells`` without the points within half of ``gap`` of another.
Two cells then lie at least ``gap`` apart wherever the boundary lies
between them, along a shared edge, around a corner and where an edge
meets the boundary, and the cut moves continuously as the cells do. Each
grown cell circumscribes its arcs, so no kept point lies closer than half
the gap. Shared edges extended across the boundary had cut cells
kilometers from any edge, and jumped as the boundary grew.
"""
if gap == 0 or len(cells) < 2:
return list(cells)
radius = 0.5 * gap / np.cos(np.pi / (4 * _GAP_QUAD_SEGMENTS))
grown = [cell.buffer(radius, quad_segs=_GAP_QUAD_SEGMENTS) for cell in cells]
cut = []
for index, cell in enumerate(cells):
near = [
other
for position, other in enumerate(grown)
if position != index and other.intersects(cell)
]
if near:
cut.extend(_subtract([cell], shapely.union_all(near), extent))
else:
cut.append(cell)
return cut
def _subtract(geometries, cutter, extent) -> list:
"""Return each geometry minus ``cutter`` by snap-rounded overlay.
The cutter grows by two grid cells first, since rounding the cutter and
rounding the cut each move a vertex by less than one, so a cut is never
narrower than asked. Raises ``ValueError`` when the overlay fails on every
grid, like the other partition failures.
"""
span = _span(extent)
for relative in _CUT_GRIDS:
grid = span * relative
widened = cutter.buffer(2.0 * grid, join_style="mitre")
try:
return [
shapely.difference(geometry, widened, grid_size=grid)
for geometry in geometries
]
except GEOSException:
continue
raise ValueError("the gap could not be cut from the partition cells")
def _cleaned_cells(cells, name: str) -> list:
"""Return ``cells`` without splinters, refusing one the gap removed or broke."""
cleaned_cells = []
for cell in cells:
polygons = _polygons(cell)
if not polygons:
raise ValueError(f"gap is too wide for a {name} cell, which it removes")
cleaned = unary_union(polygons)
if not cleaned.is_valid or cleaned.area <= 0:
raise ValueError(f"gap produces an invalid {name} cell")
cleaned_cells.append(cleaned)
return cleaned_cells
class _PartitionModifier(BoundaryModifier, ABC):
"""A boundary modifier that splits a boundary and keeps its cells apart.
``n_farms`` is the exact number of cells, ``gap`` the distance in meters
kept between neighboring cells, and ``seed`` drives any random choice.
The gap stays in meters when the boundary is resized. ``details`` holds a
JSON-compatible record of the last partition.
"""
#: Whether the partition is cut around sites, which it samples if not given.
_has_sites: bool = False
#: Whether the sites were carried from another size of the boundary.
_carried: bool = False
def __init__(
self,
n_farms: int,
*,
gap: float = 0.0,
seed: int | np.random.Generator | None = None,
) -> None:
self.n_farms = _validate_integer("n_farms", n_farms)
self.gap = _validate_nonnegative("gap", gap)
self._pinned_seed = seed is not None
self.seed = resolve_seed(seed)
self.details: dict = {}
@property
def _in_meters(self) -> bool:
return self.gap > 0.0
@property
def _samples_sites(self) -> bool:
"""Whether the partition draws its own sites."""
return self._has_sites and getattr(self, "sites", None) is None
def _carrying(self, sites: np.ndarray) -> "_PartitionModifier":
"""A copy cutting its cells around ``sites`` carried from another size.
A resized boundary samples its sites once, at the size its base is
given, and carries them with the base, so its cells grow with it.
Sampled again at each size tried, a site near an exclusion zone was
kept at one size and rejected at the next, and the area jumped. A
carried site may lie in a zone the base grew around, where it still
bounds its cell, and each joins the component nearest it.
"""
carried = copy(self)
setattr(carried, "sites", np.asarray(sites, dtype=float))
carried._carried = True
return carried
@abstractmethod
def partition(self, geometry: BaseGeometry) -> tuple[Polygon | MultiPolygon, ...]:
"""Return the ``n_farms`` cells of ``geometry`` in a deterministic order."""
def apply(self, geometry: BaseGeometry) -> BaseGeometry:
"""Return the union of the cells ``partition`` retains.
Cells cut with no gap tile ``geometry``, which is returned as it is.
Use ``partition`` to keep individual cells.
"""
cells = self.partition(geometry)
return geometry if self.gap == 0 else unary_union(cells)
def __repr__(self) -> str:
# The classes are public only through the ``Partition`` namespace.
return f"Partition.{super().__repr__()}"
[docs]
class Voronoi(_PartitionModifier):
"""Partition a boundary into Voronoi cells clipped to it.
Each cell holds the points nearer its site than any other.
Parameters
----------
n_farms:
Exact number of partition cells.
sites:
Optional finite ``(x, y)`` site per farm, fixed in the plane, so a
resized boundary must still contain them. Sites are sampled uniformly
inside the boundary when omitted.
gap:
Distance in meters kept between cells, half taken from each.
seed:
Random seed used when sampling sites.
"""
_has_sites = True
def __init__(
self,
n_farms: int,
*,
sites: Sequence[Sequence[float]] | np.ndarray | None = None,
gap: float = 0.0,
seed: int | np.random.Generator | None = None,
) -> None:
super().__init__(n_farms, gap=gap, seed=seed)
self.sites = _validate_sites(sites, self.n_farms)
@property
def _in_meters(self) -> bool:
return super()._in_meters or self.sites is not None
[docs]
def partition(self, geometry: BaseGeometry) -> tuple[Polygon | MultiPolygon, ...]:
"""Return Voronoi cells clipped to ``geometry``, one per site.
Sites in a ``MultiPolygon`` parent are allocated across components by
area, at least one each. With more components than farms the smallest
joins its nearest until each farm has one. Holes are preserved. Each
cell loses the points within ``gap / 2`` of another. Sampled sites are redrawn
a few times when the gap removes a whole cell, and given sites or a
last draw that still loses one raise ``ValueError``. Sampled sites
replay from ``seed`` and the cell order is deterministic. Sites and
the draws taken are recorded in ``details``.
"""
components = _farm_components(geometry, self.n_farms)
if not components:
raise ValueError("partition boundary must contain polygonal geometry")
counts = _allocate_components(self.n_farms, components)
rng = as_rng(self.seed)
if self.sites is not None and self._carried:
component_sites = _sites_by_component(components, self.sites)
elif self.sites is not None:
if not np.all(
shapely.intersects_xy(geometry, self.sites[:, 0], self.sites[:, 1])
):
raise ValueError("sites must lie inside the boundary")
component_sites = []
for component in components:
mask = shapely.intersects_xy(
component, self.sites[:, 0], self.sites[:, 1]
)
component_sites.append(self.sites[mask])
if [len(values) for values in component_sites] != counts.tolist():
raise ValueError(
"sites must follow the area allocation across disconnected components"
)
def cells_of(component_sites):
cells_found = []
for component, sites in zip(components, component_sites):
cells = _voronoi_cells(component, sites)
cells = _gap_cut(cells, component, self.gap)
cells_found.extend(_cleaned_cells(cells, "Voronoi"))
return cells_found
site_attempts = 1
if self.sites is not None:
all_cells = cells_of(component_sites)
else:
# Two sites drawn close together bound a small cell, which a gap
# can remove whole, so sampled sites are redrawn a few times.
for site_attempts in range(1, _MAX_SITE_FIT_ATTEMPTS + 1):
component_sites = [
_sample_points(component, int(count), rng)
for component, count in zip(components, counts)
]
try:
all_cells = cells_of(component_sites)
except ValueError:
if site_attempts == _MAX_SITE_FIT_ATTEMPTS:
raise
else:
break
all_sites = [site for sites in component_sites for site in sites.tolist()]
self.details = {
"method": "voronoi",
"n_farms": self.n_farms,
"gap": self.gap,
"seed": self.seed,
"site_attempts": site_attempts,
"sites": all_sites,
}
return tuple(all_cells)
def _clipped_convex(
vertices: np.ndarray,
labels: list[int],
normal: np.ndarray,
cutoff: float,
label: int,
) -> tuple[np.ndarray, list[int]]:
"""Clip the convex ring ``vertices`` to ``normal @ point <= cutoff``.
Sutherland-Hodgman clipping. ``labels[k]`` names the line the edge from
vertex ``k`` lies on, and an edge along the clipping line gets ``label``.
A vertex on the line is kept once, so no vertex repeats.
"""
values = vertices @ normal - cutoff
inside = values <= 0.0
if inside.all():
return vertices, labels
if not inside.any():
return vertices[:0], []
def crossing(index: int, following: int) -> np.ndarray:
weight = values[index] / (values[index] - values[following])
return vertices[index] + weight * (vertices[following] - vertices[index])
clipped = []
clipped_labels = []
for index in range(len(vertices)):
following = (index + 1) % len(vertices)
if inside[index] and inside[following]:
clipped.append(vertices[index])
clipped_labels.append(labels[index])
elif inside[index] and values[index] < 0.0:
clipped += [vertices[index], crossing(index, following)]
clipped_labels += [labels[index], label]
elif inside[index]:
clipped.append(vertices[index])
clipped_labels.append(label)
elif inside[following] and values[following] < 0.0:
clipped.append(crossing(index, following))
clipped_labels.append(labels[index])
return np.asarray(clipped), clipped_labels
def _power_regions(
geometry: Polygon, sites: np.ndarray, weights: np.ndarray
) -> list[tuple[np.ndarray, list[int]]]:
r"""Return each site's convex power region around ``geometry``.
The region of site ``i`` is a box around ``geometry`` clipped to the
half-planes where
:math:`\lVert x-p_i\rVert^2-w_i \le \lVert x-p_j\rVert^2-w_j` for every
other site ``j``. Each region is its ring and the site index ``j`` of
each edge, ``-1`` for an edge of the box. An empty region has no vertices.
"""
min_x, min_y, max_x, max_y = geometry.bounds
margin = max(max_x - min_x, max_y - min_y, 1.0)
box = np.asarray(
[
[min_x - margin, min_y - margin],
[max_x + margin, min_y - margin],
[max_x + margin, max_y + margin],
[min_x - margin, max_y + margin],
]
)
squared_norms = np.einsum("ij,ij->i", sites, sites)
regions = []
for index, site in enumerate(sites):
vertices, labels = box, [-1] * 4
for other_index, other_site in enumerate(sites):
if index == other_index:
continue
normal = 2.0 * (other_site - site)
if not normal.any():
raise ValueError("power-diagram sites must be distinct")
cutoff = (
squared_norms[other_index]
- squared_norms[index]
+ weights[index]
- weights[other_index]
)
vertices, labels = _clipped_convex(
vertices, labels, normal, float(cutoff), other_index
)
if len(vertices) < 3:
vertices, labels = vertices[:0], []
break
regions.append((vertices, labels))
return regions
def _power_cells(geometry: Polygon, regions) -> list[Polygon | MultiPolygon | None]:
"""Return each power region clipped to ``geometry``, or ``None`` where empty."""
polygons = np.asarray(
[Polygon(vertices) if len(vertices) else Polygon() for vertices, _ in regions],
dtype=object,
)
cells: list[Polygon | MultiPolygon | None] = []
for cell in shapely.intersection(geometry, polygons):
parts = _polygons(cell)
cells.append(unary_union(parts) if parts else None)
return cells
def _facet_lengths(geometry: Polygon, regions) -> np.ndarray:
"""Return the length inside ``geometry`` of the edge shared by each region pair.
Each edge is measured once, from the region of the lower site index.
"""
pairs = []
segments = []
for index, (vertices, labels) in enumerate(regions):
for position, other_index in enumerate(labels):
if other_index > index:
following = (position + 1) % len(vertices)
pairs.append((index, other_index))
segments.append(LineString([vertices[position], vertices[following]]))
lengths = np.zeros((len(regions), len(regions)))
if pairs:
rows, columns = np.asarray(pairs).T
measured = shapely.length(shapely.intersection(geometry, segments))
lengths[rows, columns] = measured
lengths[columns, rows] = measured
return lengths
def _area_fractions(cells, area: float) -> np.ndarray:
"""Return each cell's share of ``area``, counting a missing cell as zero."""
return np.asarray([0.0 if cell is None else cell.area / area for cell in cells])
def _second_moment(geometry, point) -> float:
"""The integral over ``geometry`` of the squared distance to ``point``.
Oriented rings run counterclockwise round each shell and clockwise round
each hole, so the signed sum over every edge subtracts the holes.
"""
rings = shapely.get_rings(shapely.get_parts(shapely.orient_polygons(geometry)))
coordinates, ring = shapely.get_coordinates(rings, return_index=True)
joined = ring[1:] == ring[:-1]
x0, y0 = (coordinates[:-1][joined] - point).T
x1, y1 = (coordinates[1:][joined] - point).T
cross = x0 * y1 - x1 * y0
return float(
np.sum(cross * (x0 * x0 + x0 * x1 + x1 * x1 + y0 * y0 + y0 * y1 + y1 * y1))
/ 12.0
)
def _fit_power_cells(
geometry: Polygon,
sites: np.ndarray,
target_fractions: np.ndarray,
):
"""Fit power weights so cell areas match ``target_fractions``.
Returns the cells, the weights in square meters and a JSON-compatible fit
record. Raises ``ValueError`` when a cell stays empty or an area misses its
target by more than ``_AREA_FRACTION_TOLERANCE``.
"""
if len(sites) == 1:
return (
[geometry],
np.zeros(1),
{
"success": True,
"evaluations": 0,
"max_area_fraction_error": 0.0,
},
)
# Translate to the centroid before power fitting to preserve precision.
origin = np.asarray(geometry.centroid.coords[0], dtype=float)
local_geometry = affinity.translate(
geometry, xoff=-float(origin[0]), yoff=-float(origin[1])
)
local_sites = sites - origin
span = max(
local_geometry.bounds[2] - local_geometry.bounds[0],
local_geometry.bounds[3] - local_geometry.bounds[1],
1.0,
)
weight_scale = float(span**2)
def weights_from_free(free_weights: np.ndarray) -> np.ndarray:
"""Return square-meter weights for ``free_weights``, pinning the last to zero.
The diagram is invariant to a common shift of every weight, so one site
is fixed while the rest are scaled by the squared boundary span.
"""
return np.append(np.asarray(free_weights, dtype=float), 0.0) * weight_scale
area = local_geometry.area
targets = np.asarray(target_fractions, dtype=float)
def evaluate(free_weights: np.ndarray):
"""The dual value, cell area fractions and power regions at the weights.
The dual of fitting the areas, the sum over cells of the integral of
the power distance plus each weight times its target area, is concave
in the weights for any region, also one in pieces or with holes, and
its gradient is each target area less the cell's area.
"""
weights = weights_from_free(free_weights)
regions = _power_regions(local_geometry, local_sites, weights)
cells = _power_cells(local_geometry, regions)
fractions = _area_fractions(cells, area)
moments = sum(
_second_moment(cell, site)
for cell, site in zip(cells, local_sites)
if cell is not None
)
value = (
moments / (area * weight_scale)
+ weights @ (targets - fractions) / weight_scale
)
return value, fractions, regions
distances = np.linalg.norm(local_sites[:, None] - local_sites[None], axis=-1)
np.fill_diagonal(distances, np.inf)
free = np.zeros(len(sites) - 1)
value, fractions, regions = evaluate(free)
evaluations = 1
converged = False
for _ in range(_POWER_FIT_ITERATIONS):
errors = targets - fractions
if np.max(np.abs(errors)) <= _POWER_FIT_TOLERANCE:
converged = True
break
lengths = _facet_lengths(local_geometry, regions) + _PAIR_EDGE_SHARE * span
rates = lengths / (2.0 * distances)
curvature = (np.diag(rates.sum(axis=1)) - rates)[:-1, :-1]
step = np.linalg.solve(curvature * weight_scale / area, errors[:-1])
slope = float(errors[:-1] @ step)
scale = 1.0
for _ in range(_POWER_LINE_SEARCH_STEPS):
trial = free + scale * step
trial_value, trial_fractions, trial_regions = evaluate(trial)
evaluations += 1
# Near the answer the dual gains less than its rounding, so a step
# that shrinks the largest area error without losing more than
# rounding counts as progress too. Taken while the dual fell, such
# steps had cycled through layouts with four cells empty.
if trial_value >= value + 1e-4 * scale * slope or (
trial_value >= value - _DUAL_ROUNDING * (1.0 + abs(value))
and np.max(np.abs(targets - trial_fractions)) < np.max(np.abs(errors))
):
break
scale /= 2.0
else:
break
free, value, fractions, regions = (
trial,
trial_value,
trial_fractions,
trial_regions,
)
weights = weights_from_free(free)
local_cells = _power_cells(
local_geometry, _power_regions(local_geometry, local_sites, weights)
)
fractions = _area_fractions(local_cells, local_geometry.area)
maximum_error = float(np.max(np.abs(fractions - target_fractions)))
if any(cell is None for cell in local_cells) or (
maximum_error > _AREA_FRACTION_TOLERANCE
):
raise ValueError(
"could not fit nonempty power-diagram cells to the requested area weights"
)
cells = [
None
if cell is None
else affinity.translate(cell, xoff=float(origin[0]), yoff=float(origin[1]))
for cell in local_cells
]
return (
cells,
weights,
{
"success": converged,
"evaluations": evaluations,
"max_area_fraction_error": maximum_error,
},
)
[docs]
class PowerDiagram(_PartitionModifier):
r"""Partition one boundary into cells with controlled relative areas.
Cell ``i`` contains points where its weighted squared distance
.. math:: d_i(x) = \lVert x-p_i\rVert^2-w_i
is no larger than for any other site, the power diagram of Aurenhammer
(1987, `doi:10.1137/0216006 <https://doi.org/10.1137/0216006>`_).
Additive weights ``w_i`` are fitted so cell areas follow ``area_weights``
before any gap is removed.
Parameters
----------
n_farms:
Exact number of partition cells.
area_weights:
Positive relative target areas, one per farm. Equal by default.
sites:
Optional finite ``(x, y)`` site per farm. Sites are sampled inside the
boundary when omitted.
gap:
Distance in meters kept between cells, half taken from each.
seed:
Random seed used when sampling sites.
"""
_has_sites = True
def __init__(
self,
n_farms: int,
*,
area_weights: Sequence[float] | None = None,
sites: Sequence[Sequence[float]] | np.ndarray | None = None,
gap: float = 0.0,
seed: int | np.random.Generator | None = None,
) -> None:
super().__init__(n_farms, gap=gap, seed=seed)
self.sites = _validate_sites(sites, self.n_farms)
if area_weights is None:
values = np.ones(self.n_farms, dtype=float)
else:
values = _validate_real_array("area_weights", area_weights)
if values.shape != (self.n_farms,):
raise ValueError("area_weights must contain one value per farm")
if not np.all(np.isfinite(values)) or np.any(values <= 0):
raise ValueError("area_weights must be finite and positive")
self.area_fractions = values / values.sum()
@property
def _in_meters(self) -> bool:
return super()._in_meters or self.sites is not None
[docs]
def partition(self, geometry: BaseGeometry) -> tuple[Polygon | MultiPolygon, ...]:
"""Return power cells of ``geometry``, one per site.
Holes are honored. A parent in several disconnected polygons, as a
setback or exclusion zone can leave, is fitted as one region, so the
areas stay exact and a cell may take parts of several polygons.
Unseeded sites are redrawn a few times before a failed area fit
propagates. Each cell loses the points within ``gap / 2`` of another, and a
gap that consumes a whole cell raises ``ValueError``. Cell order is
deterministic. ``details`` records the sites, target, achieved and
retained area fractions, and the fitted weights in
``details["power_weights_m2"]``.
"""
components = _polygon_components(geometry)
if not components:
raise ValueError("partition boundary must contain polygonal geometry")
# GIS data often stores one connected polygon as a MultiPolygon.
geometry = components[0] if len(components) == 1 else MultiPolygon(components)
if self.sites is None:
rng = as_rng(self.seed)
for site_attempts in range(1, _MAX_SITE_FIT_ATTEMPTS + 1):
sites = _sample_points(geometry, self.n_farms, rng)
try:
fitted_cells, weights, fit = _fit_power_cells(
geometry, sites, self.area_fractions
)
except ValueError:
if site_attempts == _MAX_SITE_FIT_ATTEMPTS:
raise
else:
break
else:
site_attempts = 1
sites = self.sites.copy()
inside = shapely.intersects_xy(geometry, sites[:, 0], sites[:, 1])
if not self._carried and not np.all(inside):
raise ValueError("sites must lie inside the boundary")
# Duplicate supplied sites are rejected by ``_validate_sites``.
fitted_cells, weights, fit = _fit_power_cells(
geometry, sites, self.area_fractions
)
cells = [cell for cell in fitted_cells if cell is not None]
# Apply one common precision grid before finding shared cell edges.
precision = _partition_precision(geometry)
cells = [shapely.set_precision(cell, precision) for cell in cells]
achieved_fractions = [float(cell.area / geometry.area) for cell in cells]
cells = _gap_cut(cells, geometry, self.gap)
cleaned_cells = _cleaned_cells(cells, "power-diagram")
self.details = {
"method": "power_diagram",
"n_farms": self.n_farms,
"gap": self.gap,
"seed": self.seed,
"site_attempts": site_attempts,
"sites": sites.tolist(),
"target_area_fractions": self.area_fractions.tolist(),
"achieved_area_fractions": achieved_fractions,
"retained_area_fractions": [
float(cell.area / geometry.area) for cell in cleaned_cells
],
"power_weights_m2": weights.tolist(),
"fit": fit,
}
return tuple(cleaned_cells)
def _normalised_guide(guide: StringGuide) -> LineString:
"""Return the guide moved and turned to run from (0, 0) to (1, 0)."""
geometry = guide.geometry
if not isinstance(geometry, LineString):
raise ValueError("channel guide must contain one continuous string")
# Which side of a self-crossing cut lies below it is undefined, and every
# cut along one failed.
if not geometry.is_simple:
raise ValueError("channel guide must not cross itself")
coordinates = np.asarray(geometry.coords, dtype=float)
chord = coordinates[-1] - coordinates[0]
length = float(np.linalg.norm(chord))
if length == 0:
raise ValueError("channel guide endpoints must be distinct")
angle = -np.degrees(np.arctan2(chord[1], chord[0]))
transformed = affinity.translate(
geometry, xoff=-coordinates[0, 0], yoff=-coordinates[0, 1]
)
transformed = affinity.rotate(transformed, angle, origin=(0, 0))
transformed = affinity.scale(transformed, 1.0 / length, 1.0 / length, origin=(0, 0))
return transformed
def _guide_frame(
template: LineString, geometry
) -> tuple[LineString, float, np.ndarray]:
"""Return the guide stretched across ``geometry``, its scale span and centroid.
The halfway point of the guide lies on the x axis, so a cut through the
centroid follows the guide where it is curved rather than its chord.
"""
span = _span(geometry)
line = affinity.scale(template, 2.5 * span, 2.5 * span, origin=(0, 0))
line = affinity.translate(line, xoff=-1.25 * span)
# An arc of more than about 140 degrees bows so far from its chord that a
# cut along the chord through the cell never meets the curve.
bow = line.interpolate(0.5, normalized=True).y
if bow != 0.0:
line = affinity.translate(line, yoff=-bow)
return line, span, np.asarray(geometry.centroid.coords[0])
def _sides(geometry, line: LineString, region: Polygon) -> list | None:
"""Return the parts of ``geometry`` inside and outside ``region``, or ``None``.
``geometry`` is split along ``line``, the edge of ``region`` within it, and
each face goes to the side holding it. Splinter faces are dropped as in
``_polygons``. A side made of several faces is a ``MultiPolygon``. Returns
``None`` when the overlay fails, a side is empty, or a side is invalid.
"""
try:
faces = shapely.get_parts(split(geometry, line))
except (GEOSException, ValueError):
return None
faces = faces[shapely.get_type_id(faces) == shapely.GeometryType.POLYGON]
areas = shapely.area(faces)
faces = faces[areas >= _SPLINTER_AREA * areas.max(initial=0.0)]
if not len(faces):
return None
inside = shapely.contains(region, shapely.point_on_surface(faces))
if inside.all() or not inside.any():
return None
sides = [shapely.union_all(faces[inside]), shapely.union_all(faces[~inside])]
if not all(side.is_valid for side in sides):
return None
return sides
def _channel(lines, gap: float):
"""Return the band ``gap`` wide along each of ``lines``, with flat ends."""
return shapely.buffer(
lines, gap / 2.0, quad_segs=16, cap_style="flat", join_style="mitre"
)
def _region_below(line: LineString, span: float) -> Polygon:
"""Return the region under ``line``, closed far below and beside it.
``line`` runs left to right in its own frame, as ``_guide_frame`` returns
it. The closing edges lie more than ``span`` from the frame origin, so
within that distance the region's edge is ``line`` alone.
"""
coordinates = np.asarray(line.coords, dtype=float)
margin = 2.0 * span
left = coordinates[:, 0].min() - margin
right = coordinates[:, 0].max() + margin
low = coordinates[:, 1].min() - margin
start, end = coordinates[0], coordinates[-1]
closing = [(right, end[1]), (right, low), (left, low), (left, start[1])]
return Polygon([*map(tuple, coordinates), *closing])
[docs]
class Channels(_PartitionModifier):
"""Partition a boundary by recursive cuts along a string guide.
Each cut halves the number of farms left in a cell. Its orientation is
drawn at random and its offset gives that share of the area exactly. The
parts of a cell on one side of the cut form one new cell, so a cut across
a bay of the boundary leaves a cell in several parts. For a guide that is
a graph over its chord the cut moves continuously with the cell. The cuts
are chosen before any gap is removed, so a resized boundary is cut alike
and the area left changes smoothly with its size. Each cell then loses a
band ``gap`` wide along every cut that split one of its parents.
Parameters
----------
n_farms:
Exact number of partition cells.
guide:
``StringGuide`` with one continuous string and distinct endpoints
that does not cross itself, stretched across the cell to shape each
cut. ``Straight()`` by default.
A random guide given no seed is redrawn from ``seed``, so every
reseeded partition bends its channels anew. A seeded guide keeps its
shape.
gap:
Width in meters of the channel each cut removes.
seed:
Random seed for the cut orientations and an unseeded guide.
Raises
------
ValueError
If ``guide`` is not one continuous string with distinct endpoints, if
it crosses itself, if no orientation gives a cut with two valid sides,
or if the gap is too wide for a cell.
"""
def __init__(
self,
n_farms: int,
*,
guide: StringGuide | None = None,
gap: float = 0.0,
seed: int | np.random.Generator | None = None,
) -> None:
super().__init__(n_farms, gap=gap, seed=seed)
if guide is None:
guide = Straight()
if not isinstance(guide, StringGuide):
raise TypeError("guide must be a StringGuide")
# Reject a guide of several strings here rather than at the first cut.
_normalised_guide(guide)
self.guide = guide
def _realized_guide(self) -> StringGuide:
"""Return the guide the cuts follow, an unseeded one drawn from ``seed``."""
# The guide's own unseeded draw would stay fixed while a sampler
# reseeds the partition, giving every farm the same bend.
return self.guide._reseeded(derive_seed(self.seed, "guide"))
def _best_cut(self, geometry, template: LineString, target_fraction: float, rng):
"""Return the two sides of the cut and its line, the first side below it.
The cut's orientation is drawn at random, and its offset along the
normal is root-found so the cell's area below the guide is exactly
``target_fraction``. The faces the cut leaves go to the side holding
them, so a side may have several parts. For a guide that is a graph
over its chord, such as ``Straight`` or ``Organic``, the root is
unique and the cut moves continuously with the cell. A later of 32
orientations is tried only when the overlay fails or leaves an empty
or invalid side.
"""
# SciPy loads only when a partition needs it, keeping ``import oowl`` light.
from scipy.optimize import brentq
base, span, centroid = _guide_frame(template, geometry)
region = _region_below(base, span)
heights = np.asarray(base.coords, dtype=float)[:, 1]
area = geometry.area
points = shapely.get_coordinates(geometry)
def excess(offset: float, below: Polygon, normal: np.ndarray) -> float:
moved = _translated_copies(below, centroid + offset * normal)[0]
return shapely.intersection(geometry, moved).area / area - target_fraction
for angle in rng.uniform(0.0, 180.0, size=32):
angle = float(angle)
radians = np.deg2rad(angle)
normal = np.asarray([-np.sin(radians), np.cos(radians)])
rotated = affinity.rotate(base, angle, origin=(0, 0), use_radians=False)
below = affinity.rotate(region, angle, origin=(0, 0), use_radians=False)
reach = (points - centroid) @ normal
try:
# Below the lowest offset the region misses the cell, above
# the highest it holds the whole cell.
offset = brentq(
excess,
reach.min() - heights.max(),
reach.max() - heights.min(),
args=(below, normal),
xtol=1e-12 * span,
)
shift = centroid + offset * normal
line = _translated_copies(rotated, shift)[0]
sides = _sides(geometry, line, _translated_copies(below, shift)[0])
except GEOSException:
sides = None
if sides is not None:
return sides, line
raise ValueError("could not find a channel that splits the boundary in two")
def _partition_component(
self, geometry, template: LineString, n_farms: int, seed: int
):
"""Return the cells and channel lines of one connected component."""
pending: list[tuple[BaseGeometry, int, int, tuple[int, ...]]] = [
(geometry, n_farms, seed, ())
]
cells = []
# The indices in ``channels`` of the cuts that split each cell's parents.
ancestry = []
channels = []
while pending:
cell, leaves, child_seed, cuts = pending.pop()
if leaves == 1:
cells.append(cell)
ancestry.append(cuts)
continue
left_leaves = leaves // 2
right_leaves = leaves - left_leaves
pieces, line = self._best_cut(
cell, template, left_leaves / leaves, as_rng(child_seed)
)
channels.append(line)
cuts = (*cuts, len(channels) - 1)
child_seeds = spawn_seeds(child_seed, 2)
pending.append((pieces[1], right_leaves, child_seeds[1], cuts))
pending.append((pieces[0], left_leaves, child_seeds[0], cuts))
if self.gap == 0 or not channels:
return cells, channels
# A band is removed only from the cells its cut split, so it never
# reaches into a neighbor across an earlier cut.
bands = _channel(channels, self.gap)
cells = [
_subtract([cell], shapely.union_all(bands[list(cuts)]), geometry)[0]
for cell, cuts in zip(cells, ancestry)
]
return _cleaned_cells(cells, "Channels"), channels
[docs]
def partition(self, geometry: BaseGeometry) -> tuple[Polygon | MultiPolygon, ...]:
"""Return the cells cut from ``geometry`` by recursive channel splits.
Cells in a ``MultiPolygon`` parent are allocated across components by
area, at least one each. With more components than farms the smallest
joins its nearest until each farm has one. Holes survive the cuts.
Cells and their order replay from
``seed``. The guide and the channel lines are recorded in ``details``.
"""
components = _farm_components(geometry, self.n_farms)
if not components:
raise ValueError("partition boundary must contain polygonal geometry")
counts = _allocate_components(self.n_farms, components)
seeds = spawn_seeds(self.seed, len(components))
guide = self._realized_guide()
template = _normalised_guide(guide)
cells = []
channels = []
for component, count, seed in zip(components, counts, seeds):
local_cells, local_channels = self._partition_component(
component, template, int(count), seed
)
cells.extend(local_cells)
channels.extend(local_channels)
self.details = {
"method": "channels",
"n_farms": self.n_farms,
"gap": self.gap,
"seed": self.seed,
"guide": guide.metadata(),
"channels": [list(map(list, line.coords)) for line in channels],
}
return tuple(cells)
class Partition:
"""Namespace for concrete boundary-partition modifiers."""
Channels = Channels
PowerDiagram = PowerDiagram
Voronoi = Voronoi
def __new__(cls, *args, **kwargs):
raise TypeError(
"use Partition.Channels, Partition.PowerDiagram, or Partition.Voronoi"
)
