Module weavingspace.tiling_utils

Functions

def as_numpy_matrix(transform: list[float]) ‑> 

Expand source code

def as_numpy_matrix(transform:list[float]) -> np.array:
  """Converts the supplied shapely affine transform list to an augmented
  affine transform matrix in numpy array form. This makes combining transforms
  much easier.

  Args:
    transform (list[float]): the transform in shapely format.

  Returns:
    np.array: the transform in numpy matrix format.
  """
  return np.array([[transform[0], transform[1], transform[4]],
                   [transform[2], transform[3], transform[5]],
                   [           0,            0,            1]])

Converts the supplied shapely affine transform list to an augmented affine transform matrix in numpy array form. This makes combining transforms much easier.

Args

transform : list[float]

the transform in shapely format.

Returns

np.array

the transform in numpy matrix format.

def as_shapely_transform(arr: np.array) ‑> list[float]

Expand source code

def as_shapely_transform(arr:np.array) -> list[float]:
  """Returns the shapely affine transform list equivalent to the supplied
  numpy matrix of a conventional augmented affine transform matrix.

  Args:
    arr (np.array): augmented affine transform matrix of the desired 
      transform.

  Returns:
    list[float]: desired shapely affine transform list of floats.
  """
  return [arr[0][0], arr[0][1], arr[1][0], arr[1][1], arr[0][2], arr[1][2]]

Returns the shapely affine transform list equivalent to the supplied numpy matrix of a conventional augmented affine transform matrix.

Args

arr : np.array

augmented affine transform matrix of the desired transform.

Returns

list[float]

desired shapely affine transform list of floats.

def centre_offset(shape: geom.Polygon, target: tuple[float] = (0, 0)) ‑> tuple[float]

Expand source code

def centre_offset(shape: geom.Polygon,
                  target:tuple[float] = (0, 0)) -> tuple[float]:
  """Returns vector required to move centroid of polygon to target.

  Args:
    shape (Polygon): polygon to move.
    target (tuple[float], optional): target to move to.
      Defaults to (0, 0).

  Returns:
    tuple[float]: tuple of x, y movement required.
  """
  shape_c = shape.centroid
  return (target[0] - shape_c.x, target[1] - shape_c.y)

Returns vector required to move centroid of polygon to target.

Args

shape : Polygon

polygon to move.

target : tuple[float], optional

target to move to. Defaults to (0, 0).

Returns

tuple[float]

tuple of x, y movement required.

def combine_transforms(transforms: list[list[float]]) ‑> list[float]

Expand source code

def combine_transforms(transforms:list[list[float]]) -> list[float]:
  """Returns a shapely affine transform list that combines the listed
  sequence of transforms applied in order.

  Args:
    transforms (list[list[float]]): sequence of transforms to combine.

  Returns:
    list[float]: a transform tuple combining the supplied transforms applied
      in order, see 
        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
  """
  result = np.identity(3)
  for t in transforms:
    result = as_numpy_matrix(t) @ result
  return as_shapely_transform(result)

Returns a shapely affine transform list that combines the listed sequence of transforms applied in order.

Args

transforms : list[list[float]]

sequence of transforms to combine.

Returns

list[float]

a transform tuple combining the supplied transforms applied in order, see https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations

def ensure_cw(shape: geom.Polygon) ‑> shapely.geometry.polygon.Polygon

Expand source code

def ensure_cw(shape:geom.Polygon) -> geom.Polygon:
  """Returns the polygon with its outer boundary vertices in clockwise order.

  It is important to note that shapely.set_precision() imposes clockwise order
  on polygons, and since it is used widely throughout theses modules, it makes
  sense to impose this order.

    Args:
    shape (geom.Polygon): the polygon.

  Returns:
    geom.Polygon: the polygon in clockwise order.
  """
  if geom.LinearRing(shape.exterior.coords).is_ccw:
    return shape.reverse()
  else:
    return shape

Returns the polygon with its outer boundary vertices in clockwise order.

It is important to note that shapely.set_precision() imposes clockwise order on polygons, and since it is used widely throughout theses modules, it makes sense to impose this order.

Args: shape (geom.Polygon): the polygon.

Returns

geom.Polygon

the polygon in clockwise order.

def geometry_matches(geom1: geom.Polygon, geom2: geom.Polygon)

Expand source code

def geometry_matches(geom1:geom.Polygon, geom2:geom.Polygon):
  a = geom1.area
  return np.isclose(a, geom1.intersection(geom2).area,
                    rtol = RESOLUTION * 100, atol=RESOLUTION * 100)
  return shapely.equals_exact(
    geom1.normalize(), geom2.normalize(), RESOLUTION * 10)

def get_angle_bisector(shape: geom.Polygon, v=0) ‑> shapely.geometry.linestring.LineString

Expand source code

def get_angle_bisector(shape:geom.Polygon, v = 0) -> geom.LineString:
  """Returns a line which is the angle bisector of the specified corner of the
  supplied polygon.

  Args:
      shape (geom.Polygon): the polygon
      v (int, optional): index of the corner whose bisector is required. 
      Defaults to 0.

  Returns:
      geom.LineString: line which bisects the specified corner.
  """
  pts = get_corners(shape, repeat_first = False)
  n = len(pts)
  p0 = pts[v % n]
  p1 = pts[(v + 1) % n]
  p2 = pts[(v - 1) % n]
  d01 = p0.distance(p1)
  d02 = p0.distance(p2)
  if d01 < d02:
    p2 = geom.LineString([p0, p2]).interpolate(d01)
  else:
    p1 = geom.LineString([p0, p1]).interpolate(d02)
  m = geom.Point((p1.x + p2.x) / 2, (p1.y + p2.y) / 2)
  # we have to scale this to much larger in case of an angle near 180
  return affine.scale(geom.LineString([p0, m]), 10, 10, origin = p0)

Returns a line which is the angle bisector of the specified corner of the supplied polygon.

Args

shape : geom.Polygon

the polygon

v : int, optional

index of the corner whose bisector is required.

Defaults to 0.

Returns

geom.LineString

line which bisects the specified corner.

def get_apothem_length(shape: geom.Polygon) ‑> float

Expand source code

def get_apothem_length(shape:geom.Polygon) -> float:
  return 2 * shape.area / geom.LineString(shape.exterior.coords).length

def get_bounding_ellipse(shapes: gpd.GeoSeries, mag: float = 1.0) ‑> geopandas.geoseries.GeoSeries

Expand source code

def get_bounding_ellipse(shapes:gpd.GeoSeries, 
                         mag:float = 1.0) -> gpd.GeoSeries:
  """Returns an ellipse containing the supplied shapes.

  The method used is to calculate the size of square that would contain
  the shapes, if they had an aspect ratio 1, then stretch the circle in
  the x, y directions according to the actual aspect ratio of the shapes.

  Args:
    shapes (gpd.GeoSeries): the shapes to be contained.
    mag (float, optional): optionally increase the size of the returned
      ellipse by this scale factor. Defaults to 1.0.

  Returns:
    gpd.GeoSeries: the set of shapes.
  """

  w, h, l, b = get_width_height_left_bottom(shapes)

  c = geom.Point(l + w / 2, b + h / 2)
  r = min(w, h) * np.sqrt(2)
  circle = [c.buffer(r)]
  return gridify(gpd.GeoSeries(circle, crs = shapes.crs).scale(
    w / r * mag / np.sqrt(2), h / r * mag / np.sqrt(2), origin = c))

