weavingspace.tiling_utils

   1"""Utility functions for the `weavingspace` module.
   2
   3Most are geometric convenience functions for commonly applied
   4operations.
   5"""
   6
   7from __future__ import annotations
   8
   9import re
  10import string
  11from typing import TYPE_CHECKING, NamedTuple
  12
  13import geopandas as gpd
  14import numpy as np
  15import pandas as pd
  16import shapely.affinity as affine
  17import shapely.geometry as geom
  18import shapely.ops
  19from shapely.algorithms import polylabel
  20
  21if TYPE_CHECKING:
  22  from collections.abc import Iterable
  23
  24PRECISION = 6
  25RESOLUTION = 1e-6
  26
  27def _parse_strand_label(s:str) -> list[str]:
  28  """Break a strand label specification in to a list of labels.
  29
  30  Args:
  31    s (str): see get_strand_ids() for details
  32
  33  Returns:
  34    list[str]: list of strand labels.
  35
  36  """
  37  # replace multiple open and close parentheses with single parentheses
  38  clean_s = re.sub("[(]+", "(", re.sub("[)]+", ")", s))
  39  is_combo = False
  40  output = []
  41  current = ""
  42  for c in clean_s:
  43    if is_combo:
  44      if c == ")":
  45        # end of a combination
  46        output.append(current)
  47        current = ""
  48        is_combo = False
  49      else:
  50        current = current + c
  51    elif c == "(":
  52      # start of a combination
  53      is_combo = True
  54    else:
  55      output.append(c)
  56  return output
  57
  58
  59def get_strand_ids(strands_spec: str) -> tuple[list[str]]:
  60  """Convert strands specification string to list of lists of strand labels.
  61
  62  String format is "a|bc|(de)f" where | separates strands in each
  63  direction and () designates combining labels into a single strand that
  64  will be sliced lengthwise. Example output:
  65
  66        "a|bc|(de)f" -> (["a"], ["b", "c"], ["de", "f"])
  67
  68  Superflous parentheses are removed, but no other error-checks are
  69  applied.
  70
  71  Args:
  72    strands_spec (str): the strands specification to be parsed.
  73
  74  Returns:
  75    tuple[str]: tuple of lists of labels for each set of strands.
  76
  77  """
  78  strand_ids = [_parse_strand_label(s) for s in strands_spec.split("|")]
  79  strand_ids = strand_ids if len(strand_ids) == 3 else [*strand_ids, [""]]
  80  return tuple(strand_ids)
  81
  82
  83def get_regular_polygon(spacing:float, n:int) -> geom.Polygon:
  84  r"""Return regular n-gon centered on (0, 0) with height specified.
  85
  86  n-gon has horizontal base and height given by spacing. Generation is based on
  87  equally spaced points around a circle with a radius calculated slightly
  88  differently, depending on whether the number of sides is even or odd.
  89
  90            Even         Odd
  91        \  /    \  /      |
  92         \/   ___\/___    |
  93         /\      /\      / \
  94        /  \    /  \    /   \
  95
  96  The returned polygon is gridified. Note also that shapely silently reorders
  97  the coordinates clockwise from the lower left corner, regardless of the CCW
  98  order in which we supply them.
  99
 100  Args:
 101    spacing (_type_): required height.
 102    n (int): number of sides.
 103
 104  Returns:
 105    geom.Polygon: required geom.Polygon.
 106
 107  """
 108  if n % 2 == 0:
 109    R = spacing / np.cos(np.pi / n) / 2
 110  else:
 111    R = spacing / (1 + np.cos(np.pi / n))
 112  angle_0 = -np.pi / 2 + np.pi / n
 113  a_diff = 2 * np.pi / n
 114  angles = [angle_0 + a_diff * i for i in range(n)]
 115  corners = [(R * np.cos(angle), R * np.sin(angle)) for angle in angles]
 116  return geom.Polygon(corners)
 117
 118
 119def get_corners(shape:geom.Polygon,
 120                repeat_first:bool = True) -> list[geom.Point]:
 121  """Return list of geom.Points on boundary of a polygon.
 122
 123  The first point is optionally repeated. No simplification is carried out
 124  (e.g. if a line segment has a 'corner' along its length, it is NOT removed;
 125  see get_clean_polygon for that). Points have precision set to the
 126  packagedefault tiling_utils.RESOLUTION.
 127
 128  Args:
 129    shape (geom.Polygon): polygon whose corners are required.
 130    repeat_first (bool, optional): if True the first corner is repeated in the
 131      returned list, if False it is omitted. Defaults to True.
 132
 133  Returns:
 134    list[geom.Point]: list of geom.Point vertices of the polygon.
 135
 136  """
 137  corners = [geom.Point(pt) for pt in shape.exterior.coords]
 138  return corners if repeat_first else corners[:-1]
 139
 140
 141def get_sides(shape:geom.Polygon) -> list[geom.LineString]:
 142  """Return polygon sides as list of geom.LineStrings.
 143
 144  Resolution is set to the package default tiling_utils.RESOLUTION. No
 145  simplification for successive colinear sides is carried out.
 146
 147  Args:
 148    shape (geom.Polygon): polygon whose edges are required.
 149
 150  Returns:
 151    list[geom.LineString]: list of geom.LineString sides of the polygon.
 152
 153  """
 154  corners = get_corners(shape)
 155  return [geom.LineString([p1, p2])
 156          for p1, p2 in zip(corners[:-1], corners[1:], strict = True)]
 157
 158
 159def get_side_lengths(shape:geom.Polygon) -> list[float]:
 160  """Return list of lengths of polygon sides.
 161
 162  No simplification for corners along sides is carried out.
 163
 164  Args:
 165    shape (geom.Polygon): polygon whose edge lengths are required.
 166
 167  Returns:
 168    list[float]: list of side lengths.
 169
 170  """
 171  corners = get_corners(shape)
 172  return [
 173    p1.distance(p2) for p1, p2 in zip(corners[:-1], corners[1:], strict = True)]
 174
 175
 176def get_side_bearings(shape:geom.Polygon) -> tuple[float,...]:
 177  """Return list of angles between sides of polygon and the positive x-axis.
 178
 179  Angle is in degrees and calculated when proceeding from the first point in
 180  each side to its end point. This should usually be CW around the polygon.
 181
 182  Args:
 183    shape (geom.Polygon): polygon whose side bearings are required.
 184
 185  Returns:
 186    tuple[float,...]: tuple of bearings of each edge.
 187
 188  """
 189  return tuple([np.degrees(np.arctan2(
 190    e.coords[1][1] - e.coords[0][1], e.coords[1][0] - e.coords[0][0]))
 191    for e in get_sides(shape)])
 192
 193
 194def get_interior_angles(shape:geom.Polygon) -> list[float]:
 195  """Return angles (in degrees) between successive edges of a polygon.
 196
 197  No polygon simplification is carried out so some angles may be 180 (i.e. a
 198  'corner' along a side, such that successive sides are colinear).
 199
 200  Args:
 201    shape (geom.Polygon): polygon whose angles are required.
 202
 203  Returns:
 204    list[float]: list of angles.
 205
 206  """
 207  corners = get_corners(shape, repeat_first = False)
 208  wrap_corners = corners[-1:] + corners + corners[:1]
 209  triples = zip(wrap_corners[:-2],
 210                wrap_corners[1:-1],
 211                wrap_corners[2:], strict = True)
 212  return [get_inner_angle(p1, p2, p3) for p1, p2, p3 in triples]
 213
 214
 215def get_clean_polygon(
 216      shape:geom.MultiPolygon|geom.Polygon|gpd.GeoSeries,
 217    ) ->  geom.MultiPolygon|geom.Polygon|gpd.GeoSeries:
 218  """Return polygon with successive co-linear or very close corners removed.
 219
 220  Particularly useful for tidying polygons weave tilings that have been
 221  assembled from multiple 'cells' in the weave grid.
 222
 223  Args:
 224    shape (geom.MultiPolygon|geom.Polygon): polygon to clean.
 225
 226  Returns:
 227    geom.Polygon: cleaned polygon.
 228
 229  """
 230  if isinstance(shape, geom.MultiPolygon):
 231    return geom.MultiPolygon([get_clean_polygon(p) for p in shape.geoms])
 232  if isinstance(shape, geom.GeometryCollection):
 233    elements = [p for p in shape.geoms if isinstance(p, geom.Polygon)]
 234    if len(elements) > 1:
 235      return get_clean_polygon(geom.MultiPolygon(elements))
 236    return get_clean_polygon(elements[0])
 237  corners = get_corners(shape, repeat_first = False)
 238  # first remove any 'near neighbour' corners
 239  distances = get_side_lengths(shape)
 240  to_remove = [np.isclose(d, 0, rtol = RESOLUTION, atol = 10 * RESOLUTION)
 241                for d in distances]
 242  corners = [c for c, r in zip(corners, to_remove, strict = True) if not r]
 243  # next remove any that are colinear
 244  p = geom.Polygon(corners)
 245  corners = get_corners(p, repeat_first = False)
 246  angles = get_interior_angles(p)
 247  to_remove = [np.isclose(a,   0, rtol = RESOLUTION, atol = RESOLUTION) or
 248               np.isclose(a, 180, rtol = RESOLUTION, atol = RESOLUTION)
 249               for a in angles]
 250  corners = [c for c, r in zip(corners, to_remove, strict = True) if not r]
 251  return gridify(geom.Polygon(corners))
 252
 253
 254def rotate_preserving_order(
 255    polygon:geom.Polygon,
 256    angle:float,
 257    centre:geom.Point = None,
 258  ) -> geom.Polygon:
 259  """Return supplied polygon rotated with order of corners preserved.
 260
 261  Order of polygon corners is not guaranteed to be preserved by shapely.
 262  affinity.rotate).
 263
 264  Args:
 265    polygon (geom.Polygon): polygon to rotate.
 266    angle (float): desired angle of rotation (in degrees).
 267    centre (geom.Point): the rotation centre (passed on to shapely.affinity.
 268      rotate).
 269
 270  Returns:
 271    geom.Polygon: rotated polygon.
 272
 273  """
 274  ctr = "center" if centre is None else centre
 275  corners = get_corners(polygon, repeat_first = False)
 276  return geom.Polygon([affine.rotate(c, angle, ctr) for c in corners])
 277
 278
 279def ensure_cw(shape:geom.Polygon) -> geom.Polygon:
 280  """Return the polygon with its outer boundary vertices in clockwise order.
 281
 282  It is important to note that shapely.set_precision() imposes clockwise order
 283  on polygons, and since it is used widely throughout theses modules, it makes
 284  sense to impose this order.
 285
 286  Args:
 287    shape (geom.Polygon): the polygon.
 288
 289  Returns:
 290    geom.Polygon: the polygon in clockwise order.
 291
 292  """
 293  if geom.LinearRing(shape.exterior.coords).is_ccw:
 294    return shape.reverse()
 295  return shape
 296
 297
 298def get_inner_angle(p1:geom.Point, p2:geom.Point, p3:geom.Point) -> float:
 299  r"""Return angle between line p1-p2 and p2-p3, i.e., the angle A below.
 300
 301            p2
 302           / \
 303          / A \
 304        p1     p3
 305
 306  Args:
 307    p1 (geom.Point): first point.
 308    p2 (geom.Point): second 'corner' point.
 309    p3 (geom.Point): third point.
 310
 311  Returns:
 312    float: angle in degrees.
 313
 314  """
 315  return np.degrees(np.arctan2(p3.y - p2.y, p3.x - p2.x) -
 316                    np.arctan2(p1.y - p2.y, p1.x - p2.x)) % 360
 317
 318
 319def is_regular_polygon(shape:geom.Polygon) -> bool:
 320  """Test if supplied polygon is regular (i.e. equal sides and angles).
 321
 322  Args:
 323    shape (geom.Polygon): polygon to test.
 324
 325  Returns:
 326    bool: True if polygon is regular, False if not.
 327
 328  """
 329  side_lengths = get_side_lengths(shape)
 330  angles = get_interior_angles(shape)
 331  return all(np.isclose(side_lengths, side_lengths[0])) \
 332     and all(np.isclose(angles, angles[0]))
 333
 334
 335def is_tangential(shape:geom.Polygon) -> bool:
 336  """Determine if supplied polygon is tangential.
 337
 338  A tangential polygon can have a circle inscribed tangential to all its sides.
 339
 340  NOTE: not currently used but retained in case incentre is restored to use
 341  'home made' code rather than polylabel.polylabel.
 342
 343  Note that this will fail for polygons with successive colinear sides, meaning
 344  that polygons should be fully simplified...
 345
 346  Args:
 347    shape (geom.Polygon): polygon to test.
 348
 349  Returns:
 350    bool: True if polygon is tangential, False if not.
 351
 352  """
 353  if is_regular_polygon(shape):
 354    return True
 355  side_lengths = get_side_lengths(shape)
 356  n = len(side_lengths)
 357  if n % 2 == 1:
 358    if n == 3: #triangles are easy!
 359      return True
 360    # odd number of sides there is a nice solvable system  of equations see
 361    # https://math.stackexchange.com/questions/4065370/tangential-polygons-conditions-on-edge-lengths
 362    #
 363    # TODO: this is not reliable for reasons I don't fully understand. There is
 364    # a corresponding crude catch exception in incentre... but this needs a
 365    # more complete investigation
 366    mat = np.identity(n) + np.roll(np.identity(n), 1, axis = 1)
 367    fractions = np.linalg.inv(mat) @ side_lengths / side_lengths
 368    if not(np.isclose(np.mean(fractions),
 369                      0.5, rtol= RESOLUTION, atol = RESOLUTION)):
 370      return False
 371    ones = [np.isclose(f, 1, rtol = RESOLUTION, atol = RESOLUTION)
 372            for f in fractions]
 373    zeros = [np.isclose(f, 0, rtol = RESOLUTION, atol = RESOLUTION)
 374             for f in fractions]
 375    negatives = [f < 0 for f in fractions]
 376    return not (any(ones) or any(zeros) or any(negatives))
 377  if n == 4:
 378    # for quadrilaterals odd and even side lengths are equal
 379    return np.isclose(sum(side_lengths[0::2]),
 380                      sum(side_lengths[1::2]), rtol = RESOLUTION)
 381  # other even numbers of sides... hard to avoid brute force
 382  bisectors = [get_angle_bisector(shape, i) for i in range(n)]
 383  intersects = [b1.intersection(b2) for b1, b2 in
 384                zip(bisectors, bisectors[1:] + bisectors[:1], strict = False)
 385                if b1.intersects(b2)]
 386  distances = [i1.distance(i2) for i1, i2 in
 387               zip(intersects, intersects[1:] + intersects[:1],
 388                   strict = False)]
 389  return all([d <= 2 * RESOLUTION for d in distances])  # noqa: C419
 390
 391
 392def get_incentre(shape:geom.Polygon) -> geom.Point:
 393  """Get a polygon incentre.
 394
 395  See https://en.wikipedia.org/wiki/Incenter.
 396
 397  This method relies on the polygon being tangential, i.e. there is an
 398  inscribed circle to which all sides of the polygon are tangent. It will work
 399  on all the polygons encountered in the Laves tilings, but is not guaranteed
 400  to work on all polygons.
 401
 402  Given that the polygon is tangential, the radius of the inscribed circle is
 403  the [apothem of the polygon](https://en.wikipedia.org/wiki/Apothem) given by
 404  2 x Area / Perimeter. We apply a parallel offset of this size to two sides of
 405  the polygon and find their intersection to determine the centre of the circle
 406  givng the incentre of the polygon.
 407
 408  Args:
 409    shape (geom.Polygon): the polygon.
 410
 411  Returns:
 412    geom.Point: the incentre of the polygon.
 413
 414  """
 415  if is_regular_polygon(shape):
 416    return shape.centroid
 417  return polylabel.polylabel(shape)
 418
 419
 420def get_incircle(shape:geom.Polygon) -> geom.Polygon:
 421  """Return incircle of supplied polygon.
 422
 423  Args:
 424    shape(geom.Polygon): polygon for which incircle is required.
 425
 426  Returns:
 427    geom.Polygon: a polygon representation of the incircle
 428
 429  """
 430  c = get_incentre(shape)
 431  r = min([c.distance(e) for e in get_sides(shape)])
 432  return c.buffer(r)
 433
 434
 435def get_apothem(shape:geom.Polygon) -> float:
 436  """Return length of the apothem of supplied polygon.
 437
 438  Args:
 439    shape(geom.Polygon): polygon for which apothem length is required.
 440
 441  Returns:
 442    float: an estimate of the apothem length.
 443
 444  """
 445  if is_regular_polygon(shape):
 446    return 2 * shape.area / geom.LineString(shape.exterior.coords).length
 447  c = polylabel.polylabel(shape, 1)
 448  return min([c.distance(e) for e in get_sides(shape)])
 449
 450
 451def get_angle_bisector(shape:geom.Polygon, v = 0) -> geom.LineString:
 452  """Return angle bisector of specified corner of supplied polygon.
 453
 454  Args:
 455    shape (geom.Polygon): the polygon
 456    v (int, optional): index of the corner whose bisector is required.
 457    Defaults to 0.
 458
 459  Returns:
 460    geom.LineString: line which bisects the specified corner.
 461
 462  """
 463  pts = get_corners(shape, repeat_first = False)
 464  n = len(pts)
 465  p0 = pts[v % n]
 466  p1 = pts[(v + 1) % n]
 467  p2 = pts[(v - 1) % n]
 468  d01 = p0.distance(p1)
 469  d02 = p0.distance(p2)
 470  if d01 < d02:
 471    p2 = geom.LineString([p0, p2]).interpolate(d01)
 472  else:
 473    p1 = geom.LineString([p0, p1]).interpolate(d02)
 474  m = geom.Point((p1.x + p2.x) / 2, (p1.y + p2.y) / 2)
 475  # we have to scale this to much larger in case of an angle near 180
 476  return affine.scale(geom.LineString([p0, m]), 10, 10, origin = p0)
 477
 478
 479def _get_interior_vertices(tiles:gpd.GeoDataFrame) -> gpd.GeoSeries:
 480  """Return points not on outer boundary of supplied set of tiles.
 481
 482  Args:
 483    tiles (gpd.GeoDataFrame): a set of polygons.
 484
 485  Returns:
 486    gpd.GeoSeries: interior vertices of the set of polygons.
 487
 488  """
 489  tiles:gpd.GeoSeries = gridify(tiles.geometry)
 490  uu:geom.Polygon = safe_union(tiles, as_polygon = True)
 491  interior_pts = set()
 492  for tile in tiles:
 493    for corner in tile.exterior.coords:
 494      if uu.contains(geom.Point(corner).buffer(1e-3, cap_style = "square")):
 495        interior_pts.add(corner)
 496  return gridify(gpd.GeoSeries([geom.Point(p) for p in interior_pts]))
 497
 498
 499def gridify(
 500    gs:gpd.GeoSeries | gpd.GeoDataFrame | shapely.Geometry,
 501  ) -> gpd.GeoSeries | gpd.GeoDataFrame | shapely.Geometry:
 502  """Return supplied GeoSeries rounded to tiling_utils.RESOLUTION.
 503
 504  Args:
 505    gs (gpd.GeoSeries): geometries to gridify.
 506    precision (int, optional): digits of precision. Defaults to 6.
 507
 508  Returns:
 509    gpd.GeoSeries: the rounded geometries.
 510
 511  """
 512  if isinstance(gs, shapely.Geometry) :
 513    return shapely.set_precision(gs, grid_size = RESOLUTION)
 514  if isinstance(gs, gpd.GeoSeries):
 515    return gpd.GeoSeries(
 516      [shapely.set_precision(g, grid_size = RESOLUTION) for g in list(gs)],
 517      crs = gs.crs)
 518  if isinstance(gs, gpd.GeoDataFrame):
 519    gs.geometry = gpd.GeoSeries(
 520      [shapely.set_precision(g, grid_size = RESOLUTION)
 521       for g in list(gs.geometry)], crs = gs.crs)
 522    return gs
 523  return None
 524
 525
 526def get_dual_tile_unit(unit) -> gpd.GeoDataFrame:
 527  """Convert supplied TileUnit to candidate GeoDataFrame of its dual TileUnit.
 528
 529  NOTE: this is complicated and not remotely guaranteed to work!
 530  a particular issue is that where to place the vertices of the faces
 531  of the dual with respect to the tiles in the original is ill-defined.
 532  This is because the dual process is topologically not metrically defined,
 533  so that exact vertex locations are ambiguous. Tiling duality is defined in
 534  Section 4.2 of Grunbaum B, Shephard G C, 1987 _Tilings and Patterns_ (W. H.
 535  Freeman and Company, New York)
 536
 537  NOTE: In general, this method will work only if all supplied tiles are
 538  regular polygons. A known exception is if the only non-regular polygons are
 539  triangles.
 540
 541  NOTE: 'clean' polygons are required. If supplied polygons have messy
 542  vertices with multiple points where there is only one proper point, bad
 543  things are likely to happen! Consider using `clean_polygon()` on the
 544  tile geometries.
 545
 546  Because of the above limitations, we only return a GeoDataFrame
 547  for inspection. However some `weavingspace.tile_unit.TileUnit` setup
 548  methods in `weavingspace.tiling_geometries` use this method, where we are
 549  confident the returned dual is valid.
 550
 551  Args:
 552    unit (TileUnit): the tiling for which the dual is required.
 553
 554  Returns:
 555    gpd.GeoDataFrame: GeoDataFrame that could be the tiles attribute for
 556      a TileUnit of the dual tiling.
 557
 558  """
 559  # get a local patch of this Tiling
 560  local_patch:gpd.GeoDataFrame = unit.get_local_patch(r = 1, include_0 = True)
 561  # Find the interior points of these tiles - these will be guaranteed
 562  # to have a sequence of surrounding tiles incident on them
 563  interior_pts = _get_interior_vertices(local_patch)
 564  # Compile a list of the polygons incident on the interior points
 565  cycles = [{poly_id for poly_id, p in enumerate(local_patch.geometry)
 566             if pt.distance(p) < RESOLUTION * 2} for pt in interior_pts]
 567  dual_faces = []
 568  for cycle in cycles:
 569    ids, pts = [], []
 570    for poly_id in cycle:
 571      ids.append(local_patch.tile_id[poly_id])
 572      poly = local_patch.geometry[poly_id]
 573      pts.append(get_incentre(poly))
 574    # sort them into CW order so they are well formed
 575    sorted_coords = _sort_cw([(pt.x, pt.y, ID)
 576                              for pt, ID in zip(pts, ids, strict = True)])
 577    dual_faces.append(
 578      (geom.Polygon([(pt_id[0], pt_id[1]) for pt_id in sorted_coords]),
 579       "".join([pt_id[2] for pt_id in sorted_coords])))
 580  # ensure the resulting face centroids are inside the original tile
 581  # displaced a little to avoid uncertainties at corners/edges
 582  # This is a guess at something that _could_ work, which turned out well!
 583  dual_faces = [(f, ID) for f, ID in dual_faces
 584                if affine.translate(
 585                  unit.get_prototile_from_vectors().loc[0, "geometry"],
 586                  RESOLUTION * 10, RESOLUTION * 10).contains(f.centroid)]
 587  gdf = gpd.GeoDataFrame(
 588    data = {"tile_id": [f[1] for f in dual_faces]}, crs = unit.crs,
 589    geometry = gpd.GeoSeries([f[0] for f in dual_faces]))
 590  # ensure no duplicates
 591  gdf = gdf.dissolve(by = "tile_id", as_index = False).explode(
 592    index_parts = False, ignore_index = True)
 593  gdf.tile_id = _relabel(gdf.tile_id)
 594  return gdf
 595
 596
 597def _relabel(data:Iterable) -> list:
 598  """Return supplied data reassigned with unique values from ascii_letters.
 599
 600  Args:
 601    data (Iterable): the data to relabel
 602
 603  Returns:
 604    list: the reassigned data
 605
 606  """
 607  new_data = {}
 608  d_count = 0
 609  for d in data:
 610    if d not in new_data:
 611      new_data[d] = string.ascii_letters[d_count]
 612      d_count = d_count + 1
 613  return [new_data[d] for d in data]
 614
 615
 616def _sort_cw(
 617  pts_ids:list[tuple[float, float, str]]) -> list[tuple[float, float, str]]:
 618  """Sort supplied x,y,ID tuple into CW order relative to their mean centre.
 619
 620  Args:
 621    pts_ids (list[tuple[float, float, str]]): A tuple of a pair of
 622      floats and a string.
 623
 624  Returns:
 625    list: a list in the same format as supplied sorted into
 626      clockwise order of the point locations.
 627
 628  """
 629  x = [p[0] for p in pts_ids]
 630  y = [p[1] for p in pts_ids]
 631  cx, cy = np.mean(x), np.mean(y)
 632  dx = [_ - cx for _ in x]
 633  dy = [_ - cy for _ in y]
 634  angles = [np.arctan2(dy, dx) for dx, dy in zip(dx, dy, strict = True)]
 635  d = dict(zip(angles, pts_ids, strict = True))
 636  return [pt_id for angle, pt_id in sorted(d.items(), reverse = True)]
 637
 638
 639def write_map_to_layers(gdf:gpd.GeoDataFrame,
 640                        fname:str = "output.gpkg",
 641                        tile_var:str = "tile_id") -> None:
 642  """Write supplied GeoDataFrame to a GPKG file.
 643
 644  Geopackage file layers will be based on the `tile_var` attribute.
 645
 646  Args:
 647    gdf (gpd.GeoDataFrame): the GeoDataFrame.
 648    fname (str, optional): filename to write.
 649    tile_var (str, optional): the attribute to use to separate
 650      output file into layers. Defaults to "tile_id".
 651
 652  """
 653  grouped = gdf.groupby(tile_var, as_index = False)
 654  for e in pd.Series.unique(gdf[tile_var]):
 655    grouped.get_group(e).to_file(fname, layer = e, driver = "GPKG")
 656
 657
 658def get_collapse_distance(geometry:geom.Polygon) -> float:
 659  """Return distance under which polygon will disappear if negatively buffered.
 660
 661  Performs a binary search between an upper bound based on the radius of circle
 662  of equal area to the polygon, and 0.
 663
 664  Args:
 665    geometry (geom.Polygon): the polygon.
 666
 667  Returns:
 668    float: its collapse distance.
 669
 670  """
 671  if is_regular_polygon(geometry):
 672    return get_apothem(geometry)
 673  radius = np.sqrt(geometry.area / np.pi)
 674  lower = 0
 675  upper = radius
 676  delta = 1e12
 677  stop_delta = radius / 1000
 678  while delta > stop_delta:
 679    new_r = (lower + upper) / 2
 680    if geometry.buffer(-new_r).area > 0:
 681      lower = new_r
 682      delta = upper - new_r
 683    else:
 684      upper = new_r
 685      delta = new_r - lower
 686  return new_r
 687
 688
 689def get_largest_polygon(polygons:gpd.GeoSeries) -> gpd.GeoSeries:
 690  """Return largest polygon in a GeoSeries as a GeoSeries of one polygon.
 691
 692  NOT CURRENTLY USED ANYWHERE IN MODULE.
 693
 694  Args:
 695    polygons (gpd.GeoSeries): the set of polygons to pick from.
 696
 697  Returns:
 698    gpd.GeoSeries: the largest polygon.
 699
 700  """
 701  actual_polygons = [p for p in polygons.geometry
 702                     if isinstance(p, (geom.Polygon, geom.MultiPolygon))]
 703  areas = sorted(actual_polygons, key = lambda x: x.area, reverse = True)
 704  return gpd.GeoSeries(areas[:1])
 705
 706
 707def touch_along_an_edge(p1:geom.Polygon, p2:geom.Polygon) -> bool:
 708  """Test if two polygons touch along an edge.
 709
 710  Checks that the intersection area of the two polygons buffered by
 711  a small amount is large enough to indicate that they neighbour at more
 712  than a corner.
 713
 714  Args:
 715    p1 (geom.Polygon): First polygon
 716    p2 (geom.Polygon): Second polygon
 717
 718  Returns:
 719    bool: True if they neighbour along an edge
 720
 721  """
 722  return (
 723    p1.buffer(RESOLUTION, join_style = "mitre", cap_style = "square")
 724      .intersection(p2.buffer(RESOLUTION, join_style = "mitre",
 725                                          cap_style = "square"))
 726      .area > 16 * RESOLUTION ** 2)
 727
 728
 729def get_width_height_left_bottom(
 730    gs:gpd.GeoSeries,
 731  ) -> tuple[float,float,float,float]:
 732  """Return width, height, left and bottom limits of a GeoSeries.
 733
 734  Args:
 735    gs (geopandas.GeoSeries): GeoSeries for which limits are required.
 736
 737  Returns:
 738    tuple: four float values of width, height, left and bottom of gs.
 739
 740  """
 741  bb = gs.total_bounds
 742  return (bb[2] - bb[0], bb[3] - bb[1], bb[0], bb[1])
 743
 744
 745def get_bounding_ellipse(
 746    shapes:gpd.GeoSeries,
 747    mag:float = 1.0,
 748  ) -> gpd.GeoSeries:
 749  """Return an ellipse containing the supplied shapes.
 750
 751  The method used is to calculate the size of square that would contain
 752  the shapes, if they had an aspect ratio 1, then stretch the circle in
 753  the x, y directions according to the actual aspect ratio of the shapes.
 754
 755  Args:
 756    shapes (gpd.GeoSeries): the shapes to be contained.
 757    mag (float, optional): optionally increase the size of the returned
 758      ellipse by this scale factor. Defaults to 1.0.
 759
 760  Returns:
 761    gpd.GeoSeries: the set of shapes.
 762
 763  """
 764  width, height, left, bottom = get_width_height_left_bottom(shapes)
 765  c = geom.Point(left + width / 2, bottom + height / 2)
 766  r = min(width, height) * np.sqrt(2)
 767  circle = [c.buffer(r)]
 768  return gpd.GeoSeries(circle, crs = shapes.crs).scale(
 769    width / r * mag / np.sqrt(2), height / r * mag / np.sqrt(2), origin = c)
 770
 771
 772def get_tiling_edges(tiles:gpd.GeoSeries) -> gpd.GeoSeries:
 773  """Return linestring GeoSeries from supplied polygon GeoSeries.
 774
 775  This is used to allow display of edges of tiles in legend when they are
 776  masked by an ellipse (if we instead clip polygons then the ellipse edge
 777  will also show in the result.)
 778
 779  Args:
 780    tiles (gpd.GeoSeries): Polygons to convert.
 781
 782  Returns:
 783    gpd.GeoSeries: LineStrings from the supplied Polygons.
 784
 785  """
 786  tiles = tiles.explode(ignore_index = True)
 787  edges = [geom.LineString(p.exterior.coords) for p in tiles]
 788  return gpd.GeoSeries(edges, crs = tiles.crs)
 789
 790
 791def get_polygon_sector(
 792    shape:geom.Polygon,
 793    start:float = 0.0,
 794    end:float = 1.0,
 795  ) -> geom.Polygon:
 796  """Return a sector of the provided Polygon.
 797
 798  The returned sector is a section of the polygon boundary between the
 799  normalized start and end positions, and including the polygon centroid.
 800  Should (probably) only be applied to convex polygons.
 801
 802  Args:
 803    shape (geom.Polygon): the Polygon.
 804    start (float): normalized start position along the boundary. Defaults to
 805      0.
 806    end (float): normalized start position along the boundary. Defaults to
 807      1.
 808
 809  Returns:
 810    geom.Polygon: the requested polygon sector.
 811
 812  """
 813  if start == end:
 814    # must return a null polygon since the calling context
 815    # expects to get something back... which most likely
 816    # is needed to align with other data
 817    return geom.Polygon()
 818  if start * end < 0:
 819    # either side of 0/1 so assume required sector includes 0
 820    e1, e2 = min(start, end), max(start, end)
 821    arc1 = shapely.ops.substring(geom.LineString(shape.exterior.coords),
 822                                 np.mod(e1, 1), 1, normalized = True)
 823    arc2 = shapely.ops.substring(geom.LineString(shape.exterior.coords),
 824                                 0, e2, normalized = True)
 825    sector = geom.Polygon(
 826      [shape.centroid] + list(arc1.coords) + list(arc2.coords)[1:])
 827  else:
 828    arc = shapely.ops.substring(geom.LineString(shape.exterior.coords),
 829                                start, end, normalized = True)
 830    sector = geom.Polygon([shape.centroid, *list(arc.coords)])
 831  return sector
 832
 833
 834def repair_polygon(
 835    polygon:geom.Polygon|geom.MultiPolygon|gpd.GeoSeries,
 836    shrink_then_grow:bool = True,
 837  ) -> geom.Polygon|gpd.GeoSeries:
 838  """Repair a polyon or GeoSeries by double buffering.
 839
 840  This is method is often unofficially recommended (on stackexchange etc.)
 841  even in the shapely docs, to resolve topology issues and extraneous
 842  additional vertices appearing when spatial operations are repeatedly
 843  applied.
 844
 845  Optionally the buffer may be applied in the opposite order (i.e. grow then
 846  shrink). This operation may also convert a MultiPolygon that has some 'stray'
 847  parts to a Polygon.
 848
 849  Args:
 850    polygon(geom.Polygon|gpd.GeoSeries): Polygon or GeoSeries to clean.
 851    res (float, optional): buffer size to use. Defaults to 1e-3.
 852    shrink_then_grow (bool, optional): if True the negative buffer is
 853      applied first, otherwise the buffer operations are applied in
 854      reverse. Defaults to True.
 855
 856  Returns:
 857    geom.Polygon|gpd.GeoSeries: the cleaned Polygon or GeoSeries.
 858
 859  """
 860  if type(polygon) is gpd.GeoSeries:
 861    return gpd.GeoSeries([repair_polygon(p, shrink_then_grow)
 862                          for p in polygon]).set_crs(polygon.crs)
 863  if polygon is None:
 864    return None
 865  if shrink_then_grow:
 866    return polygon.buffer(
 867      -RESOLUTION * 10, join_style = "mitre", cap_style = "square").buffer(
 868       RESOLUTION * 10, join_style = "mitre", cap_style = "square")
 869  return polygon.buffer(
 870     RESOLUTION * 10, join_style = "mitre", cap_style = "square").buffer(
 871    -RESOLUTION * 10, join_style = "mitre", cap_style = "square")
 872
 873
 874def safe_union(
 875    gs:gpd.GeoSeries,
 876    as_polygon:bool = False,
 877  ) -> gpd.GeoSeries|geom.Polygon:
 878  """Union supplied GeoSeries buffering to avoid gaps and internal edges.
 879
 880  Optionally returns a Polygon or a GeoSeries.
 881
 882  NOTE: APPEARS NOT TO BE VERY SAFE... find workarounds for the few occasions
 883  on which it is used.
 884
 885  Frequently when unioning polygons that are ostensibly adjacent 'rogue'
 886  internal boundaries remain in the result. We can avoid this by buffering the
 887  polygons before unioning them, then reversing the buffer on the unioned
 888  shape.
 889
 890  Args:
 891    gs (gpd.GeoSeries): the Polygons to union.
 892    res (float, optional): size of the buffer to use. Defaults to 1e-3.
 893    as_polygon (bool, optional): if True returns a Polygon, otherwise
 894      returns a one Polygon GeoSeries. Defaults to False.
 895
 896  Returns:
 897    gpd.GeoSeries|geom.Polygon: the resulting union of supplied
 898      polygons.
 899
 900  """
 901  r = RESOLUTION * 10
 902  union = gpd.GeoSeries(
 903    [p.buffer(r, join_style = "mitre", cap_style = "square") for p in gs],
 904    crs = gs.crs).union_all().buffer(-r, join_style = "mitre",
 905                                         cap_style = "square")
 906  if as_polygon:
 907    return union
 908  return gpd.GeoSeries([union], crs = gs.crs)
 909
 910
 911def get_translation_transform(dx:float, dy:float) -> tuple[float,...]:
 912  """Return shapely affine transform tuple for a translation.
 913
 914  Args:
 915      dx (float): translation distance in x direction.
 916      dy (float): translation distance in y direction.
 917
 918  Returns:
 919    list[float]: a six item list of floats, per the shapely.affinity.
 920    affine_transform method, see
 921      https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
 922
 923  """
 924  return (1, 0, 0, 1, dx, dy)
 925
 926
 927def get_rotation_transform(
 928    angle:float, centre:tuple[float,...]|None = None) -> tuple[float,...]:
 929  """Return the shapely affine transform tuple for a rotation.
 930
 931  Rotation is optionally about a supplied centre point.
 932
 933  Args:
 934    angle (float): the angle of rotation (in degrees).
 935    centre (tuple[float], optional): An option centre location. Defaults to
 936    None, which will in turn be converted to (0, 0).
 937
 938  Returns:
 939    list[float]: a six item list of floats, per the shapely.affinity.
 940      affine_transform method, see
 941        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
 942
 943  """
 944  if centre is None or np.allclose((0, 0), centre, atol = RESOLUTION):
 945    a = np.radians(angle)
 946    return (np.cos(a), -np.sin(a), np.sin(a), np.cos(a), 0, 0)
 947  dx, dy = centre
 948  t1 = get_translation_transform(-dx, -dy)
 949  r = get_rotation_transform(angle)
 950  t2 = get_translation_transform(dx, dy)
 951  return combine_transforms([t1, r, t2])
 952
 953
 954def get_reflection_transform(
 955    angle:float, centre:tuple[float,...]|None = None) -> tuple[float,...]:
 956  """Return a shapely affine reflection transform.
 957
 958  Reflection is in a line at the specified angle, optionally through a
 959  specified centre point.
 960
 961  Args:
 962    angle (float): angle to the x-axis of the line of reflection.
 963    centre (tuple[float], optional): point through which the line of
 964      reflection passes. Defaults to None, which
 965      will in turn be converted to (0, 0).
 966
 967  Returns:
 968    list[float]: a six item list of floats, per the shapely.affinity.
 969      affine_transform method, see
 970        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
 971
 972  """
 973  if centre is None or np.allclose((0, 0), centre, atol = RESOLUTION):
 974    A = 2 * np.radians(angle)
 975    return (np.cos(A), np.sin(A), np.sin(A), -np.cos(A), 0, 0)
 976  dx, dy = centre
 977  t1 = get_translation_transform(-dx, -dy)
 978  r = get_reflection_transform(angle)
 979  t2 = get_translation_transform(dx, dy)
 980  return combine_transforms([t1, r, t2])
 981
 982
 983def combine_transforms(transforms:list[tuple[float,...]]) -> tuple[float,...]:
 984  """Return shapely affine transform that combines supplied list of transforms.
 985
 986  Transforms are applied in the order of the supplied list
 987
 988  Args:
 989    transforms (list[list[float]]): sequence of transforms to combine.
 990
 991  Returns:
 992    list[float]: a transform tuple combining the supplied transforms applied
 993      in order, see
 994        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
 995
 996  """
 997  result = np.identity(3)
 998  for t in transforms:
 999    result = as_numpy_matrix(t) @ result
