weavingspace.tileable

  1#!/usr/bin/env python
  2# coding: utf-8
  3
  4"""Implements `weavingspace.tileable.TileShape` and
  5`weavingspace.tileable.Tileable` the base classes for
  6`weavingspace.tile_unit.TileUnit` and `weavingspace.weave_unit.WeaveUnit`.
  7
  8`Tileable` should not be called directly, but is instead accessed from the
  9`weavingspace.tile_unit.TileUnit` or `weavingspace.weave_unit.WeaveUnit`
 10constructor.
 11
 12Several methods of `weavingspace.tileable.Tileable` are generally useful and
 13can be accessed through its subclasses.
 14"""
 15
 16from enum import Enum
 17from typing import Union
 18from dataclasses import dataclass
 19import copy
 20
 21import matplotlib.pyplot as pyplot
 22
 23import geopandas as gpd
 24import numpy as np
 25import shapely.geometry as geom
 26import shapely.affinity as affine
 27
 28import weavingspace.tiling_utils as tiling_utils
 29
 30
 31class TileShape(Enum):
 32  """The available base tile shapes.
 33
 34  NOTE: the TRIANGLE type does not persist, but should be converted to a
 35  DIAMOND or HEXAGON type during `Tileable` construction.
 36  """
 37
 38  RECTANGLE = "rectangle"
 39  HEXAGON = "hexagon"
 40  TRIANGLE = "triangle"
 41  DIAMOND = "diamond"
 42
 43
 44@dataclass
 45class Tileable:
 46  """Class to represent a tileable set of tile geometries.
 47  """
 48
 49  tiles: gpd.GeoDataFrame = None
 50  """the geometries with associated `title_id` attribute encoding their 
 51  different colouring."""
 52  prototile: gpd.GeoDataFrame = None
 53  """the tileable polygon (rectangle, hexagon or diamond)"""
 54  spacing: float = 1000.0
 55  """the tile spacing effectively the resolution of the tiling. Defaults to
 56  1000"""
 57  base_shape: TileShape = TileShape.RECTANGLE
 58  """the tile shape. Defaults to 'RECTANGLE'"""
 59  vectors: dict[tuple[int], tuple[float]] = None
 60  """translation vector symmetries of the tiling"""
 61  regularised_prototile: gpd.GeoDataFrame = None
 62  """polygon containing the tiles of this tileable, usually a union of its
 63  tile polygons"""
 64  crs: int = 3857
 65  """coordinate reference system of the tile. Most often an ESPG code but
 66  any valid geopandas CRS specification is valid. Defaults to 3857 (i.e. Web
 67  Mercator)."""
 68  rotation: float = 0.0
 69  """cumulative rotation of the tileable."""
 70  debug: bool = False
 71  """if True prints debug messages. Defaults to False."""
 72
 73  # Tileable constructor called by subclasses - should not be used directly
 74  def __init__(self, **kwargs):
 75    for k, v in kwargs.items():
 76      self.__dict__[k] = v
 77    if self.debug:
 78      print(
 79        f"""Debugging messages enabled for Tileable (but there aren't
 80          any at the moment...)"""
 81      )
 82    self._setup_tiles()
 83    self.setup_vectors()
 84    self._setup_regularised_prototile()
 85    return
 86
 87
 88  def setup_vectors(self) -> None:
 89    """Sets up the symmetry translation vectors as floating point pairs
 90    indexed by integer tuples with respect to either a rectangular or
 91    triangular grid location.
 92
 93    Derived from the size and shape of the tile attribute. These are not
 94    the minimal translation vectors, but the 'face to face' vectors of the
 95    tile, such that a hexagonal tile will have 3 vectors, not the minimal
 96    parallelogram pair. Also supplies the inverse vectors.
 97
 98    The vectors are stored in a dictionary indexed by their
 99    coordinates, e.g.
100
101      {( 1,  0): ( 100, 0), ( 0,  1): (0,  100),
102       (-1,  0): (-100, 0), ( 0, -1): (0, -100)}
103
104    For a tileable of type `TileShape.HEXAGON`, the indexing tuples
105    have three components. See https://www.redblobgames.com/grids/hexagons/
106    """
107    t = self.prototile.geometry[0]
108    pts = [p for p in t.exterior.coords][:-1]
109    n_pts = len(pts)
110    vec_dict = {}
111    if n_pts == 4:
112      vecs = [(q[0] - p[0], q[1] - p[1])
113          for p, q in zip(pts, pts[1:] + pts[:1])]
114      i = [1, 0, -1,  0]
115      j = [0, 1,  0, -1]
116      vec_dict = {(i, j): v for i, j, v in zip(i, j, vecs)}
117    elif n_pts == 6:
118      vecs = [(q[0] - p[0], q[1] - p[1])
119          for p, q in zip(pts, pts[2:] + pts[:2])]
120      # hex grid coordinates associated with each of the vectors
121      i = [ 0,  1,  1,  0, -1, -1]
122      j = [ 1,  0, -1, -1,  0,  1]
123      k = [-1, -1,  0,  1,  1,  0]
124      vec_dict = {(i, j, k): v for i, j, k, v in zip(i, j, k, vecs)}
125    self.vectors = vec_dict
126
127
128  def get_vectors(
129      self, as_dict: bool = False
130    ) -> Union[dict[tuple[int], tuple[float]], list[tuple[float]]]:
131    """
132    Returns symmetry translation vectors as floating point pairs.
133    Optionally returns the vectors in a dictionary indexed by their
134    coordinates, e.g.
135
136      {( 1,  0): ( 100, 0), ( 0,  1): (0,  100),
137       (-1,  0): (-100, 0), ( 0, -1): (0, -100)}
138
139    Returns:
140      Union[ dict[tuple[int],tuple[float]], list[tuple[float]] ]:
141        either the vectors as a list of float tuples, or a dictionary
142        of those vectors indexed by integer coordinate tuples.
143    """
144    if as_dict:
145      return self.vectors
146    else:
147      return list(self.vectors.values())
148
149
150  # Make up a regularised tile by carefully unioning the tiles
151  def setup_regularised_prototile_from_tiles(self) -> None:
152    """Sets the regularised tile to a union of the tiles."""
153    self.regularised_prototile = copy.deepcopy(self.prototile)
154    self.regularised_prototile.geometry = [tiling_utils.safe_union(
155      self.tiles.geometry, as_polygon = True)]
156    # This simplification seems very crude but fixes all kinds of issues...
157    # particularly with the triaxial weave units... where intersection 
158    # operations are prone to creating spurious vertices, etc.
159    # self.regularised_prototile.geometry[0] = \
160    #   self.regularised_prototile.geometry[0].simplify(
161    #     self.spacing * tiling_utils.RESOLUTION)
162    return
163
164
165  def merge_fragments(self, fragments:list[geom.Polygon]) -> list[geom.Polygon]:
166    """
167    Merges a set of polygons based on testing if they touch when subjected
168    to the translation vectors provided by `get_vectors()`.
169
170    Called by `regularise_tiles()` to combine tiles in a tile unit that
171    may be fragmented as supplied but will combine after tiling into single
172    tiles. This step makes for more efficient implementation of the
173    tiling of map regions.
174
175    Args:
176      fragments (list[geom.Polygon]): A set of polygons to merge.
177
178    Returns:
179      list[geom.Polygon]: A minimal list of merged polygons.
180    """
181    if len(fragments) == 1:
182      return [f for f in fragments if not f.is_empty]
183    fragments = [f for f in fragments if not f.is_empty]
184    prototile = self.prototile.geometry[0]
185    reg_prototile = copy.deepcopy(self.regularised_prototile.geometry[0])
186    changes_made = True
187    while changes_made:
188      changes_made = False
189      for v in self.vectors.values():
190        # empty list to collect the new fragments
191        # assembled in this iteration
192        next_frags = []
193        t_frags = [affine.translate(f, v[0], v[1]) for f in fragments]
194        # build a set of any near matching pairs of
195        # fragments and their translated copies
196        matches = set()
197        for i, f1 in enumerate(fragments):
198          for j, f2, in enumerate(t_frags):
199            if i < j and tiling_utils.touch_along_an_edge(f1, f2):
200              matches.add((i, j))
201        # determine which of these when unioned has the larger area in common # with the prototile
202        frags_to_remove = set()
203        for i, j in matches:
204          f1, f2 = fragments[i], t_frags[j]
205          u1 = f1.buffer(tiling_utils.RESOLUTION, join_style = 2, cap_style = 3).union(
206            f2.buffer(tiling_utils.RESOLUTION, join_style = 2, cap_style = 3))
207          u2 = affine.translate(u1, -v[0], -v[1])
208          if prototile.intersection(u1).area > prototile.intersection(u2).area:
209            u1 = u1.buffer(-tiling_utils.RESOLUTION, join_style = 2, cap_style = 3)
210            u2 = u2.buffer(-tiling_utils.RESOLUTION, join_style = 2, cap_style = 3)
211            next_frags.append(u1)
212            reg_prototile = reg_prototile.union(u1).difference(u2)
213          else:
214            u1 = u1.buffer(-tiling_utils.RESOLUTION, join_style = 2, cap_style = 3)
215            u2 = u2.buffer(-tiling_utils.RESOLUTION, join_style = 2, cap_style = 3)
216            next_frags.append(u2)
217            reg_prototile = reg_prototile.union(u2).difference(u1)
218          changes_made = True
219          frags_to_remove.add(i)
220          frags_to_remove.add(j)
221        fragments = [f for i, f in enumerate(fragments)
222                     if not (i in frags_to_remove)] + next_frags
223    self.regularised_prototile.loc[0, "geometry"] = reg_prototile
224    # self.regularised_prototile.geometry[0] = reg_prototile
225    return [f for f in fragments if not f.is_empty] # don't return any duds
226
227
228  def reattach_tiles(self) -> None:
229    """Move tiles that are outside the regularised prototile main polygon
230    back inside it adjusting regularised prototile if needed.
231    """
232    reg_prototile = self.regularised_prototile.geometry[0]
233    new_reg_prototile = copy.deepcopy(reg_prototile)
234    new_tiles = list(self.tiles.geometry)
235    for i, p in enumerate(self.tiles.geometry):
236      if np.isclose(reg_prototile.intersection(p).area, p.area):
237        new_tiles[i] = p
238        continue
239      for v in self.vectors.values():
240        t_p = affine.translate(p, v[0], v[1])
241        if reg_prototile.intersects(t_p):
242          new_reg_prototile = new_reg_prototile.union(t_p)
243          new_tiles[i] = t_p
244    self.tiles.geometry = gpd.GeoSeries(new_tiles)
245    self.regularised_prototile.loc[0, "geometry"] = new_reg_prototile
246    # self.regularised_prototile.geometry[0] = new_reg_prototile
247    return None
248
249
250  def regularise_tiles(self) -> None:
251    """Combines separate tiles that share a tile_id value into
252    single tiles, if they would end up touching after tiling.
253
254    Also adjusts the `Tileable.regularised_prototile`
255    attribute accordingly.
256    """
257    self.regularised_prototile = copy.deepcopy(self.prototile)
258    # This preserves order while finding uniques, unlike list(set()).
259    # Reordering ids might cause confusion when colour palettes
260    # are not assigned explicitly to each id, but in the order
261    # encountered in the tile_id Series of the GeoDataFrame.
262    tiles, tile_ids = [], []
263    ids = list(self.tiles.tile_id.unique())
264    for id in ids:
265      fragment_set = list(
266        self.tiles[self.tiles.tile_id == id].geometry)
267      merge_result = self.merge_fragments(fragment_set)
268      tiles.extend(merge_result)
269      tile_ids.extend([id] * len(merge_result))
270
271    self.tiles = gpd.GeoDataFrame(
272      data = {"tile_id": tile_ids},
273      crs = self.crs,
274      geometry = gpd.GeoSeries([tiling_utils.get_clean_polygon(t) 
275                                for t in tiles]))
276
277    self.regularised_prototile = \
278      self.regularised_prototile.explode(ignore_index = True)
279    if self.regularised_prototile.shape[0] > 1:
280      self.regularised_prototile.geometry = tiling_utils.get_largest_polygon(
281        self.regularised_prototile.geometry)
282    return None
283
284
285  def get_local_patch(self, r: int = 1,
286                      include_0: bool = False) -> gpd.GeoDataFrame:
287    """Returns a GeoDataFrame with translated copies of the Tileable.
288
289    The geodataframe takes the same form as the `Tileable.tile` attribute.
290
291    Args:
292      r (int, optional): the number of 'layers' out from the unit to
293        which the translate copies will extendt. Defaults to `1`.
294      include_0 (bool, optional): If True includes the Tileable itself at
295        (0, 0). Defaults to `False`.
296
297    Returns:
298      gpd.GeoDataFrame: A GeoDataframe of the tiles extended by a number 
299        of 'layers'.
300    """
301    # a dictionary of all the vectors we need, starting with (0, 0)
302    vecs = (
303      {(0, 0, 0): (0, 0)}
304      if self.base_shape in (TileShape.HEXAGON,)
305      else {(0, 0): (0, 0)}
306    )
307    steps = r if self.base_shape in (TileShape.HEXAGON,) else r * 2
308    # a dictionary of the last 'layer' of added vectors
309    last_vecs = copy.deepcopy(vecs)
310    # get the translation vectors in a dictionary indexed by coordinates
311    # we keep track of the sum of vectors using the (integer) coordinates
312    # to avoid duplication of moves due to floating point inaccuracies
313    vectors = self.get_vectors(as_dict = True)
314    for i in range(steps):
315      new_vecs = {}
316      for k1, v1 in last_vecs.items():
317        for k2, v2 in vectors.items():
318          # add the coordinates to make a new key...
319          new_key = tuple([k1[i] + k2[i] for i in range(len(k1))])
320          # ... and the vector components to make a new value
321          new_val = (v1[0] + v2[0], v1[1] + v2[1])
322          # if we haven't reached here before store it
323          if not new_key in vecs:
324            new_vecs[new_key] = new_val
325      # extend the vectors and set the last layer to the set just added
326      vecs = vecs | new_vecs
327      last_vecs = new_vecs
328    if not include_0:  # throw away the identity vector
329      vecs.pop((0, 0, 0) if self.base_shape in (TileShape.HEXAGON,) else (0, 0))
330    ids, tiles = [], []
331    # we need to add the translated prototiles in order of their distance from # tile 0, esp. in the square case, i.e. something like this:
332    #
333    #      5 4 3 4 5
334    #      4 2 1 2 4
335    #      3 1 0 1 3
336    #      4 2 1 2 4
337    #      5 4 3 4 5
338    #
339    # this is important for topology detection, where filtering back to the
340    # local patch of radius 1 is greatly eased if prototiles have been added in 
341    # this order. We use the vector index tuples not the euclidean distances
342    # because this may be more resistant to odd effects for non-convex tiles
343    extent = self.prototile.geometry.scale(
344      2 * r + tiling_utils.RESOLUTION, 2 * r + tiling_utils.RESOLUTION,
345      origin = self.prototile.geometry[0].centroid)[0]
346    vector_lengths = {index: np.sqrt(sum([_ ** 2 for _ in index]))
347                      for index in vecs.keys()}
348    ordered_vector_keys = [k for k, v in sorted(vector_lengths.items(), 
349                                                key = lambda item: item[1])]
350    for k in ordered_vector_keys:
351      v = vecs[k]
352      if geom.Point(v[0], v[1]).within(extent):
353        ids.extend(self.tiles.tile_id)
354        tiles.extend(
355          self.tiles.geometry.apply(affine.translate, xoff = v[0], yoff = v[1]))
356    return gpd.GeoDataFrame(
357      data = {"tile_id": ids}, crs=self.crs,
358      geometry = tiling_utils.gridify(gpd.GeoSeries(tiles))
359    )
360
361
362  def fit_tiles_to_prototile(self, centre_tile: int = 0) -> None:
363    """Fits the tiles so they sit inside the prototile boundary.
364
365    If tiles project outside the boundaries of the prototile, this
366    method will clip them so that they don't. This may result in
367    'fragmented' tiles, i.e. pieces that would form a single tile
368    after tiling which are separated into fragments.
369
370    Args:
371      centre_tile (int, optional): the index position of the central
372        tile. Defaults to `0`.
373    """
374    dxy = self.tiles.geometry[centre_tile].centroid
375    self.tiles.geometry = self.tiles.translate(-dxy.x, -dxy.y)
376    # use r = 2 because rectangular tiles may need diagonal neighbours
377    patch = (
378      self.get_local_patch(r=2, include_0=True)
379      if self.base_shape in (TileShape.RECTANGLE,)
380      else self.get_local_patch(r=1, include_0=True)
381    )
382    self.tiles = patch.clip(self.prototile)
383    # repair any weirdness...
384    self.tiles.geometry = tiling_utils.repair_polygon(self.tiles.geometry)
385    self.tiles = self.tiles[self.tiles.geometry.area > 0]
386    self.regularised_prototile = copy.deepcopy(self.prototile)
387    return None
388
389
390  # applicable to both TileUnits and WeaveUnits
391  def inset_tiles(self, inset: float = 0) -> "Tileable":
392    """Returns a new Tileable with an inset applied around the tiles.
393
394    Works by applying a negative buffer of specfied size to all tiles.
395    Tiles that collapse to zero area are removed and the tile_id
396    attribute updated accordingly.
397
398    NOTE: this method is likely to not preserve the relative area of tiles.
399
400    Args:
401      inset (float, optional): The distance to inset. Defaults to `0`.
402
403    Returns:
404      "Tileable": the new inset Tileable.
405    """
406    inset_tiles, inset_ids = [], []
407    for p, id in zip(self.tiles.geometry, self.tiles.tile_id):
408      b = p.buffer(-inset, join_style = 2, cap_style = 3)
409      if not b.area <= 0:
410        inset_tiles.append(b)
411        inset_ids.append(id)
412    result = copy.deepcopy(self)
413    result.tiles = gpd.GeoDataFrame(
414      data={"tile_id": inset_ids},
415      crs=self.crs,
416      geometry=gpd.GeoSeries(inset_tiles),
417    )
418    return result
419
420
421  def plot(self, ax = None, show_prototile: bool = True, 
422    show_reg_prototile: bool = True, show_ids: str = "tile_id",
423    show_vectors: bool = False, r: int = 0, prototile_edgecolour: str = "k", 
424    reg_prototile_edgecolour: str = "r", r_alpha: float = 0.3, 
425    cmap: list[str] = None, figsize: tuple[float] = (8, 8), **kwargs) -> pyplot.axes:
426    """Plots a representation of the Tileable on the supplied axis. **kwargs
427    are passed on to matplotlib.plot()
428
429    Args:
430      ax (_type_, optional): matplotlib axis to draw to. Defaults to None.
431      show_prototile (bool, optional): if `True` show the tile outline.
432        Defaults to `True`.
433      show_reg_prototile (bool, optional): if `True` show the regularised tile
434        outline. Defaults to `True`.
435      show_ids (str, optional): if `tile_id` show the tile_ids. If
436        `id` show index number. If None or `''` don't label tiles.
437        Defaults to `tile_id`.
438      show_vectors (bool, optional): if `True` show the translation
439        vectors (not the minimal pair, but those used by
440        `get_local_patch()`). Defaults to `False`.
441      r (int, optional): passed to `get_local_patch()` to show context if
442        greater than 0. Defaults to `0`.
443      r_alpha (float, optional): alpha setting for units other than the
444        central one. Defaults to 0.3.
445      prototile_edgecolour (str, optional): outline colour for the tile.
446        Defaults to `"k"`.
447      reg_prototile_edgecolour (str, optional): outline colour for the
448        regularised. Defaults to `"r"`.
449      cmap (list[str], optional): colour map to apply to the central
450        tiles. Defaults to `None`.
451      figsize (tuple[float], optional): size of the figure.
452        Defaults to `(8, 8)`.
453    
454    Returns:
455      pyplot.axes: to which calling context may add things.
456    """
457    w = self.prototile.geometry[0].bounds[2] - \
458      self.prototile.geometry[0].bounds[0]
459    n_cols = len(set(self.tiles.tile_id))
460    if cmap is None:
461      cm = "Dark2" if n_cols <= 8 else "Paired"
462    else:
463      cm = cmap
464    if ax is None:
465      ax = self.tiles.plot(
466        column="tile_id", cmap=cm, figsize=figsize, **kwargs)
467    else:
468      self.tiles.plot(
469        ax=ax, column="tile_id", cmap=cm, figsize=figsize, **kwargs)
470    if show_ids != None and show_ids != "":
471      do_label = True
472      if show_ids == "tile_id" or show_ids == True:
473        labels = self.tiles.tile_id
474      elif show_ids == "id":
475        labels = [str(i) for i in range(self.tiles.shape[0])]
476      else:
477        do_label = False
478      if do_label:
479        for id, tile in zip(labels, self.tiles.geometry):
480          ax.annotate(id, (tile.centroid.x, tile.centroid.y),
481            ha = "center", va = "center", bbox = {"lw": 0, "fc": "#ffffff40"})
482    if r > 0:
483      self.get_local_patch(r=r).plot(
484        ax = ax, column = "tile_id", alpha = r_alpha, cmap = cm, **kwargs)
485    if show_prototile:
486      self.prototile.plot(ax = ax, ec = prototile_edgecolour, lw = 0.5, 
487                          fc = "#00000000", **kwargs)
488    if show_vectors:  # note that arrows in mpl are dimensioned in plotspace
489      vecs = self.get_vectors()
490      for v in vecs[: len(vecs) // 2]:
491        ax.arrow(0, 0, v[0], v[1], color = "k", width = w * 0.002,
492          head_width = w * 0.05, length_includes_head = True, zorder = 3)
493    if show_reg_prototile:
494      self.regularised_prototile.plot(
495        ax = ax, ec = reg_prototile_edgecolour, fc = "#00000000", 
496        lw = 1.5, zorder = 2, **kwargs)
497    return ax
498
499
500  def _get_legend_tiles(self):
501    """Returns the tiles augmented by a rotation column.
502
503    This base implementation may be overridden by specific tile unit types.
504    In particular see
505    `weavingspace.weave_unit.WeaveUnit._get_legend_tiles()`.
506    """
507    tiles = copy.deepcopy(self.tiles)
508    tiles["rotation"] = 0
509    return tiles
510
511
512  def transform_scale(self, xscale: float = 1.0, yscale: float = 1.0) -> "Tileable":
513    """Transforms tileable by scaling.
514
515    Args:
516      xscale (float, optional): x scale factor. Defaults to 1.0.
517      yscale (float, optional): y scale factor. Defaults to 1.0.
518
519    Returns:
520      Tileable: the transformed Tileable.
521    """
522    result = copy.deepcopy(self)
523    result.tiles.geometry = tiling_utils.gridify(
524      self.tiles.geometry.scale(xscale, yscale, origin=(0, 0)))
525    result.prototile.geometry = tiling_utils.gridify(
526      self.prototile.geometry.scale(xscale, yscale, origin=(0, 0)))
527    result.regularised_prototile.geometry = tiling_utils.gridify(
528      self.regularised_prototile.geometry.scale(xscale, yscale, origin=(0, 0)))
529    result.setup_vectors()
530    return result
531
532
533  def transform_rotate(self, angle: float = 0.0) -> "Tileable":
534    """Transforms tiling by rotation.
535
536    Args:
537      angle (float, optional): angle to rotate by. Defaults to 0.0.
538
539    Returns:
540      Tileable: the transformed Tileable.
541    """
542    result = copy.deepcopy(self)
543    result.tiles.geometry = tiling_utils.gridify(
544      self.tiles.geometry.rotate(angle, origin=(0, 0)))
545    result.prototile.geometry = tiling_utils.gridify(
546      self.prototile.geometry.rotate(angle, origin=(0, 0)))
547    result.regularised_prototile.geometry = tiling_utils.gridify(
548      self.regularised_prototile.geometry.rotate(angle, origin=(0, 0)))
549    result.setup_vectors()
550    result.rotation = result.rotation + angle
551    return result
552
553
554  def transform_skew(self, xa: float = 0.0, ya: float = 0.0) -> "Tileable":
555    """Transforms tiling by skewing
556
557    Args:
558      xa (float, optional): x direction skew. Defaults to 0.0.
559      ya (float, optional): y direction skew. Defaults to 0.0.
560
561    Returns:
562      Tileable: the transformed Tileable.
563    """
564    result = copy.deepcopy(self)
565    result.tiles.geometry = tiling_utils.gridify(
566      self.tiles.geometry.skew(xa, ya, origin=(0, 0)))
567    result.prototile.geometry = tiling_utils.gridify(
568      self.prototile.geometry.skew(xa, ya, origin=(0, 0)))
569    result.regularised_prototile.geometry = tiling_utils.gridify(
570      self.regularised_prototile.geometry.skew(xa, ya, origin=(0, 0)))
571    result.setup_vectors()
572    return result