Returns an ellipse containing the supplied shapes.

The method used is to calculate the size of square that would contain the shapes, if they had an aspect ratio 1, then stretch the circle in the x, y directions according to the actual aspect ratio of the shapes.

Args

shapes : gpd.GeoSeries

the shapes to be contained.

mag : float, optional

optionally increase the size of the returned ellipse by this scale factor. Defaults to 1.0.

Returns

gpd.GeoSeries

the set of shapes.

def get_clean_polygon(shape: Union[geom.MultiPolygon, geom.Polygon]) ‑> shapely.geometry.polygon.Polygon

Expand source code

def get_clean_polygon(
    shape:Union[geom.MultiPolygon,geom.Polygon]) -> geom.Polygon:
  """Returns polygon with any successive corners that are co-linear along a 
  side or very close to one another removed. 
  
  Particularly useful for tidying polygons weave tilings that have been 
  assembled from multiple 'cells' in the weave grid.

  Args:
    shape (Union[geom.MultiPolygon,geom.Polygon]): polygon to clean.

  Returns:
    geom.Polygon: cleaned polygon.
  """
  if isinstance(shape, geom.MultiPolygon):
    return geom.MultiPolygon([get_clean_polygon(p) for p in shape.geoms])
  else:
    corners = get_corners(shape, repeat_first = False)
    # first remove any 'near neighbour' corners
    distances = get_side_lengths(shape)
    to_remove = [np.isclose(d, 0, rtol = RESOLUTION, atol = 10 * RESOLUTION) 
                 for d in distances]
    corners = [c for c, r in zip(corners, to_remove) if not r]
    # next remove any that are colinear
    p = geom.Polygon(corners)
    corners = get_corners(p, repeat_first = False)
    angles = get_interior_angles(p)
    to_remove = [np.isclose(a, 0, rtol = RESOLUTION, atol = RESOLUTION) or
                 np.isclose(a, 180, rtol = RESOLUTION, atol = RESOLUTION)
                 for a in angles]
    corners = [c for c, r in zip(corners, to_remove) if not r]
    return gridify(geom.Polygon(corners))

Returns polygon with any successive corners that are co-linear along a side or very close to one another removed.

Particularly useful for tidying polygons weave tilings that have been assembled from multiple 'cells' in the weave grid.

Args

shape : Union[geom.MultiPolygon,geom.Polygon]

polygon to clean.

Returns

geom.Polygon

cleaned polygon.

def get_collapse_distance(geometry: geom.Polygon) ‑> float

Expand source code

def get_collapse_distance(geometry:geom.Polygon) -> float:
  """Returns the distance under which the supplied polygon will shrink
  to nothing if negatively buffered by that distance.

  Performs a binary search between an upper bound based on the radius of
  the circle of equal area to the polygon, and 0.

  Args:
    geometry (geom.Polygon): the polygon.

  Returns:
    float: its collapse distance.
  """
  if is_convex(geometry):
    return get_apothem_length(geometry)
  radius = np.sqrt(geometry.area / np.pi)
  lower = 0
  upper = radius
  delta = 1e12
  stop_delta = radius / 1000
  while delta > stop_delta:
    new_r = (lower + upper) / 2
    if geometry.buffer(-new_r).area > 0:
      lower = new_r
      delta = upper - new_r
    else:
      upper = new_r
      delta = new_r - lower
  return new_r

Returns the distance under which the supplied polygon will shrink to nothing if negatively buffered by that distance.

Performs a binary search between an upper bound based on the radius of the circle of equal area to the polygon, and 0.

Args

geometry : geom.Polygon

the polygon.

Returns

float

its collapse distance.

def get_corners(shape: geom.Polygon, repeat_first: bool = True) ‑> list[shapely.geometry.point.Point]

Expand source code

def get_corners(shape:geom.Polygon,
                repeat_first:bool = True) -> list[geom.Point]:
  """Returns a list of geom.Points around the boundary of a polygon, optionally
  repeating the first. Does no simplification (e.g. if a line segment has a 
  'corner' along its length, it is NOT removed; see get_clean_polygon for 
  that). Points have precision set to the package default tiling_utils.
  RESOLUTION.

  Args:
    shape (geom.Polygon): polygon whose corners are required.
    repeat_first (bool, optional): if True the first corner is repeated in the 
      returned list, if False it is omitted. Defaults to True.

  Returns:
    list[geom.Point]: list of geom.Point vertices of the polygon.
  """
  corners = [gridify(geom.Point(pt)) for pt in shape.exterior.coords]
  if repeat_first:
    return corners
  else:
    return corners[:-1]

Returns a list of geom.Points around the boundary of a polygon, optionally repeating the first. Does no simplification (e.g. if a line segment has a 'corner' along its length, it is NOT removed; see get_clean_polygon for that). Points have precision set to the package default tiling_utils. RESOLUTION.

Args

shape : geom.Polygon

polygon whose corners are required.

repeat_first : bool, optional

if True the first corner is repeated in the returned list, if False it is omitted. Defaults to True.

Returns

list[geom.Point]

list of geom.Point vertices of the polygon.

def get_dual_tile_unit(unit: TileUnit) ‑> gpd.GeoDataFrame

Expand source code

