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    mat = np.identity(n) + np.roll(np.identity(n), 1, axis = 1)
 373    fractions = np.linalg.inv(mat) @ side_lengths / side_lengths
 374    if not(np.isclose(np.mean(fractions), 
 375                      0.5, rtol= RESOLUTION, atol = RESOLUTION)):
 376      return False
 377    ones = [np.isclose(f, 1, rtol = RESOLUTION, atol = RESOLUTION) 
 378            for f in fractions]
 379    zeros = [np.isclose(f, 0, rtol = RESOLUTION, atol = RESOLUTION) 
 380             for f in fractions]
 381    negatives = [f < 0 for f in fractions]
 382    return not (any(ones) or any(zeros) or any(negatives))
 383  elif n == 4:
 384    # for quadrilaterals odd and even side lengths are equal
 385    return np.isclose(sum(side_lengths[0::2]), 
 386                      sum(side_lengths[1::2]), rtol = RESOLUTION)
 387  else:
 388    # other even numbers of sides... hard to avoid brute force
 389    bisectors = [get_angle_bisector(shape, i) for i in range(n)]
 390    intersects = [b1.intersection(b2) for b1, b2 in
 391                  zip(bisectors, bisectors[1:] + bisectors[:1])
 392                  if b1.intersects(b2)]
 393    distances = [i1.distance(i2) for i1, i2 in
 394                 zip(intersects, intersects[1:] + intersects[:1])]
 395    return all([d <= 2 * RESOLUTION for d in distances])
 396
 397
 398def incentre(shape:geom.Polygon) -> geom.Point:
 399  """A different polygon centre, which produces better results for some
 400  dual tilings where tiles are not regular polygons... see
 401  https://en.wikipedia.org/wiki/Incenter
 402
 403  This method relies on the polygon being tangential, i.e. there is an
 404  inscribed circle to which all sides of the polygon are tangent. It will
 405  work on all the polygons encountered in the Laves tilings, but is not
 406  guaranteed to work on all polygons.
 407
 408  Given that the polygon is tangential, the radius of the inscribed circle is
 409  the [apothem of the polygon](https://en.wikipedia.org/wiki/Apothem) given
 410  by 2 x Area / Perimeter. We apply a parallel offset of this size to two
 411  sides of the polygon and find their intersection to determine the centre of
 412  the circle, givng the incentre of the polygon.
 413
 414    Args:
 415    shape (geom.Polygon): the polygon.
 416
 417  Returns:
 418    geom.Point: the incentre of the polygon.
 419  """
 420  if is_regular_polygon(shape):
 421    return shape.centroid
 422  shape = ensure_cw(shape)
 423  if is_convex(shape):
 424    if is_tangential(shape):  # find the incentre
 425      corners = get_corners(shape, repeat_first = False)
 426      r = get_apothem_length(shape)
 427      e1 = geom.LineString(corners[:2]).parallel_offset(r, side = "right")
 428      e2 = geom.LineString(corners[1:3]).parallel_offset(r, side = "right")
 429      return e1.intersection(e2)  # TODO: add exception if no intersection
 430  return shape.centroid
 431
 432
 433def get_apothem_length(shape:geom.Polygon) -> float:
 434  return 2 * shape.area / geom.LineString(shape.exterior.coords).length
 435
 436
 437def ensure_cw(shape:geom.Polygon) -> geom.Polygon:
 438  """Returns the polygon with its outer boundary vertices in clockwise order.
 439
 440  It is important to note that shapely.set_precision() imposes clockwise order
 441  on polygons, and since it is used widely throughout theses modules, it makes
 442  sense to impose this order.
 443
 444    Args:
 445    shape (geom.Polygon): the polygon.
 446
 447  Returns:
 448    geom.Polygon: the polygon in clockwise order.
 449  """
 450  if geom.LinearRing(shape.exterior.coords).is_ccw:
 451    return shape.reverse()
 452  else:
 453    return shape
 454
 455
 456def get_angle_bisector(shape:geom.Polygon, v = 0) -> geom.LineString:
 457  """Returns a line which is the angle bisector of the specified corner of the
 458  supplied polygon.
 459
 460  Args:
 461      shape (geom.Polygon): the polygon
 462      v (int, optional): index of the corner whose bisector is required. 
 463      Defaults to 0.
 464
 465  Returns:
 466      geom.LineString: line which bisects the specified corner.
 467  """
 468  pts = get_corners(shape, repeat_first = False)
 469  n = len(pts)
 470  p0 = pts[v % n]
 471  p1 = pts[(v + 1) % n]
 472  p2 = pts[(v - 1) % n]
 473  d01 = p0.distance(p1)
 474  d02 = p0.distance(p2)
 475  if d01 < d02:
 476    p2 = geom.LineString([p0, p2]).interpolate(d01)
 477  else:
 478    p1 = geom.LineString([p0, p1]).interpolate(d02)
 479  m = geom.Point((p1.x + p2.x) / 2, (p1.y + p2.y) / 2)
 480  # we have to scale this to much larger in case of an angle near 180
 481  return affine.scale(geom.LineString([p0, m]), 10, 10, origin = p0)
 482
 483
 484def get_side_bisector(shape:geom.Polygon, i:int = 0) -> geom.LineString:
 485  return affine.scale(affine.rotate(get_sides(shape)[i], 90), 100, 100)
 486
 487
 488def _get_interior_vertices(tiles:gpd.GeoDataFrame) -> gpd.GeoSeries:
 489  """Returns points not on the outer boundary of the supplied set
 490  of tiles.
 491
 492  Args:
 493    shapes (gpd.GeoDataFrame): a set of polygons.
 494
 495  Returns:
 496    gpd.GeoSeries: interior vertices of the set of polygons.
 497  """
 498  tiles = gridify(tiles.geometry)
 499  uu = safe_union(tiles, as_polygon = True)
 500  interior_pts = set()
 501  for tile in tiles:
 502    for corner in tile.exterior.coords:
 503      if uu.contains(geom.Point(corner).buffer(1e-3, cap_style = 3)):
 504        interior_pts.add(corner)
 505  return gridify(gpd.GeoSeries([geom.Point(p) for p in interior_pts]))
 506
 507
 508def gridify(
 509      gs:Union[gpd.GeoSeries, gpd.GeoDataFrame, geom.Polygon]
 510    ) -> Union[gpd.GeoSeries, gpd.GeoDataFrame, geom.Polygon]:
 511  """Returns the supplied GeoSeries rounded to the specified precision.
 512
 513  Args:
 514    gs (gpd.GeoSeries): geometries to gridify.
 515    precision (int, optional): digits of precision. Defaults to 6.
 516
 517  Returns:
 518    gpd.GeoSeries: the rounded geometries.
 519  """
 520  # return gpd.GeoSeries(
 521  #   list(gs.apply(
 522  #     wkt.dumps, rounding_precision = precision).apply(wkt.loads)))
 523  if isinstance(gs, (geom.Polygon, geom.Point, geom.MultiPoint,
 524                     geom.MultiPolygon, geom.LineString)) :
 525    return shapely.set_precision(gs, grid_size = RESOLUTION)
 526  if isinstance(gs, gpd.GeoSeries):
 527    return gpd.GeoSeries(
 528      [shapely.set_precision(g, grid_size = RESOLUTION) for g in list(gs)],
 529      crs = gs.crs)
 530  if isinstance(gs, gpd.GeoDataFrame):
 531    gs.geometry = gridify(gs.geometry)
 532    return gs
 533
 534
 535def get_dual_tile_unit(unit: TileUnit) -> gpd.GeoDataFrame:
 536  """Converts supplied TileUnit to a candidate GeoDataFrame of its dual
 537  TileUnit.
 538
 539  NOTE: this is complicated and not remotely guaranteed to work!
 540  a particular issue is that where to place the vertices of the faces
 541  of the dual with respect to the tiles in the original is ill-defined.
 542  This is because the dual process is topologically not metrically defined,
 543  so that exact vertex locations are ambiguous. Tiling duality is defined in
 544  Section 4.2 of Grunbaum B, Shephard G C, 1987 _Tilings and Patterns_ (W. H.
 545  Freeman and Company, New York)
 546
 547  NOTE: In general, this method will work only if all supplied tiles are
 548  regular polygons. A known exception is if the only non-regular polygons are
 549  triangles.
 550
 551  NOTE: 'clean' polygons are required. If supplied polygons have messy
 552  vertices with multiple points where there is only one proper point, bad
 553  things are likely to happen! Consider using `clean_polygon()` on the
 554  tile geometries.
 555
 556  Because of the above limitations, we only return a GeoDataFrame
 557  for inspection. However some `weavingspace.tile_unit.TileUnit` setup
 558  methods in `weavingspace.tiling_geometries` use this method, where we are
 559  confident the returned dual is valid.
 560
 561  Args:
 562    unit (TileUnit): the tiling for which the dual is required.
 563
 564  Returns:
 565    gpd.GeoDataFrame: GeoDataFrame that could be the tiles attribute for
 566      a TileUnit of the dual tiling.
 567  """
 568  # get a local patch of this Tiling
 569  local_patch = unit.get_local_patch(r = 1, include_0 = True)
 570  # Find the interior points of these tiles - these will be guaranteed
 571  # to have a sequence of surrounding tiles incident on them
 572  interior_pts = _get_interior_vertices(local_patch)
 573  # Compile a list of the polygons incident on the interior points
 574  cycles = []
 575  for pt in interior_pts:
 576    cycles.append(
 577      set([poly_id for poly_id, p in enumerate(local_patch.geometry)
 578           if pt.distance(p) < RESOLUTION * 2]))
 579  # convert the polygon ID sequences to (centroid.x, centroid.y, ID) tuples
 580  dual_faces = []
 581  for cycle in cycles:
 582    ids, pts = [], []
 583    for poly_id in cycle:
 584      ids.append(local_patch.tile_id[poly_id])
 585      poly = local_patch.geometry[poly_id]
 586      pts.append(incentre(poly))
 587    # sort them into CW order so they are well formed
 588    sorted_coords = sort_cw([(pt.x, pt.y, id)
 589                              for pt, id in zip(pts, ids)])
 590    dual_faces.append(
 591      (geom.Polygon([(pt_id[0], pt_id[1]) for pt_id in sorted_coords]),
 592       "".join([pt_id[2] for pt_id in sorted_coords])))
 593  # ensure the resulting face centroids are inside the original tile
 594  # displaced a little to avoid uncertainties at corners/edges
 595  # TODO: Check  the logic of this - it seems like dumb luck that it works...
 596  dual_faces = [(f, id) for f, id in dual_faces
 597                if affine.translate(
 598                  unit.prototile.geometry[0], RESOLUTION * 10, 
 599                  RESOLUTION * 10).contains(f.centroid)]
 600  gdf = gpd.GeoDataFrame(
 601    data = {"tile_id": [f[1] for f in dual_faces]}, crs = unit.crs,
 602    geometry = gridify(gpd.GeoSeries([f[0] for f in dual_faces])))
 603  # ensure no duplicates
 604  gdf = gdf.dissolve(by = "tile_id", as_index = False).explode(
 605    index_parts = False, ignore_index = True)
 606
 607  gdf.tile_id = relabel(gdf.tile_id)
 608  return gdf
 609
 610
 611def relabel(data:Iterable) -> list:
 612  """Returns supplied data reassigned with unique values from
 613  string.ascii_letters.
 614
 615  Args:
 616    data (Iterable): the data to relabel
 617
 618  Returns:
 619    list: the reassigned data
 620  """
 621  new_data = {}
 622  d_count = 0
 623  for d in data:
 624    if not d in new_data:
 625      new_data[d] = string.ascii_letters[d_count]
 626      d_count = d_count + 1
 627  return [new_data[d] for d in data]
 628
 629
 630def sort_cw(pts_ids:list[tuple[float, float, str]]):
 631  """Sorts supplied tuple of x, y, ID into clockwise order relative to their
 632  mean centre.
 633
 634  Args:
 635    pts_ids (list[tuple[float, float, str]]): A tuple of a pair of
 636      floats and a string.
 637
 638  Returns:
 639    list: a list in the same format as supplied sorted into
 640      clockwise order of the point locations.
 641  """
 642  x = [p[0] for p in pts_ids]
 643  y = [p[1] for p in pts_ids]
 644  cx, cy = np.mean(x), np.mean(y)
 645  dx = [_ - cx for _ in x]
 646  dy = [_ - cy for _ in y]
 647  angles = [np.arctan2(dy, dx) for dx, dy in zip(dx, dy)]
 648  d = dict(zip(angles, pts_ids))
 649  return [pt_id for angle, pt_id in reversed(sorted(d.items()))]
 650
 651
 652def order_of_pts_cw_around_centre(pts:list[geom.Point], centre:geom.Point):
 653  """Returns the order of the supplied points clockwise relative to supplied 
 654  centre point, i.e. a list of the indices in clockwise order.
 655
 656  Args:
 657      pts (list[geom.Point]): list of points to order.
 658      centre (geom.Point): centre relative to which CW order is determined.
 659
 660  Returns:
 661      _type_: list of reordered points.
 662  """
 663  dx = [p.x - centre.x for p in pts]
 664  dy = [p.y - centre.y for p in pts]
 665  angles = [np.arctan2(dy, dx) for dx, dy in zip(dx, dy)]
 666  d = dict(zip(angles, range(len(pts))))
 667  return [i for angle, i in reversed(sorted(d.items()))]
 668
 669
 670def write_map_to_layers(gdf:gpd.GeoDataFrame, fname:str = "output.gpkg",
 671                        tile_var:str = "tile_id") -> None:
 672  """Writes supplied GeoDataFrame to a GPKG file with layers based on
 673  the tile_var attribute.
 674
 675  Args:
 676    gdf (gpd.GeoDataFrame): the GeoDataFrame.
 677    fname (str, optional): filename to write.
 678    tile_var (str, optional): the attribute to use to separate
 679      output file into layers. Defaults to "tile_id".
 680  """
 681  grouped = gdf.groupby(tile_var, as_index = False)
 682  for e in pd.Series.unique(gdf[tile_var]):
 683    grouped.get_group(e).to_file(fname, layer = e, driver = "GPKG")
 684
 685
 686def get_collapse_distance(geometry:geom.Polygon) -> float:
 687  """Returns the distance under which the supplied polygon will shrink
 688  to nothing if negatively buffered by that distance.
 689
 690  Performs a binary search between an upper bound based on the radius of
 691  the circle of equal area to the polygon, and 0.
 692
 693  Args:
 694    geometry (geom.Polygon): the polygon.
 695
 696  Returns:
 697    float: its collapse distance.
 698  """
 699  if is_convex(geometry):
 700    return get_apothem_length(geometry)
 701  radius = np.sqrt(geometry.area / np.pi)
 702  lower = 0
 703  upper = radius
 704  delta = 1e12
 705  stop_delta = radius / 1000
 706  while delta > stop_delta:
 707    new_r = (lower + upper) / 2
 708    if geometry.buffer(-new_r).area > 0:
 709      lower = new_r
 710      delta = upper - new_r
 711    else:
 712      upper = new_r
 713      delta = new_r - lower
 714  return new_r
 715
 716
 717def get_largest_polygon(polygons:gpd.GeoSeries) -> gpd.GeoSeries:
 718  """Returns the largest polygon in a GeoSeries as a GeoSeries of one polygon.
 719
 720  Args:
 721    polygons (gpd.GeoSeries): the set of polygons to pick from.
 722
 723  Returns:
 724    gpd.GeoSeries: the largest polygon.
 725  """
 726  actual_polygons = [p for p in polygons.geometry
 727                     if isinstance(p, (geom.Polygon, geom.MultiPolygon))]
 728  areas = [p.area for p in actual_polygons]
 729  max_area = max(areas)
 730  return gpd.GeoSeries([actual_polygons[areas.index(max_area)]])
 731
 732
 733def touch_along_an_edge(p1:geom.Polygon, p2:geom.Polygon) -> bool:
 734  """Tests if two polygons touch along an edge.
 735
 736  Checks that the intersection area of the two polygons buffered by
 737  a small amount is large enough to indicate that they neighbour at more
 738  than a corner.
 739
 740  Args:
 741    p1 (geom.Polygon): First polygon
 742    p2 (geom.Polygon): Second polygon
 743
 744  Returns:
 745    bool: True if they neighbour along an edge
 746  """
 747  return p1.buffer(RESOLUTION, join_style = 2, cap_style = 3) \
 748    .intersection(p2.buffer(RESOLUTION, join_style = 2, cap_style = 3)) \
 749      .area > 16 * RESOLUTION ** 2
 750
 751
 752def get_width_height_left_bottom(gs:gpd.GeoSeries) -> tuple[float]:
 753  """Returns width, height, left and bottom limits of a GeoSeries
 754
 755  Args:
 756    gs (geopandas.GeoSeries): GeoSeries for which limits are required.
 757
 758  Returns:
 759    tuple: four float values of width, height, left and bottom of gs.
 760  """
 761  extent = gs.total_bounds
 762  return (extent[2] - extent[0], extent[3] - extent[1],
 763      extent[0], extent[1])
 764
 765
 766def get_bounding_ellipse(shapes:gpd.GeoSeries, 
 767                         mag:float = 1.0) -> gpd.GeoSeries:
 768  """Returns an ellipse containing the supplied shapes.
 769
 770  The method used is to calculate the size of square that would contain
 771  the shapes, if they had an aspect ratio 1, then stretch the circle in
 772  the x, y directions according to the actual aspect ratio of the shapes.
 773
 774  Args:
 775    shapes (gpd.GeoSeries): the shapes to be contained.
 776    mag (float, optional): optionally increase the size of the returned
 777      ellipse by this scale factor. Defaults to 1.0.
 778
 779  Returns:
 780    gpd.GeoSeries: the set of shapes.
 781  """
 782
 783  w, h, l, b = get_width_height_left_bottom(shapes)
 784
 785  c = geom.Point(l + w / 2, b + h / 2)
 786  r = min(w, h) * np.sqrt(2)
 787  circle = [c.buffer(r)]
 788  return gridify(gpd.GeoSeries(circle, crs = shapes.crs).scale(
 789    w / r * mag / np.sqrt(2), h / r * mag / np.sqrt(2), origin = c))
 790
 791
 792def get_tiling_edges(tiles:gpd.GeoSeries) -> gpd.GeoSeries:
 793  """Returns linestring GeoSeries from supplied polygon GeoSeries.
 794
 795  This is used to allow display of edges of tiles in legend when they are
 796  masked by an ellipse (if we instead clip polygons then the ellipse edge
 797  will also show in the result.)
 798
 799  Args:
 800    shapes (gpd.GeoSeries): Polygons to convert.
 801
 802  Returns:
 803    gpd.GeoSeries: LineStrings from the supplied Polygons.
 804  """
 805  tiles = tiles.explode(ignore_index = True)
 806  edges = [geom.LineString(p.exterior.coords) for p in tiles]
 807  return gpd.GeoSeries(edges, crs = tiles.crs)
 808
 809
 810def get_polygon_sector(polygon:geom.Polygon, start:float = 0.0,
 811                       end:float = 1.0) -> geom.Polygon:
 812  """Returns a sector of the provided Polygon.
 813
 814  The returned sector is a section of the polygon boundary between the
 815  normalized start and end positions, and including the polygon centroid.
 816  Should (probably) only be applied to convex polygons.
 817
 818  Args:
 819    shape (geom.Polygon): the Polygon.
 820    start (float): normalized start position along the boundary. Defaults to
 821      0.
 822    end (float): normalized start position along the boundary. Defaults to
 823      1.
 824
 825  Returns:
 826    geom.Polygon: the requested polygon sector.
 827  """
 828  if start == end:
 829    # must return a null polygon since the calling context
 830    # expects to get something back... which most likely
 831    # is needed to align with other data
 832    return geom.Polygon()
 833  if start * end < 0:
 834    # either side of 0/1 so assume required sector includes 0
 835    e1, e2 = min(start, end), max(start, end)
 836    arc1 = shapely.ops.substring(geom.LineString(polygon.exterior.coords),
 837                                 np.mod(e1, 1), 1, normalized = True)
 838    arc2 = shapely.ops.substring(geom.LineString(polygon.exterior.coords),
 839                                 0, e2, normalized = True)
 840    sector = geom.Polygon([polygon.centroid] +
 841               list(arc1.coords) +
 842               list(arc2.coords)[1:])
 843  else:
 844    arc = shapely.ops.substring(geom.LineString(polygon.exterior.coords),
 845                                start, end, normalized = True)
 846    sector = geom.Polygon([polygon.centroid] + list(arc.coords))
 847  return gridify(sector)
 848
 849
 850def repair_polygon(
 851    polygon:Union[geom.Polygon, geom.MultiPolygon, gpd.GeoSeries],
 852    shrink_then_grow:bool = True) -> Union[geom.Polygon, gpd.GeoSeries]:
 853  """Convenience function to 'repair' a shapely polyon or GeoSeries by applying
 854  a negative buffer then the same positive buffer.
 855
 856  Optionally the buffer may be applied in the opposite order (i.e. grow then
 857  shrink). This operation may also convert a MultiPolygon that has some 'stray'
 858  parts to a Polygon.
 859
 860  This is method is often unofficially recommended (on stackexchange etc.)
 861  even in the shapely docs, to resolve topology issues and extraneous
 862  additional vertices appearing when spatial operations are repeatedly
 863  applied.
 864
 865  Args:
 866    p (Union[geom.Polygon, gpd.GeoSeries]): Polygon or GeoSeries to clean.
 867    res (float, optional): buffer size to use. Defaults to 1e-3.
 868    shrink_then_grow (bool, optional): if True the negative buffer is
 869      applied first, otherwise the buffer operations are applied in
 870      reverse. Defaults to True.
 871
 872  Returns:
 873    Union[geom.Polygon, gpd.GeoSeries]: the cleaned Polygon or GeoSeries.
 874  """
 875  if shrink_then_grow:
 876    return polygon.buffer(
 877      -RESOLUTION * 10, join_style = 2, cap_style = 3).buffer(
 878      RESOLUTION * 10, join_style = 2, cap_style = 3)
 879  else:
 880    return polygon.buffer(
 881      RESOLUTION * 10, join_style = 2, cap_style = 3).buffer(
 882      -RESOLUTION * 10, join_style = 2, cap_style = 3)
 883
 884
 885def safe_union(gs:gpd.GeoSeries,
 886               as_polygon:bool = False) -> Union[gpd.GeoSeries, geom.Polygon]:
 887  """Unions the supplied GeoSeries of Polygons while buffering them to avoid
 888  gaps and odd internal floating edges. Optionally returns a Polygon or a
 889  GeoSeries.
 890
 891  Frequently when unioning polygons that are ostensibly adjacent 'rogue'
 892  internal boundaries remain in the result. We can avoid this by buffering the
 893  polygons before unioning them, then reversing the buffer on the unioned
 894  shape.
 895
 896  Args:
 897    gs (gpd.GeoSeries): the Polygons to union.
 898    res (float, optional): size of the buffer to use. Defaults to 1e-3.
 899    as_polygon (bool, optional): if True returns a Polygon, otherwise
 900      returns a one Polygon GeoSeries. Defaults to False.
 901
 902  Returns:
 903    Union[gpd.GeoSeries, geom.Polygon]: the resulting union of supplied
 904      polygons.
 905  """
 906  union = gs.buffer(RESOLUTION * 10, join_style = 2, cap_style = 3) \
 907    .unary_union.buffer(-RESOLUTION * 10, join_style = 2, cap_style = 3)
 908  if as_polygon:
 909    return gridify(union)
 910  else:
 911    return gridify(gpd.GeoSeries([union], crs = gs.crs))
 912
 913
 914def zigzag_between_points(
 915    p0:geom.Point, p1: geom.Point, n:int, h:float = 1.0,  smoothness: int = 0):
 916
 917  template_steps = n * 2 + 1
 918  r = p0.distance(p1)
 919  
 920  x = np.linspace(0, n * np.pi, template_steps, endpoint = True)
 921  y = [np.sin(x) for x in x]
 922  s = interpolate.InterpolatedUnivariateSpline(x, y, k = 2)
 923
 924  spline_steps = (n + smoothness) * 2 + 1
 925  xs = np.linspace(0, n * np.pi, spline_steps, endpoint = True)
 926  ys = s(xs)
 927  
 928  sfx = 1 / max(x) * r
 929  sfy = h * r / 2
 930  theta = np.arctan2(p1.y - p0.y, p1.x - p0.x)
 931
 932  ls = geom.LineString([geom.Point(x, y) for x, y in zip(xs, ys)])
 933  ls = affine.translate(ls, 0, -(ls.bounds[1] + ls.bounds[3]) / 2)
 934  ls = affine.scale(ls, xfact = sfx, yfact = sfy, origin = (0, 0))
 935  ls = affine.rotate(ls, theta, (0, 0), use_radians = True)
 936  x0, y0 = list(ls.coords)[0]
 937  return affine.translate(ls, p0.x - x0, p0.y - y0)
 938
 939
 940def get_translation_transform(dx:float, dy:float) -> list[float]:
 941  """Returns the shapely affine transform tuple for a translation.
 942
 943  Args:
 944      dx (float): translation distance in x direction.
 945      dy (float): translation distance in y direction.
 946
 947  Returns:
 948    list[float]: a six item list of floats, per the shapely.affinity.
 949    affine_transform method, see 
 950      https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
 951  """
 952  return [1, 0, 0, 1, dx, dy]
 953
 954
 955def get_rotation_transform(
 956    angle:float, centre:tuple[float] = None) -> list[float]:
 957  """Returns the shapely affine transform tuple for a rotation, optionally
 958  about a supplied centre point.
 959
 960  Args:
 961      angle (float): the angle of rotation (in degrees).
 962      centre (tuple[float], optional): An option centre location. Defaults to 
 963      None, which will in turn be converted to (0, 0).
 964
 965  Returns:
 966    list[float]: a six item list of floats, per the shapely.affinity.
 967      affine_transform method, see 
 968        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
 969  """
 970  if centre is None or np.allclose((0, 0), centre, atol = RESOLUTION):
 971    a = np.radians(angle)
 972    return [np.cos(a), -np.sin(a), np.sin(a), np.cos(a), 0, 0]
 973  
 974  dx, dy = centre
 975  t1 = get_translation_transform(-dx, -dy)
 976  r = get_rotation_transform(angle)
 977  t2 = get_translation_transform(dx, dy)
 978  return combine_transforms([t1, r, t2])
 979
 980
 981def get_reflection_transform(
 982    angle:float, centre:tuple[float] = None) -> list[float]:
 983  """Returns a shapely affine transform tuple that will reflect a shape
 984  in a line at the specified angle, optionally through a specified centre
 985  point.
 986
 987  Args:
 988    angle (float): angle to the x-axis of the line of reflection.
 989    centre (tuple[float], optional): point through which the line of 
 990      reflection passes. Defaults to None, which
 991      will in turn be converted to (0, 0).
 992
 993  Returns:
 994    list[float]: a six item list of floats, per the shapely.affinity.
 995      affine_transform method, see 
 996        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
 997  """
 998  if centre is None or np.allclose((0, 0), centre, atol = RESOLUTION):
 999    A = 2 * np.radians(angle)