1000  return as_shapely_transform(result)
1001
1002
1003def reverse_transform(transform:tuple[float,...]) -> tuple[float,...]:
1004  """Return inverse shapely affine transform of the supplied transform.
1005
1006  Args:
1007    transform (list[float]): the transform for which the inverse is desired.
1008
1009  Returns:
1010    list[float]: shapely affine transform tuple that will invert the supplied
1011      transform.
1012
1013  """
1014  return as_shapely_transform(
1015    np.linalg.inv(as_numpy_matrix(transform)))
1016
1017
1018def as_shapely_transform(arr:np.ndarray) -> tuple[float,...]:
1019  """Return shapely affine transform equivalent to supplied numpy matrix.
1020
1021  The numpy matrix is a conventional augmented (i.e. 3x3) affine transform
1022  matrix.
1023
1024  Args:
1025    arr (np.array): augmented affine transform matrix of the desired
1026      transform.
1027
1028  Returns:
1029    list[float]: desired shapely affine transform list of floats.
1030
1031  """
1032  return (arr[0][0], arr[0][1], arr[1][0], arr[1][1], arr[0][2], arr[1][2])
1033
1034
1035def as_numpy_matrix(transform:tuple[float,...]) -> np.ndarray:
1036  """Convert shapely affine transform to augmented affine transform matrix.
1037
1038  Args:
1039    transform (list[float]): the transform in shapely format.
1040
1041  Returns:
1042    np.array: the transform in numpy matrix format.
1043
1044  """
1045  return np.array([[transform[0], transform[1], transform[4]],
1046                   [transform[2], transform[3], transform[5]],
1047                   [           0,            0,            1]])
1048
1049
1050class StraightLine(NamedTuple):
1051  """Straight line in canonical Ax+By+C=0 form."""
1052
1053  A: float
1054  B: float
1055  C: float
1056
1057
1058def get_straight_line(
1059    p1:geom.Point, p2:geom.Point,
1060    perpendicular:bool = False,
1061  ) -> StraightLine:
1062  """Return StraightLine p1-p2, optionally perpendicular through the midpoint.
1063
1064  Args:
1065    p1 (geom.Point): First point.
1066    p2 (geom.Point): Second point
1067    perpendicular (bool, optional): If True the perpendicular bisector is
1068      returned. Defaults to False.
1069
1070  Returns:
1071    StraightLine: _description_
1072
1073  """
1074  if perpendicular:
1075    ls = affine.rotate(geom.LineString([p1, p2]), 90)
1076    pts = list(ls.coords)
1077    x1, y1 = pts[0]
1078    x2, y2 = pts[1]
1079  else:
1080    x1, y1 = p1.x, p1.y
1081    x2, y2 = p2.x, p2.y
1082  return StraightLine(y1 - y2, x2 - x1, x1 * y2 - x2 * y1)
1083
1084
1085def get_intersection(line1:StraightLine,
1086                     line2:StraightLine) -> geom.Point|None:
1087  """Return point of intersection of straight lines.
1088
1089  Args:
1090    line1 (StraightLine): First straight line.
1091    line2 (StraightLine): Second straight line.
1092
1093  Returns:
1094    geom.Point|None: Point is returned if the lines intersect, otherwise None.
1095
1096  """
1097  x_set, y_set = False, False
1098  denominator = line1.A * line2.B - line2.A * line1.B
1099  if np.isclose(line1.A, 0, atol = 1e-4, rtol = 1e-4):
1100    y = -line1.C / line1.B
1101    y_set = True
1102  elif np.isclose(line2.A, 0, atol = 1e-4, rtol = 1e-4):
1103    y = -line2.C / line2.B
1104    y_set = True
1105  if np.isclose(line1.B, 0, atol = 1e-4, rtol = 1e-4):
1106    x = -line1.C / line1.A
1107    x_set = True
1108  elif np.isclose(line2.B, 0, atol = 1e-4, rtol = 1e-4):
1109    x = -line2.C / line2.A
1110    x_set = True
1111  if np.isclose(denominator, 0, atol = 1e-4, rtol = 1e-4):
1112    return None
1113  x = x if x_set else (line1.B * line2.C - line2.B * line1.C) / denominator
1114  y = y if y_set else (line1.C * line2.A - line2.C * line1.A) / denominator
1115  return geom.Point(x, y)