def get_dual_tile_unit(unit: TileUnit) -> gpd.GeoDataFrame:
  """Converts supplied TileUnit to a candidate GeoDataFrame of its dual
  TileUnit.

  NOTE: this is complicated and not remotely guaranteed to work!
  a particular issue is that where to place the vertices of the faces
  of the dual with respect to the tiles in the original is ill-defined.
  This is because the dual process is topologically not metrically defined,
  so that exact vertex locations are ambiguous. Tiling duality is defined in
  Section 4.2 of Grunbaum B, Shephard G C, 1987 _Tilings and Patterns_ (W. H.
  Freeman and Company, New York)

  NOTE: In general, this method will work only if all supplied tiles are
  regular polygons. A known exception is if the only non-regular polygons are
  triangles.

  NOTE: 'clean' polygons are required. If supplied polygons have messy
  vertices with multiple points where there is only one proper point, bad
  things are likely to happen! Consider using `clean_polygon()` on the
  tile geometries.

  Because of the above limitations, we only return a GeoDataFrame
  for inspection. However some `weavingspace.tile_unit.TileUnit` setup
  methods in `weavingspace.tiling_geometries` use this method, where we are
  confident the returned dual is valid.

  Args:
    unit (TileUnit): the tiling for which the dual is required.

  Returns:
    gpd.GeoDataFrame: GeoDataFrame that could be the tiles attribute for
      a TileUnit of the dual tiling.
  """
  # get a local patch of this Tiling
  local_patch = unit.get_local_patch(r = 1, include_0 = True)
  # Find the interior points of these tiles - these will be guaranteed
  # to have a sequence of surrounding tiles incident on them
  interior_pts = _get_interior_vertices(local_patch)
  # Compile a list of the polygons incident on the interior points
  cycles = []
  for pt in interior_pts:
    cycles.append(
      set([poly_id for poly_id, p in enumerate(local_patch.geometry)
           if pt.distance(p) < RESOLUTION * 2]))
  # convert the polygon ID sequences to (centroid.x, centroid.y, ID) tuples
  dual_faces = []
  for cycle in cycles:
    ids, pts = [], []
    for poly_id in cycle:
      ids.append(local_patch.tile_id[poly_id])
      poly = local_patch.geometry[poly_id]
      pts.append(incentre(poly))
    # sort them into CW order so they are well formed
    sorted_coords = sort_cw([(pt.x, pt.y, id)
                              for pt, id in zip(pts, ids)])
    dual_faces.append(
      (geom.Polygon([(pt_id[0], pt_id[1]) for pt_id in sorted_coords]),
       "".join([pt_id[2] for pt_id in sorted_coords])))
  # ensure the resulting face centroids are inside the original tile
  # displaced a little to avoid uncertainties at corners/edges
  # TODO: Check  the logic of this - it seems like dumb luck that it works...
  dual_faces = [(f, id) for f, id in dual_faces
                if affine.translate(
                  unit.prototile.geometry[0], RESOLUTION * 10, 
                  RESOLUTION * 10).contains(f.centroid)]
  gdf = gpd.GeoDataFrame(
    data = {"tile_id": [f[1] for f in dual_faces]}, crs = unit.crs,
    geometry = gridify(gpd.GeoSeries([f[0] for f in dual_faces])))
  # ensure no duplicates
  gdf = gdf.dissolve(by = "tile_id", as_index = False).explode(
    index_parts = False, ignore_index = True)

  gdf.tile_id = relabel(gdf.tile_id)
  return gdf

Converts supplied TileUnit to a candidate GeoDataFrame of its dual TileUnit.

NOTE: this is complicated and not remotely guaranteed to work! a particular issue is that where to place the vertices of the faces of the dual with respect to the tiles in the original is ill-defined. This is because the dual process is topologically not metrically defined, so that exact vertex locations are ambiguous. Tiling duality is defined in Section 4.2 of Grunbaum B, Shephard G C, 1987 Tilings and Patterns (W. H. Freeman and Company, New York)

NOTE: In general, this method will work only if all supplied tiles are regular polygons. A known exception is if the only non-regular polygons are triangles.

NOTE: 'clean' polygons are required. If supplied polygons have messy vertices with multiple points where there is only one proper point, bad things are likely to happen! Consider using clean_polygon() on the tile geometries.

Because of the above limitations, we only return a GeoDataFrame for inspection. However some TileUnit setup methods in weavingspace.tiling_geometries use this method, where we are confident the returned dual is valid.

Args

unit : TileUnit

the tiling for which the dual is required.

Returns

gpd.GeoDataFrame

GeoDataFrame that could be the tiles attribute for a TileUnit of the dual tiling.

def get_inner_angle(p1: geom.Point, p2: geom.Point, p3: geom.Point) ‑> float

Expand source code

def get_inner_angle(p1:geom.Point, p2:geom.Point, p3:geom.Point) -> float:
  r"""Returns the angle (in degrees) between line p1-p2 and p2-p3, i.e., the 
  angle A below
  
            p2
           / \ 
          / A \ 
        p1     p3
  
  Args:
    p1 (geom.Point): first point.
    p2 (geom.Point): second 'corner' point.
    p3 (geom.Point): third point.

  Returns:
      float: angle in degrees.
  """
  return np.degrees(np.arctan2(p3.y - p2.y, p3.x - p2.x) -
                    np.arctan2(p1.y - p2.y, p1.x - p2.x)) % 360

Returns the angle (in degrees) between line p1-p2 and p2-p3, i.e., the angle A below

      p2
     / \ 
    / A \ 
  p1     p3

Args

p1 : geom.Point

first point.

p2 : geom.Point

second 'corner' point.

p3 : geom.Point

third point.

Returns

float

angle in degrees.

def get_interior_angles(shape: geom.Polygon) ‑> list[float]

Expand source code

def get_interior_angles(shape:geom.Polygon) -> list[float]:
  """Returns angles (in degrees) between successive edges of a polygon. No 
  polygon simplification is carried out so some angles may be 180 (i.e. a 
  'corner' along a side, such that successive sides are colinear). 

  Args:
    shape (geom.Polygon): polygon whose angles are required.

  Returns:
    list[float]: list of angles.
  """
  corners = get_corners(shape, repeat_first = False)
  wrap_corners = corners[-1:] + corners + corners[:1]
  triples = zip(wrap_corners[:-2], wrap_corners[1:-1], wrap_corners[2:])
  return [get_inner_angle(p1, p2, p3) for p1, p2, p3 in triples]

Returns angles (in degrees) between successive edges of a polygon. No polygon simplification is carried out so some angles may be 180 (i.e. a 'corner' along a side, such that successive sides are colinear).

Args

shape : geom.Polygon

polygon whose angles are required.

Returns

list[float]

list of angles.

def get_largest_polygon(polygons: gpd.GeoSeries) ‑> geopandas.geoseries.GeoSeries

Expand source code

def get_largest_polygon(polygons:gpd.GeoSeries) -> gpd.GeoSeries:
  """Returns the largest polygon in a GeoSeries as a GeoSeries of one polygon.

  Args:
    polygons (gpd.GeoSeries): the set of polygons to pick from.

  Returns:
    gpd.GeoSeries: the largest polygon.
  """
  actual_polygons = [p for p in polygons.geometry
                     if isinstance(p, (geom.Polygon, geom.MultiPolygon))]
  areas = [p.area for p in actual_polygons]
  max_area = max(areas)
  return gpd.GeoSeries([actual_polygons[areas.index(max_area)]])

Returns the largest polygon in a GeoSeries as a GeoSeries of one polygon.

Args

polygons : gpd.GeoSeries

the set of polygons to pick from.

Returns

gpd.GeoSeries

the largest polygon.

def get_outer_angle(p1: geom.Point, p2: geom.Point, p3: geom.Point) ‑> float

Expand source code

def get_outer_angle(p1:geom.Point, p2:geom.Point, p3:geom.Point) -> float:
  r"""Returns outer angle (in degrees) between lines p1-p2 and p2-p3, i.e., the 
  angle A below
                
              /
             /
            p2 A
           / \ 
          /   \ 
        p1     p3
  
  Args:
    p1 (geom.Point): first point.
    p2 (geom.Point): second 'corner' point.
    p3 (geom.Point): third point.

  Returns:
    float: angle in degrees.
  """
  return 180 - get_inner_angle(p1, p2, p3)