1000    return (np.cos(A), np.sin(A), np.sin(A), -np.cos(A), 0, 0)
1001  dx, dy = centre
1002  t1 = get_translation_transform(-dx, -dy)
1003  r = get_reflection_transform(angle)
1004  t2 = get_translation_transform(dx, dy)
1005  return combine_transforms([t1, r, t2])
1006
1007
1008def combine_transforms(transforms:list[list[float]]) -> list[float]:
1009  """Returns a shapely affine transform list that combines the listed
1010  sequence of transforms applied in order.
1011
1012  Args:
1013    transforms (list[list[float]]): sequence of transforms to combine.
1014
1015  Returns:
1016    list[float]: a transform tuple combining the supplied transforms applied
1017      in order, see 
1018        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
1019  """
1020  result = np.identity(3)
1021  for t in transforms:
1022    result = as_numpy_matrix(t) @ result
1023  return as_shapely_transform(result)
1024
1025
1026def reverse_transform(transform:list[float]) -> list[float]:
1027  """Returns the inverse shapely affine transform of the supplied transform.
1028
1029  Args:
1030    transform (list[float]): the transform for which the inverse is desired.
1031
1032  Returns:
1033    list[float]: shapely affine transform tuple that will invert the supplied
1034      transform.
1035  """
1036  return as_shapely_transform(
1037    np.linalg.inv(as_numpy_matrix(transform)))
1038
1039
1040def as_shapely_transform(arr:np.array) -> list[float]:
1041  """Returns the shapely affine transform list equivalent to the supplied
1042  numpy matrix of a conventional augmented affine transform matrix.
1043
1044  Args:
1045    arr (np.array): augmented affine transform matrix of the desired 
1046      transform.
1047
1048  Returns:
1049    list[float]: desired shapely affine transform list of floats.
1050  """
1051  return [arr[0][0], arr[0][1], arr[1][0], arr[1][1], arr[0][2], arr[1][2]]
1052
1053
1054def as_numpy_matrix(transform:list[float]) -> np.array:
1055  """Converts the supplied shapely affine transform list to an augmented
1056  affine transform matrix in numpy array form. This makes combining transforms
1057  much easier.
1058
1059  Args:
1060    transform (list[float]): the transform in shapely format.
1061
1062  Returns:
1063    np.array: the transform in numpy matrix format.
1064  """
1065  return np.array([[transform[0], transform[1], transform[4]],
1066                   [transform[2], transform[3], transform[5]],
1067                   [           0,            0,            1]])