weavingspace.tiling_utils

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