Returns outer angle (in degrees) between lines p1-p2 and p2-p3, i.e., the angle A below

        /
       /
      p2 A
     / \ 
    /   \ 
  p1     p3

Args

p1 : geom.Point

first point.

p2 : geom.Point

second 'corner' point.

p3 : geom.Point

third point.

Returns

float

angle in degrees.

def get_polygon_sector(polygon: geom.Polygon, start: float = 0.0, end: float = 1.0) ‑> shapely.geometry.polygon.Polygon

Expand source code

def get_polygon_sector(polygon:geom.Polygon, start:float = 0.0,
                       end:float = 1.0) -> geom.Polygon:
  """Returns a sector of the provided Polygon.

  The returned sector is a section of the polygon boundary between the
  normalized start and end positions, and including the polygon centroid.
  Should (probably) only be applied to convex polygons.

  Args:
    shape (geom.Polygon): the Polygon.
    start (float): normalized start position along the boundary. Defaults to
      0.
    end (float): normalized start position along the boundary. Defaults to
      1.

  Returns:
    geom.Polygon: the requested polygon sector.
  """
  if start == end:
    # must return a null polygon since the calling context
    # expects to get something back... which most likely
    # is needed to align with other data
    return geom.Polygon()
  if start * end < 0:
    # either side of 0/1 so assume required sector includes 0
    e1, e2 = min(start, end), max(start, end)
    arc1 = shapely.ops.substring(geom.LineString(polygon.exterior.coords),
                                 np.mod(e1, 1), 1, normalized = True)
    arc2 = shapely.ops.substring(geom.LineString(polygon.exterior.coords),
                                 0, e2, normalized = True)
    sector = geom.Polygon([polygon.centroid] +
               list(arc1.coords) +
               list(arc2.coords)[1:])
  else:
    arc = shapely.ops.substring(geom.LineString(polygon.exterior.coords),
                                start, end, normalized = True)
    sector = geom.Polygon([polygon.centroid] + list(arc.coords))
  return gridify(sector)

Returns a sector of the provided Polygon.

The returned sector is a section of the polygon boundary between the normalized start and end positions, and including the polygon centroid. Should (probably) only be applied to convex polygons.

Args

shape : geom.Polygon

the Polygon.

start : float

normalized start position along the boundary. Defaults to 0.

end : float

normalized start position along the boundary. Defaults to 1.

Returns

geom.Polygon

the requested polygon sector.

def get_reflection_transform(angle: float, centre: tuple[float] = None) ‑> list[float]

Expand source code

def get_reflection_transform(
    angle:float, centre:tuple[float] = None) -> list[float]:
  """Returns a shapely affine transform tuple that will reflect a shape
  in a line at the specified angle, optionally through a specified centre
  point.

  Args:
    angle (float): angle to the x-axis of the line of reflection.
    centre (tuple[float], optional): point through which the line of 
      reflection passes. Defaults to None, which
      will in turn be converted to (0, 0).

  Returns:
    list[float]: a six item list of floats, per the shapely.affinity.
      affine_transform method, see 
        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
  """
  if centre is None or np.allclose((0, 0), centre, atol = RESOLUTION):
    A = 2 * np.radians(angle)
    return (np.cos(A), np.sin(A), np.sin(A), -np.cos(A), 0, 0)
  dx, dy = centre
  t1 = get_translation_transform(-dx, -dy)
  r = get_reflection_transform(angle)
  t2 = get_translation_transform(dx, dy)
  return combine_transforms([t1, r, t2])

Returns a shapely affine transform tuple that will reflect a shape in a line at the specified angle, optionally through a specified centre point.

Args

angle : float

angle to the x-axis of the line of reflection.

centre : tuple[float], optional

point through which the line of reflection passes. Defaults to None, which will in turn be converted to (0, 0).

Returns

list[float]

a six item list of floats, per the shapely.affinity. affine_transform method, see https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations

def get_regular_polygon(spacing, n: int) ‑> shapely.geometry.polygon.Polygon

Expand source code

def get_regular_polygon(spacing, n:int) -> geom.Polygon:
  """Returns regular polygon with n sides centered on (0, 0) with a horizontal base, and height given by spacing.

  Args:
    spacing (_type_): required height.
    n (int): number of sides.

  Returns:
    geom.Polygon: required geom.Polygon.
  """
  R = spacing / np.cos(np.radians(180 / n)) / 2
  a0 = -90 + 180 / n
  a_diff = 360 / n
  angles = [a0 + a * a_diff for a in range(n)]
  corners = [(R * np.cos(np.radians(a)),
              R * np.sin(np.radians(a))) for a in angles]
  return gridify(geom.Polygon(corners))

Returns regular polygon with n sides centered on (0, 0) with a horizontal base, and height given by spacing.

Args

spacing : _type_

required height.

n : int

number of sides.

Returns

geom.Polygon

required geom.Polygon.

def get_rotation_transform(angle: float, centre: tuple[float] = None) ‑> list[float]

Expand source code

def get_rotation_transform(
    angle:float, centre:tuple[float] = None) -> list[float]:
  """Returns the shapely affine transform tuple for a rotation, optionally
  about a supplied centre point.

  Args:
      angle (float): the angle of rotation (in degrees).
      centre (tuple[float], optional): An option centre location. Defaults to 
      None, which will in turn be converted to (0, 0).

  Returns:
    list[float]: a six item list of floats, per the shapely.affinity.
      affine_transform method, see 
        https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
  """
  if centre is None or np.allclose((0, 0), centre, atol = RESOLUTION):
    a = np.radians(angle)
    return [np.cos(a), -np.sin(a), np.sin(a), np.cos(a), 0, 0]
  
  dx, dy = centre
  t1 = get_translation_transform(-dx, -dy)
  r = get_rotation_transform(angle)
  t2 = get_translation_transform(dx, dy)
  return combine_transforms([t1, r, t2])

Returns the shapely affine transform tuple for a rotation, optionally about a supplied centre point.

Args

angle : float

the angle of rotation (in degrees).

centre : tuple[float], optional

An option centre location. Defaults to

None, which will in turn be converted to (0, 0).

Returns

list[float]

a six item list of floats, per the shapely.affinity. affine_transform method, see https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations

def get_side_bearings(shape: geom.Polygon) ‑> tuple[tuple[float]]

Expand source code

def get_side_bearings(shape:geom.Polygon) -> tuple[tuple[float]]:
  """Returns a list of angles (in degrees) between the sides of a polygon and
  the positive x-axis, when proceeding from the first point in each side to its
  end point. This should usually be CW around the polygon.

  Args:
    shape (geom.Polygon): polygon whose side bearings are required.

  Returns:
    tuple[tuple[float]]: tuple of bearings of each edge.
  """
  return tuple([np.degrees(np.arctan2(
    e.coords[1][1] - e.coords[0][1], e.coords[1][0] - e.coords[0][0])) 
    for e in get_sides(shape)])

Returns a list of angles (in degrees) between the sides of a polygon and the positive x-axis, when proceeding from the first point in each side to its end point. This should usually be CW around the polygon.

Args

shape : geom.Polygon

polygon whose side bearings are required.

Returns

tuple[tuple[float]]

tuple of bearings of each edge.

def get_side_bisector(shape: geom.Polygon, i: int = 0) ‑> shapely.geometry.linestring.LineString

Expand source code

def get_side_bisector(shape:geom.Polygon, i:int = 0) -> geom.LineString:
  return affine.scale(affine.rotate(get_sides(shape)[i], 90), 100, 100)

def get_side_lengths(shape: geom.Polygon) ‑> list[float]

Expand source code

def get_side_lengths(shape:geom.Polygon) -> list[float]:
  """Returns list of lengths of polygon sides. No simplification for corners
  along sides is carried out.

  Args:
    shape (geom.Polygon): polygon whose edge lengths are required.

  Returns:
    list[float]: list of side lengths.
  """
  corners = get_corners(shape)
  return [p1.distance(p2) for p1, p2 in zip(corners[:-1], corners[1:])]

Returns list of lengths of polygon sides. No simplification for corners along sides is carried out.

Args

shape : geom.Polygon

polygon whose edge lengths are required.

Returns

list[float]

list of side lengths.

def get_sides(shape: geom.Polygon) ‑> list[shapely.geometry.linestring.LineString]

Expand source code

def get_sides(shape:geom.Polygon) -> list[geom.LineString]:
  """Returns polygon sides as a list of geom.LineStrings, with resolution set
  to the package default tiling_utils.RESOLUTION. No simplification for 
  successive colinear sides is carried out.

  Args:
    shape (geom.Polygon): polygon whose edges are required.

  Returns:
    list[geom.LineString]: list of geom.LineString sides of the polygon.
  """
  corners = get_corners(shape)
  return [gridify(geom.LineString([p1, p2]))
          for p1, p2 in zip(corners[:-1], corners[1:])]

Returns polygon sides as a list of geom.LineStrings, with resolution set to the package default tiling_utils.RESOLUTION. No simplification for successive colinear sides is carried out.

Args

shape : geom.Polygon

polygon whose edges are required.

Returns

list[geom.LineString]

list of geom.LineString sides of the polygon.

def get_strand_ids(strands_spec: str) ‑> tuple[list[str]]

Expand source code

def get_strand_ids(strands_spec: str) -> tuple[list[str]]:
  """Conversts a strands specification string to a list of lists of strand
  labels.

  Args:
    strands_spec (str): string format "a|bc|(de)f" | separates strands in
      each direction and () designates combining labels into a single
      strand that will be sliced lengthwise. Example output:

        "a|bc|(de)f" -> (["a"], ["b", "c"], ["de", "f"])

    Superflous parentheses are removed, but no other error-checks are
      applied.

  Returns:
    tuple[str]: tuple of lists of labels for each set of strands.
  """
  strand_ids = [_parse_strand_label(s) for s in strands_spec.split("|")]
  strand_ids = strand_ids if len(strand_ids) == 3 else strand_ids + [[""]]
  return tuple(strand_ids)

Conversts a strands specification string to a list of lists of strand labels.

Args

strands_spec : str

string format "a|bc|(de)f" | separates strands in each direction and () designates combining labels into a single strand that will be sliced lengthwise. Example output:

"a|bc|(de)f" -> (["a"], ["b", "c"], ["de", "f"])

Superflous parentheses are removed, but no other error-checks are applied.

Returns

tuple[str]

tuple of lists of labels for each set of strands.

def get_tiling_edges(tiles: gpd.GeoSeries) ‑> geopandas.geoseries.GeoSeries

Expand source code

def get_tiling_edges(tiles:gpd.GeoSeries) -> gpd.GeoSeries:
  """Returns linestring GeoSeries from supplied polygon GeoSeries.

  This is used to allow display of edges of tiles in legend when they are
  masked by an ellipse (if we instead clip polygons then the ellipse edge
  will also show in the result.)

  Args:
    shapes (gpd.GeoSeries): Polygons to convert.

  Returns:
    gpd.GeoSeries: LineStrings from the supplied Polygons.
  """
  tiles = tiles.explode(ignore_index = True)
  edges = [geom.LineString(p.exterior.coords) for p in tiles]
  return gpd.GeoSeries(edges, crs = tiles.crs)

Returns linestring GeoSeries from supplied polygon GeoSeries.

This is used to allow display of edges of tiles in legend when they are masked by an ellipse (if we instead clip polygons then the ellipse edge will also show in the result.)

Args

shapes : gpd.GeoSeries

Polygons to convert.

Returns

gpd.GeoSeries

LineStrings from the supplied Polygons.

def get_translation_transform(dx: float, dy: float) ‑> list[float]

Expand source code

def get_translation_transform(dx:float, dy:float) -> list[float]:
  """Returns the shapely affine transform tuple for a translation.

  Args:
      dx (float): translation distance in x direction.
      dy (float): translation distance in y direction.

  Returns:
    list[float]: a six item list of floats, per the shapely.affinity.
    affine_transform method, see 
      https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations
  """
  return [1, 0, 0, 1, dx, dy]

Returns the shapely affine transform tuple for a translation.

Args

dx : float

translation distance in x direction.

dy : float

translation distance in y direction.

Returns

list[float]

a six item list of floats, per the shapely.affinity.

affine_transform method, see https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations

def get_width_height_left_bottom(gs: gpd.GeoSeries) ‑> tuple[float]

Expand source code

def get_width_height_left_bottom(gs:gpd.GeoSeries) -> tuple[float]:
  """Returns width, height, left and bottom limits of a GeoSeries

  Args:
    gs (geopandas.GeoSeries): GeoSeries for which limits are required.

  Returns:
    tuple: four float values of width, height, left and bottom of gs.
  """
  extent = gs.total_bounds
  return (extent[2] - extent[0], extent[3] - extent[1],
      extent[0], extent[1])

Returns width, height, left and bottom limits of a GeoSeries

Args

gs : geopandas.GeoSeries

GeoSeries for which limits are required.

Returns

tuple

four float values of width, height, left and bottom of gs.

def gridify(gs: Union[gpd.GeoSeries, gpd.GeoDataFrame, geom.Polygon]) ‑> geopandas.geoseries.GeoSeries | geopandas.geodataframe.GeoDataFrame | shapely.geometry.polygon.Polygon

Expand source code

def gridify(
      gs:Union[gpd.GeoSeries, gpd.GeoDataFrame, geom.Polygon]
    ) -> Union[gpd.GeoSeries, gpd.GeoDataFrame, geom.Polygon]:
  """Returns the supplied GeoSeries rounded to the specified precision.

  Args:
    gs (gpd.GeoSeries): geometries to gridify.
    precision (int, optional): digits of precision. Defaults to 6.

  Returns:
    gpd.GeoSeries: the rounded geometries.
  """
  # return gpd.GeoSeries(
  #   list(gs.apply(
  #     wkt.dumps, rounding_precision = precision).apply(wkt.loads)))
  if isinstance(gs, (geom.Polygon, geom.Point, geom.MultiPoint,
                     geom.MultiPolygon, geom.LineString)) :
    return shapely.set_precision(gs, grid_size = RESOLUTION)
  if isinstance(gs, gpd.GeoSeries):
    return gpd.GeoSeries(
      [shapely.set_precision(g, grid_size = RESOLUTION) for g in list(gs)],
      crs = gs.crs)
  if isinstance(gs, gpd.GeoDataFrame):
    gs.geometry = gridify(gs.geometry)
    return gs

Returns the supplied GeoSeries rounded to the specified precision.

Args

gs : gpd.GeoSeries

geometries to gridify.

precision : int, optional

digits of precision. Defaults to 6.

Returns

gpd.GeoSeries

the rounded geometries.

def incentre(shape: geom.Polygon) ‑> shapely.geometry.point.Point

Expand source code

def incentre(shape:geom.Polygon) -> geom.Point:
  """A different polygon centre, which produces better results for some
  dual tilings where tiles are not regular polygons... see
  https://en.wikipedia.org/wiki/Incenter

  This method relies on the polygon being tangential, i.e. there is an
  inscribed circle to which all sides of the polygon are tangent. It will
  work on all the polygons encountered in the Laves tilings, but is not
  guaranteed to work on all polygons.

  Given that the polygon is tangential, the radius of the inscribed circle is
  the [apothem of the polygon](https://en.wikipedia.org/wiki/Apothem) given
  by 2 x Area / Perimeter. We apply a parallel offset of this size to two
  sides of the polygon and find their intersection to determine the centre of
  the circle, givng the incentre of the polygon.

    Args:
    shape (geom.Polygon): the polygon.

  Returns:
    geom.Point: the incentre of the polygon.
  """
  if is_regular_polygon(shape):
    return shape.centroid
  shape = ensure_cw(shape)
  if is_convex(shape):
    if is_tangential(shape):  # find the incentre
      corners = get_corners(shape, repeat_first = False)
      r = get_apothem_length(shape)
      e1 = geom.LineString(corners[:2]).parallel_offset(r, side = "right")
      e2 = geom.LineString(corners[1:3]).parallel_offset(r, side = "right")
      c = e1.intersection(e2)
      # is_tangential is unreliable, and sometimes we get to here and do not
      # get a Point, but a LineString, or even no intersection, so...
      return c if isinstance(c, geom.Point) else shape.centroid
  return shape.centroid

A different polygon centre, which produces better results for some dual tilings where tiles are not regular polygons… see https://en.wikipedia.org/wiki/Incenter

This method relies on the polygon being tangential, i.e. there is an inscribed circle to which all sides of the polygon are tangent. It will work on all the polygons encountered in the Laves tilings, but is not guaranteed to work on all polygons.

Given that the polygon is tangential, the radius of the inscribed circle is the apothem of the polygon given by 2 x Area / Perimeter. We apply a parallel offset of this size to two sides of the polygon and find their intersection to determine the centre of the circle, givng the incentre of the polygon.

Args: shape (geom.Polygon): the polygon.

Returns

geom.Point

the incentre of the polygon.

def is_convex(shape: geom.Polygon) ‑> bool

Expand source code

def is_convex(shape:geom.Polygon) -> bool:
  """Tests for shape convexity. There are better ways to do this, like
  e.g. using something like the get_interior_angles() function, but simply
  checking if the area is close to that of its convex hull works too!

  Args:
    shape (geom.Polygon): polygon to check

  Returns:
    bool: True if the polygon is convex, False otherwise.
  """
  return np.isclose(shape.area, shape.convex_hull.area, rtol = RESOLUTION)

Tests for shape convexity. There are better ways to do this, like e.g. using something like the get_interior_angles() function, but simply checking if the area is close to that of its convex hull works too!

Args

shape : geom.Polygon

polygon to check

Returns

bool

True if the polygon is convex, False otherwise.

def is_regular_polygon(shape: geom.Polygon) ‑> bool

Expand source code

def is_regular_polygon(shape:geom.Polygon) -> bool:
  """Tests if supplied polygon is regular (i.e. equal sides and angles).

  Args:
      shape (geom.Polygon): polygon to test.

  Returns:
      bool: True if polygon is regular, False if not.
  """
  side_lengths = get_side_lengths(shape)
  angles = get_interior_angles(shape)
  return all(np.isclose(side_lengths, side_lengths[0])) \
     and all(np.isclose(angles, angles[0]))

Tests if supplied polygon is regular (i.e. equal sides and angles).

Args

shape : geom.Polygon

polygon to test.

Returns

bool

True if polygon is regular, False if not.

def is_tangential(shape: geom.Polygon) ‑> bool

Expand source code

def is_tangential(shape:geom.Polygon) -> bool:
  """Determines if the supplied polygon is tangential i.e., it can have
  circle inscribed tangential to all its sides. 
  
  Note that this will fail for polygons with successive colinear sides,
  meaning that polygons should be fully simplified...
  """
  if is_regular_polygon(shape):
    return True
  side_lengths = get_side_lengths(shape)
  n = len(side_lengths)
  if n % 2 == 1:
    if n == 3: #triangles are easy!
      return True
    # odd number of sides there is a nice solvable system  of equations
    # see https://math.stackexchange.com/questions/4065370/tangential-polygons-conditions-on-edge-lengths
    #
    # TODO: this is not reliable for reasons I don't fully understand
    # there is a corresponding crude catch exception in incentre... but
    # this needs a more complete investigation
    mat = np.identity(n) + np.roll(np.identity(n), 1, axis = 1)
    fractions = np.linalg.inv(mat) @ side_lengths / side_lengths
    if not(np.isclose(np.mean(fractions), 
                      0.5, rtol= RESOLUTION, atol = RESOLUTION)):
      return False
    ones = [np.isclose(f, 1, rtol = RESOLUTION, atol = RESOLUTION) 
            for f in fractions]
    zeros = [np.isclose(f, 0, rtol = RESOLUTION, atol = RESOLUTION) 
             for f in fractions]
    negatives = [f < 0 for f in fractions]
    return not (any(ones) or any(zeros) or any(negatives))
  elif n == 4:
    # for quadrilaterals odd and even side lengths are equal
    return np.isclose(sum(side_lengths[0::2]), 
                      sum(side_lengths[1::2]), rtol = RESOLUTION)
  else:
    # other even numbers of sides... hard to avoid brute force
    bisectors = [get_angle_bisector(shape, i) for i in range(n)]
    intersects = [b1.intersection(b2) for b1, b2 in
                  zip(bisectors, bisectors[1:] + bisectors[:1])
                  if b1.intersects(b2)]
    distances = [i1.distance(i2) for i1, i2 in
                 zip(intersects, intersects[1:] + intersects[:1])]
    return all([d <= 2 * RESOLUTION for d in distances])

Determines if the supplied polygon is tangential i.e., it can have circle inscribed tangential to all its sides.

Note that this will fail for polygons with successive colinear sides, meaning that polygons should be fully simplified…

def offset_polygon_corners(polygon: geom.Polygon, offset: int) ‑> shapely.geometry.polygon.Polygon

Expand source code

def offset_polygon_corners(polygon:geom.Polygon, 
                          offset:int) -> geom.Polygon:
  """Returns this polygon but with its first corner offset from its
  original position in the coordinate sequence. The returned polygon will
  be identical but stored differently internally.

  Args:
    polygon (geom.Polygon): the polygon to reorder.
    offset (int): the number of corner positions by which to shift the
      sequence.

  Returns:
      geom.Polygon: the reordered polygon.
  """
  corners = get_corners(polygon, repeat_first = False)
  return geom.Polygon([c for c in corners[offset:] + corners[:offset]]) 

Returns this polygon but with its first corner offset from its original position in the coordinate sequence. The returned polygon will be identical but stored differently internally.

Args

polygon : geom.Polygon

the polygon to reorder.

offset : int

the number of corner positions by which to shift the sequence.

Returns

geom.Polygon

the reordered polygon.

def order_of_pts_cw_around_centre(pts: list[geom.Point], centre: geom.Point)

Expand source code

def order_of_pts_cw_around_centre(pts:list[geom.Point], centre:geom.Point):
  """Returns the order of the supplied points clockwise relative to supplied 
  centre point, i.e. a list of the indices in clockwise order.

  Args:
      pts (list[geom.Point]): list of points to order.
      centre (geom.Point): centre relative to which CW order is determined.

  Returns:
      _type_: list of reordered points.
  """
  dx = [p.x - centre.x for p in pts]
  dy = [p.y - centre.y for p in pts]
  angles = [np.arctan2(dy, dx) for dx, dy in zip(dx, dy)]
  d = dict(zip(angles, range(len(pts))))
  return [i for angle, i in reversed(sorted(d.items()))]

Returns the order of the supplied points clockwise relative to supplied centre point, i.e. a list of the indices in clockwise order.

Args

pts : list[geom.Point]

list of points to order.

centre : geom.Point

centre relative to which CW order is determined.

Returns

_type_

list of reordered points.

def relabel(data: Iterable) ‑> list

Expand source code

def relabel(data:Iterable) -> list:
  """Returns supplied data reassigned with unique values from
  string.ascii_letters.

  Args:
    data (Iterable): the data to relabel

  Returns:
    list: the reassigned data
  """
  new_data = {}
  d_count = 0
  for d in data:
    if not d in new_data:
      new_data[d] = string.ascii_letters[d_count]
      d_count = d_count + 1
  return [new_data[d] for d in data]

Returns supplied data reassigned with unique values from string.ascii_letters.

Args

data : Iterable

the data to relabel

Returns

list

the reassigned data

def repair_polygon(polygon: Union[geom.Polygon, geom.MultiPolygon, gpd.GeoSeries],
shrink_then_grow: bool = True) ‑> shapely.geometry.polygon.Polygon | geopandas.geoseries.GeoSeries

Expand source code

def repair_polygon(
    polygon:Union[geom.Polygon, geom.MultiPolygon, gpd.GeoSeries],
    shrink_then_grow:bool = True) -> Union[geom.Polygon, gpd.GeoSeries]:
  """Convenience function to 'repair' a shapely polyon or GeoSeries by applying
  a negative buffer then the same positive buffer.

  Optionally the buffer may be applied in the opposite order (i.e. grow then
  shrink). This operation may also convert a MultiPolygon that has some 'stray'
  parts to a Polygon.

  This is method is often unofficially recommended (on stackexchange etc.)
  even in the shapely docs, to resolve topology issues and extraneous
  additional vertices appearing when spatial operations are repeatedly
  applied.

  Args:
    p (Union[geom.Polygon, gpd.GeoSeries]): Polygon or GeoSeries to clean.
    res (float, optional): buffer size to use. Defaults to 1e-3.
    shrink_then_grow (bool, optional): if True the negative buffer is
      applied first, otherwise the buffer operations are applied in
      reverse. Defaults to True.

  Returns:
    Union[geom.Polygon, gpd.GeoSeries]: the cleaned Polygon or GeoSeries.
  """
  if shrink_then_grow:
    return polygon.buffer(
      -RESOLUTION * 10, join_style = 2, cap_style = 3).buffer(
      RESOLUTION * 10, join_style = 2, cap_style = 3)
  else:
    return polygon.buffer(
      RESOLUTION * 10, join_style = 2, cap_style = 3).buffer(
      -RESOLUTION * 10, join_style = 2, cap_style = 3)

Convenience function to 'repair' a shapely polyon or GeoSeries by applying a negative buffer then the same positive buffer.

Optionally the buffer may be applied in the opposite order (i.e. grow then shrink). This operation may also convert a MultiPolygon that has some 'stray' parts to a Polygon.

This is method is often unofficially recommended (on stackexchange etc.) even in the shapely docs, to resolve topology issues and extraneous additional vertices appearing when spatial operations are repeatedly applied.

Args

p : Union[geom.Polygon, gpd.GeoSeries]

Polygon or GeoSeries to clean.

res : float, optional

buffer size to use. Defaults to 1e-3.

shrink_then_grow : bool, optional

if True the negative buffer is applied first, otherwise the buffer operations are applied in reverse. Defaults to True.

Returns

Union[geom.Polygon, gpd.GeoSeries]

the cleaned Polygon or GeoSeries.

def reverse_transform(transform: list[float]) ‑> list[float]

Expand source code

def reverse_transform(transform:list[float]) -> list[float]:
  """Returns the inverse shapely affine transform of the supplied transform.

  Args:
    transform (list[float]): the transform for which the inverse is desired.

  Returns:
    list[float]: shapely affine transform tuple that will invert the supplied
      transform.
  """
  return as_shapely_transform(
    np.linalg.inv(as_numpy_matrix(transform)))

Returns the inverse shapely affine transform of the supplied transform.

Args

transform : list[float]

the transform for which the inverse is desired.

Returns

list[float]

shapely affine transform tuple that will invert the supplied transform.

def rotate_preserving_order(polygon: geom.Polygon, angle: float, centre: geom.Point) ‑> shapely.geometry.polygon.Polygon

Expand source code

def rotate_preserving_order(polygon:geom.Polygon, angle:float,
                            centre:geom.Point) -> geom.Polygon:
  """Returns the supplied polygon rotated with the order of its corner points
  preserved (not guaranteed by shapely.affinity.rotate).

  Args:
      polygon (geom.Polygon): polygon to rotate.
      angle (float): desired angle of rotation (in degrees).
      centre (geom.Point): the rotation centre (passed on to shapely.affinity.
        rotate).

  Returns:
      geom.Polygon: rotated polygon.
  """
  corners = get_corners(polygon, repeat_first = False)
  return geom.Polygon([affine.rotate(c, angle, centre) for c in corners])

Returns the supplied polygon rotated with the order of its corner points preserved (not guaranteed by shapely.affinity.rotate).

Args

polygon : geom.Polygon

polygon to rotate.

angle : float

desired angle of rotation (in degrees).

centre : geom.Point

the rotation centre (passed on to shapely.affinity. rotate).

Returns

geom.Polygon

rotated polygon.

def safe_union(gs: gpd.GeoSeries, as_polygon: bool = False) ‑> shapely.geometry.polygon.Polygon | geopandas.geoseries.GeoSeries

Expand source code

def safe_union(gs:gpd.GeoSeries,
               as_polygon:bool = False) -> Union[gpd.GeoSeries, geom.Polygon]:
  """Unions the supplied GeoSeries of Polygons while buffering them to avoid
  gaps and odd internal floating edges. Optionally returns a Polygon or a
  GeoSeries.

  Frequently when unioning polygons that are ostensibly adjacent 'rogue'
  internal boundaries remain in the result. We can avoid this by buffering the
  polygons before unioning them, then reversing the buffer on the unioned
  shape.

  Args:
    gs (gpd.GeoSeries): the Polygons to union.
    res (float, optional): size of the buffer to use. Defaults to 1e-3.
    as_polygon (bool, optional): if True returns a Polygon, otherwise
      returns a one Polygon GeoSeries. Defaults to False.

  Returns:
    Union[gpd.GeoSeries, geom.Polygon]: the resulting union of supplied
      polygons.
  """
  union = gs.buffer(RESOLUTION * 10, join_style = 2, cap_style = 3) \
    .unary_union.buffer(-RESOLUTION * 10, join_style = 2, cap_style = 3)
  if as_polygon:
    return gridify(union)
  else:
    return gridify(gpd.GeoSeries([union], crs = gs.crs))

Unions the supplied GeoSeries of Polygons while buffering them to avoid gaps and odd internal floating edges. Optionally returns a Polygon or a GeoSeries.

Frequently when unioning polygons that are ostensibly adjacent 'rogue' internal boundaries remain in the result. We can avoid this by buffering the polygons before unioning them, then reversing the buffer on the unioned shape.

Args

gs : gpd.GeoSeries

the Polygons to union.

res : float, optional

size of the buffer to use. Defaults to 1e-3.

as_polygon : bool, optional

if True returns a Polygon, otherwise returns a one Polygon GeoSeries. Defaults to False.

Returns

Union[gpd.GeoSeries, geom.Polygon]

the resulting union of supplied polygons.

def sort_cw(pts_ids: list[tuple[float, float, str]])

Expand source code

def sort_cw(pts_ids:list[tuple[float, float, str]]):
  """Sorts supplied tuple of x, y, ID into clockwise order relative to their
  mean centre.

  Args:
    pts_ids (list[tuple[float, float, str]]): A tuple of a pair of
      floats and a string.

  Returns:
    list: a list in the same format as supplied sorted into
      clockwise order of the point locations.
  """
  x = [p[0] for p in pts_ids]
  y = [p[1] for p in pts_ids]
  cx, cy = np.mean(x), np.mean(y)
  dx = [_ - cx for _ in x]
  dy = [_ - cy for _ in y]
  angles = [np.arctan2(dy, dx) for dx, dy in zip(dx, dy)]
  d = dict(zip(angles, pts_ids))
  return [pt_id for angle, pt_id in reversed(sorted(d.items()))]

Sorts supplied tuple of x, y, ID into clockwise order relative to their mean centre.

Args

pts_ids : list[tuple[float, float, str]]

A tuple of a pair of floats and a string.

Returns

list

a list in the same format as supplied sorted into clockwise order of the point locations.

def touch_along_an_edge(p1: geom.Polygon, p2: geom.Polygon) ‑> bool

Expand source code

def touch_along_an_edge(p1:geom.Polygon, p2:geom.Polygon) -> bool:
  """Tests if two polygons touch along an edge.

  Checks that the intersection area of the two polygons buffered by
  a small amount is large enough to indicate that they neighbour at more
  than a corner.

  Args:
    p1 (geom.Polygon): First polygon
    p2 (geom.Polygon): Second polygon

  Returns:
    bool: True if they neighbour along an edge
  """
  return p1.buffer(RESOLUTION, join_style = 2, cap_style = 3) \
    .intersection(p2.buffer(RESOLUTION, join_style = 2, cap_style = 3)) \
      .area > 16 * RESOLUTION ** 2

Tests if two polygons touch along an edge.

Checks that the intersection area of the two polygons buffered by a small amount is large enough to indicate that they neighbour at more than a corner.

Args

p1 : geom.Polygon

First polygon

p2 : geom.Polygon

Second polygon

Returns

bool

True if they neighbour along an edge

def write_map_to_layers(gdf: gpd.GeoDataFrame, fname: str = 'output.gpkg', tile_var: str = 'tile_id') ‑> None

Expand source code

def write_map_to_layers(gdf:gpd.GeoDataFrame, fname:str = "output.gpkg",
                        tile_var:str = "tile_id") -> None:
  """Writes supplied GeoDataFrame to a GPKG file with layers based on
  the tile_var attribute.

  Args:
    gdf (gpd.GeoDataFrame): the GeoDataFrame.
    fname (str, optional): filename to write.
    tile_var (str, optional): the attribute to use to separate
      output file into layers. Defaults to "tile_id".
  """
  grouped = gdf.groupby(tile_var, as_index = False)
  for e in pd.Series.unique(gdf[tile_var]):
    grouped.get_group(e).to_file(fname, layer = e, driver = "GPKG")

Writes supplied GeoDataFrame to a GPKG file with layers based on the tile_var attribute.

Args

gdf : gpd.GeoDataFrame

the GeoDataFrame.

fname : str, optional

filename to write.

tile_var : str, optional

the attribute to use to separate output file into layers. Defaults to "tile_id".

def zigzag_between_points(p0: geom.Point, p1: geom.Point, n: int, h: float = 1.0, smoothness: int = 0)

Expand source code

def zigzag_between_points(
    p0:geom.Point, p1: geom.Point, n:int, h:float = 1.0,  smoothness: int = 0):

  template_steps = n * 2 + 1
  r = p0.distance(p1)
  
  x = np.linspace(0, n * np.pi, template_steps, endpoint = True)
  y = [np.sin(x) for x in x]
  s = interpolate.InterpolatedUnivariateSpline(x, y, k = 2)

  spline_steps = (n + smoothness) * 2 + 1
  xs = np.linspace(0, n * np.pi, spline_steps, endpoint = True)
  ys = s(xs)
  
  sfx = 1 / max(x) * r
  sfy = h * r / 2
  theta = np.arctan2(p1.y - p0.y, p1.x - p0.x)

  ls = geom.LineString([geom.Point(x, y) for x, y in zip(xs, ys)])
  ls = affine.translate(ls, 0, -(ls.bounds[1] + ls.bounds[3]) / 2)
  ls = affine.scale(ls, xfact = sfx, yfact = sfy, origin = (0, 0))
  ls = affine.rotate(ls, theta, (0, 0), use_radians = True)
  x0, y0 = list(ls.coords)[0]
  return affine.translate(ls, p0.x - x0, p0.y - y0)