weavingspace.tile_map

  1#!/usr/bin/env python
  2# coding: utf-8
  3
  4"""Classes for tiling maps. `weavingspace.tile_map.Tiling` and
  5`weavingspace.tile_map.TiledMap` are exposed in the  public API and
  6respectively enable creation of a tiling and plotting of the tiling as a
  7multivariate map.
  8"""
  9
 10from dataclasses import dataclass
 11from typing import Union
 12import itertools
 13import copy
 14
 15import numpy as np
 16import geopandas as gpd
 17import pandas as pd
 18
 19from matplotlib.figure import Figure
 20import matplotlib.colors
 21import matplotlib.pyplot as pyplot
 22
 23import shapely.geometry as geom
 24import shapely.affinity as affine
 25import shapely.ops
 26
 27from weavingspace.tileable import Tileable
 28from weavingspace.tileable import TileShape
 29from weavingspace.tile_unit import TileUnit
 30
 31from weavingspace import tiling_utils
 32
 33from time import perf_counter
 34
 35@dataclass
 36class _TileGrid():
 37  """A class to represent the translation centres of a tiling.
 38
 39  We store the grid as a GeoSeries of Point objects to make it
 40  simple to plot it in map views if required.
 41
 42  Implementation relies on transforming the translation vectors into a square
 43  space where the tile spacing is unit squares, then transforming this back
 44  into the original map space. Some member variables of the class are in
 45  the transformed grid generation space, some in the map space.
 46  """
 47  tile:TileUnit = None
 48  """the base tile in map space."""
 49  to_tile:gpd.GeoSeries = None
 50  """geometry of the region to be tiled in map space."""
 51  transform:tuple[float] = None
 52  """the forward transformation from map space to tiling grid generation
 53  space, stored as a shapely.affinity transform tuple of 6 floats."""
 54  inverse_transform:tuple[float] = None
 55  """the inverse transform from tiling grid space to map space."""
 56  centre:tuple[float] = None
 57  """centre point of the extent in map sapce."""
 58  points:gpd.GeoSeries = None
 59  """geom.Points recording the translation vectors of the tiling in map space.
 60  """
 61  _extent:gpd.GeoSeries = None
 62  """geometry of the circular extent of the tiling transformed into tiling 
 63  grid generation space."""
 64  _at_centroids:bool = False
 65  """if True the grid will consist of the centroids of the spatial units in the 
 66  to_tile region, allowing a simple way to use a tile unit as a point symbol."""
 67
 68  def __init__(self, tile:TileUnit, to_tile:gpd.GeoSeries,
 69               at_centroids:bool = False):
 70    self.tile = tile
 71    self.to_tile = self._get_area_to_tile(to_tile)
 72    self.inverse_transform, self.transform = self._get_transforms()
 73    self.extent, self.centre = self._get_extent()
 74    self._at_centroids = at_centroids
 75    if self._at_centroids:
 76      self.points = to_tile.centroid
 77    else:
 78      self.points = self._get_grid()
 79    self.points.crs = self.tile.crs
 80    self.points = tiling_utils.gridify(self.points)
 81
 82
 83  def _get_area_to_tile(self, to_tile) -> geom.Polygon:
 84    bb = to_tile.total_bounds
 85    poly = tiling_utils.gridify(
 86      geom.Polygon(((bb[0], bb[1]), (bb[2], bb[1]),
 87                    (bb[2], bb[3]), (bb[0], bb[3]))))
 88    return gpd.GeoSeries([poly])
 89
 90
 91  def _get_extent(self) -> tuple[gpd.GeoSeries, geom.Point]:
 92    """Returns the extent and centre of the grid.
 93
 94    Extent is in the grid-generation space.
 95
 96    Returns:
 97      tuple[gpd.GeoSeries, geom.Point]: the extent of the grid and its
 98        centre.
 99    """
100
101    # TODO: the minimum_rotate_rectangle seems to throw an error?
102    mrr = self.to_tile[0].minimum_rotated_rectangle
103    mrr_centre = geom.Point(mrr.centroid.coords[0])
104    mrr_corner = geom.Point(mrr.exterior.coords[0])
105    radius = mrr_centre.distance(mrr_corner)
106    ext = tiling_utils.gridify(
107      affine.affine_transform(mrr_centre.buffer(radius), self.transform))
108    return gpd.GeoSeries([ext]), (mrr_centre.x, mrr_centre.y)
109
110
111  def _get_transforms(self) -> tuple[float]:
112    """Returns the forward and inverse transforms from map space to
113    grid generation space.
114
115    In grid generation space the translation vectors are (1, 0) and (0, 1)
116    so we can simply form a matrix from two translation vectors and invert
117    it to get the forward transform. The inverse transform is the matrix
118    formed from the vectors.
119
120    Returns:
121      tuple[float]: shapely affine transform tuples (a, b, d, e, dx, dy).
122        See https://shapely.readthedocs.io/en/stable/manual.html#affine-transformations for details.
123    """
124    v = self.tile.get_vectors()
125    vector_array = np.array([[v[0][0], v[1][0]],
126                             [v[0][1], v[1][1]]])
127    inv_tfm = np.linalg.inv(vector_array)
128    return (self._np_to_shapely_transform(vector_array),
129        self._np_to_shapely_transform(inv_tfm))
130
131
132  def _get_grid(self) -> gpd.GeoSeries:
133    """Generates the grid transformed into map space
134
135    Obtain dimensions of the transformed region, then set down a uniform
136    grid.
137
138    Returns:
139      gpd.GeoSeries: the grid as a collection of geom.Points.
140    """
141    _w, _h, _l, _b = tiling_utils.get_width_height_left_bottom(self.extent)
142    w = int(np.ceil(_w))
143    h = int(np.ceil(_h))
144    l = _l - (w - _w) / 2
145    b = _b - (h - _h) / 2
146    mesh = np.array(np.meshgrid(np.arange(w) + l,
147                  np.arange(h) + b)).reshape((2, w * h)).T
148    pts = [tiling_utils.gridify(geom.Point(p[0], p[1])) for p in mesh]
149    return gpd.GeoSeries([p for p in pts if p.within(self.extent[0])]) \
150      .affine_transform(self.inverse_transform)
151
152
153  def _np_to_shapely_transform(self, mat:np.ndarray) -> tuple[float]:
154    """Converts a numpy affine transform matrix to shapely format.
155
156    Args:
157      mat (np.ndarray): numpy array to convert.
158
159    Returns:
160      tuple[float]: shapely affine transform tuple
161    """
162    return (mat[0][0], mat[0][1], mat[1][0], mat[1][1], 0, 0)
163
164
165@dataclass
166class Tiling:
167  """Class that applies a `Tileable` object to a region to be mapped.
168
169  The result of the tiling procedure is stored in the `tiles` variable and
170  covers a region sufficient that the tiling can be rotated to any desired
171  angle.
172  """
173  tile_unit:Tileable = None
174  """tileable on which the tiling is based."""
175  tile_shape:TileShape = None
176  """base shape of the tileable."""
177  region:gpd.GeoDataFrame = None
178  """the region to be tiled."""
179  region_union: geom.Polygon = None
180  grid:_TileGrid = None
181  """the grid which will be used to apply the tiling."""
182  tiles:gpd.GeoDataFrame = None
183  """the tiles after tiling has been carried out."""
184  prototiles:gpd.GeoDataFrame = None
185  """the prototiles after tiling has been carried out."""
186  rotation:float = 0.0
187  """the cumulative rotation already applied to the tiling."""
188
189  def __init__(self, unit:Tileable, region:gpd.GeoDataFrame, id_var = None,
190         prototile_margin:float = 0, tiles_sf:float = 1,
191         tiles_margin:float = 0, as_icons:bool = False) -> None:
192    """Class to persist a tiling by filling an area relative to
193    a region sufficient to apply the tiling at any rotation.
194
195    The Tiling constructor allows a number of adjustments to the supplied
196    `weavingspace.tileable.Tileable` object:
197
198    + `prototile_margin` values greater than 0 will introduce spacing of
199    the specified distance between tiles on the boundary of each tile
200    by applying the `TileUnit.inset_prototile()` method. Note that this
201    operation does not make sense for `WeaveUnit` objects,
202    and may not preserve the equality of tile areas.
203    + `tiles_sf` values less than one scale down tiles by applying the 
204    `TileUnit.scale_tiles()` method. Does not make sense for `WeaveUnit` 
205    objects.
206    + `tiles_margin` values greater than one apply a negative buffer of
207    the specified distance to every tile in the tiling by applying the
208    `Tileable.inset_tiles()` method. This option is applicable to both
209    `WeaveUnit` and `TileUnit` objects.
210
211    Args:
212      unit (Tileable): the tile_unit to use.
213      region (gpd.GeoDataFrame): the region to be tiled.
214      prototile_margin (float, optional): values greater than 0 apply an
215        inset margin to the tile unit. Defaults to 0.
216      tiles_sf (float, optional): scales the tiles. Defaults to 1.
217      tiles_margin (float, optional): applies a negative buffer to
218        the tiles. Defaults to 0.
219      as_icons (bool, optional): if True prototiles will only be placed at
220        the region's zone centroids, one per zone. Defaults to
221        False.
222    """
223    self.tile_unit = unit
224    self.rotation = self.tile_unit.rotation
225    if tiles_margin > 0:
226      self.tile_unit = self.tile_unit.inset_tiles(tiles_margin)
227    if tiles_sf != 1:
228      if isinstance(self.tile_unit, TileUnit):
229        self.tile_unit = self.tile_unit.scale_tiles(tiles_sf)
230      else:
231        print(f"""Applying scaling to tiles of a WeaveUnit does not make sense. 
232              Ignoring tiles_sf setting of {tiles_sf}.""")
233    if prototile_margin > 0:
234      if isinstance(self.tile_unit, TileUnit):
235        self.tile_unit = self.tile_unit.inset_prototile(prototile_margin)
236      else:
237        print(f"""Applying a prototile margin to a WeaveUnit does 
238              not make sense. Ignoring prototile_margin setting of
239              {prototile_margin}.""")
240    self.region = region
241    self.region.sindex
242    self.region_union = self.region.geometry.unary_union
243    if id_var != None:
244      print("""id_var is no longer required and will be deprecated soon.
245            A temporary unique index attribute is added and removed when 
246            generating the tiled map.""")
247    if as_icons:
248      self.grid = _TileGrid(self.tile_unit, self.region.geometry, True)
249    else:
250      self.grid = _TileGrid(self.tile_unit, self.region.geometry)
251    self.tiles, self.prototiles = self.make_tiling()
252    self.tiles.sindex
253
254
255  def get_tiled_map(self, rotation:float = 0.,
256                    join_on_prototiles:bool = True,
257                    prioritise_tiles:bool = True,
258                    ragged_edges:bool = True,
259                    use_centroid_lookup_approximation = False,
260                    debug = False) -> "TiledMap":
261    """Returns a `TiledMap` filling a region at the requested rotation.
262
263    HERE BE DRAGONS! This function took a lot of trial and error to get
264    right, so modify with CAUTION!
265
266    The `proritise_tiles = True` option means that the tiling will not
267    break up the tiles in `TileUnit`s at the boundaries between areas
268    in the mapped region, but will instead ensure that tiles remain
269    complete, picking up their data from the region zone which they overlap
270    the most.
271    
272    The exact order in which operations are performed affects performance.
273    For example, the final clipping to self.region when ragged_edges =
274    False is _much_ slower if it is carried out before the dissolving of
275    tiles into the region zones. So... again... modify CAREFULLY!
276
277    Args:
278      rotation (float, optional): An optional rotation to apply. Defaults
279        to 0.
280      join_on_prototiles (bool, optional): if True data from the region
281        dataset are joined to tiles based on the prototile to which they
282        belong. If False the join is based on the tiles in relation to the
283        region areas. For weave-based tilings False is probably to be
284        preferred. Defaults to True.
285      prioritise_tiles (bool, optional): if True tiles will not be
286        broken at boundaries in the region dataset. Defaults to True.
287      ragged_edges (bool, optional): if True tiles at the edge of the
288        region will not be cut by the region extent - ignored if
289        prioritise_tiles is False when edges will always be clipped to
290        the region extent. Defaults to True.
291      use_centroid_lookup_approximation (bool, optional): if True use
292        tile centroids for lookup of region data - ignored if
293        prioritise_tiles is False when it is irrelevant. Defaults to
294        False.
295      debug (bool, optional): if True prints timing messages. Defaults
296        to False.
297
298    Returns:
299      TiledMap: a TiledMap of the source region.
300    """
301    if debug:
302      t1 = perf_counter()
303
304    id_var = self._setup_region_DZID()
305    if join_on_prototiles:
306      tiled_map, join_layer = self.rotated(rotation)
307      tiled_map["joinUID"] = self.tiles["prototile_id"]
308    else:
309      tiled_map = self.rotated(rotation)[0]
310      tiled_map["joinUID"] = self.tiles["tile_id"]
311      join_layer = tiled_map
312    join_layer["joinUID"] = list(range(join_layer.shape[0]))
313
314    # compile a list of the variable names we are NOT going to change
315    # i.e. everything except the geometry and the id_var
316    region_vars = list(self.region.columns)
317    region_vars.remove("geometry")
318    region_vars.remove(id_var)
319
320    if debug:
321      t2 = perf_counter()
322      print(f"STEP 1: prep data (rotation if requested): {t2 - t1:.3f}")
323
324    if prioritise_tiles:  # maintain tile continuity across zone boundaries
325      # select only tiles inside a spacing buffer of the region
326      # make column with unique ID for every tile in the tiling
327      # the join ID is unique per tile
328      # if join_on_prototiles:
329      #   tiled_map["joinUID"] = self.tiles["prototile_id"]
330      # else:
331      #   tiled_map["joinUID"] = self.tiles["tile_id"]
332 
333      if use_centroid_lookup_approximation:
334        t5 = perf_counter()
335        tile_pts = copy.deepcopy(join_layer)
336        tile_pts.geometry = tile_pts.centroid
337        lookup = tile_pts.sjoin(
338          self.region, how = "inner")[["joinUID", id_var]]
339      else:
340        # determine areas of overlapping tiles and drop the data we join the 
341        # data back later, so dropping makes that easier overlaying in region.
342        # overlay(tiles) seems to be faster??
343        # TODO: also... this part is performance-critical, think about fixes -- 
344        # possibly including the above centroid-based approx
345        overlaps = self.region.overlay(join_layer, make_valid = False)
346        # overlaps = self.region.overlay(tiled_map, make_valid = False)
347        if debug:
348          t3 = perf_counter()
349          print(f"STEP A2: overlay zones with tiling: {t3 - t2:.3f}")
350        overlaps["area"] = overlaps.geometry.area
351        if debug:
352          t4 = perf_counter()
353          print(f"STEP A3: calculate areas: {t4 - t3:.3f}")
354        overlaps.drop(columns = region_vars, inplace = True)
355        if debug:
356          t5 = perf_counter()
357          print(f"STEP A4: drop columns prior to join: {t5 - t4:.3f}")
358        # make a lookup by largest area tile to region id
359        lookup = overlaps \
360          .iloc[overlaps.groupby("joinUID")["area"] \
361          .agg(pd.Series.idxmax)][["joinUID", id_var]]
362      # now join the lookup and from there the region data
363      if debug:
364        t6 = perf_counter()
365        print(f"STEP A5: build lookup for join: {t6 - t5:.3f}")
366      tiled_map = tiled_map \
367        .merge(lookup, on = "joinUID") \
368        .merge(self.region.drop(columns = ["geometry"]), on = id_var)
369      if debug:
370        t7 = perf_counter()
371        print(f"STEP A6: perform lookup join: {t7 - t6:.3f}")
372      tiled_map.drop(columns = ["joinUID"], inplace = True)
373
374    else:  # here we overlay
375      tiled_map = self.region.overlay(tiled_map)
376      t7 = perf_counter()
377      if debug:
378        print(f"STEP B2: overlay tiling with zones: {t7 - t2:.3f}")
379
380    if join_on_prototiles:
381      tiled_map = tiled_map.loc[
382        shapely.intersects(self.region_union, np.array(tiled_map.geometry)), :]
383
384    tiled_map.drop(columns = [id_var], inplace = True)
385    self.region.drop(columns = [id_var], inplace = True)
386
387    # if we've retained tiles and want 'clean' edges, then clip
388    # note that this step is slow: geopandas unary_unions the clip layer
389    if prioritise_tiles and not ragged_edges:
390      tiled_map.sindex
391      tiled_map = tiled_map.clip(self.region)
392      if debug:
393        print(f"""STEP A7/B3: clip map to region: {perf_counter() - t7:.3f}""")
394
395    tm = TiledMap()
396    tm.tiling = self
397    tm.map = tiled_map
398    return tm
399
400
401  def _setup_region_DZID(self) -> str:
402    """Creates a new guaranteed-unique attribute in the self.region
403    dataframe, and returns its name.
404
405    Avoids a name clash with any existing attribute in the dataframe.
406
407    Returns:
408      str: name of the added attribute.
409    """
410    dzid = "DZID"
411    i = 0
412    while dzid in self.region.columns:
413      dzid = "DZID" + str(i)
414      i = i + 1
415    self.region[dzid] = list(range(self.region.shape[0]))
416    return dzid
417
418
419  def _rotate_gdf_to_geoseries(
420      self, gdf:gpd.GeoDataFrame,
421      angle:float, centre:tuple = (0, 0)
422    ) -> tuple[gpd.GeoSeries, tuple[float]]:
423    """Rotates the geometries in a GeoDataFrame as a single collection.
424
425    Rotation is about the supplied centre or about the centroid of the
426    GeoDataFrame (if not). This allows for reversal of  a rotation. [Note
427    that this might not be a required precaution!]
428
429    Args:
430      gdf (geopandas.GeoDataFrame): GeoDataFrame to rotate
431      angle (float): angle of rotation (degrees).
432      centre (tuple, optional): desired centre of rotation. Defaults
433        to (0, 0).
434
435    Returns:
436      tuple: a geopandas.GeoSeries and a tuple (point) of the centre of
437        the rotation.
438    """
439    centre = (
440      gdf.geometry.unary_union.centroid.coords[0]
441      if centre is None
442      else centre)
443    return gdf.geometry.rotate(angle, origin = centre), centre
444
445
446  def make_tiling(self) -> gpd.GeoDataFrame:
447    """Tiles the region with a tile unit, returning a GeoDataFrame
448
449    Returns:
450      geopandas.GeoDataFrame: a GeoDataFrame of the region tiled with the
451        tile unit.
452    """
453    # we assume the geometry column is called geometry so make it so...
454    if self.region.geometry.name != "geometry":
455      self.region.rename_geometry("geometry", inplace = True)
456
457    # chain list of lists of GeoSeries geometries to list of geometries
458    tiles = itertools.chain(*[
459      self.tile_unit.tiles.geometry.translate(p.x, p.y)
460      for p in self.grid.points])
461    prototiles = itertools.chain(*[
462      self.tile_unit.prototile.geometry.translate(p.x, p.y)
463      for p in self.grid.points])
464    # replicate the tile ids
465    prototile_ids = list(range(len(self.grid.points)))
466    tile_ids = list(self.tile_unit.tiles.tile_id) * len(self.grid.points)
467    tile_prototile_ids = sorted(prototile_ids * self.tile_unit.tiles.shape[0])
468    tiles_gs = gpd.GeoSeries(tiles)
469    prototiles_gs = gpd.GeoSeries(prototiles)
470    # assemble and return as GeoDataFrames
471    tiles_gdf = gpd.GeoDataFrame(
472      data = {"tile_id": tile_ids, "prototile_id": tile_prototile_ids},
473      geometry = tiles_gs, crs = self.tile_unit.crs)
474    prototiles_gdf = gpd.GeoDataFrame(
475      data = {"prototile_id": prototile_ids},
476      geometry = prototiles_gs, crs = self.tile_unit.crs)
477    # unclear if we need the gridify or not...
478    return (tiling_utils.gridify(tiles_gdf),
479            tiling_utils.gridify(prototiles_gdf))
480
481
482  def rotated(self, rotation:float = None) -> gpd.GeoDataFrame:
483    """Returns the stored tiling rotated.
484
485    Args:
486      rotation (float, optional): Rotation angle in degrees.
487        Defaults to None.
488
489    Returns:
490      gpd.GeoDataFrame: Rotated tiling.
491    """
492    if self.tiles is None:
493      self.tiles = self.make_tiling()
494    self.rotation = rotation
495    if self.rotation == 0:
496      return self.tiles, self.prototiles
497    tiles = gpd.GeoDataFrame(
498      data = {"tile_id": self.tiles.tile_id,
499              "prototile_id": self.tiles.tile_id},
500      crs = self.tiles.crs,
501      geometry = tiling_utils.gridify(
502        self.tiles.geometry.rotate(rotation, origin = self.grid.centre)))
503    prototiles = gpd.GeoDataFrame(
504      data = {"prototile_id": self.prototiles.prototile_id},
505      crs = self.prototiles.crs,
506      geometry = tiling_utils.gridify(
507        self.prototiles.geometry.rotate(rotation, origin = self.grid.centre)))
508    return tiles, prototiles
509
510
511@dataclass
512class TiledMap:
513  """Class representing a tiled map. Should not be accessed directly, but
514  will be created by calling `Tiling.get_tiled_map()`. After creation the
515  variables and colourmaps attributes can be set, and then
516  `TiledMap.render()` called to make a map. Settable attributes are explained
517  in documentation of the `TiledMap.render()` method.
518
519  Examples:
520    Recommended usage is as follows. First, make a `TiledMap` from a `Tiling` object.
521
522      tm = tiling.get_tiled_map(...)
523
524    Some options in the `Tiling` constructor affect the map appearance. See
525    `Tiling` for details.
526
527    Once a `TiledMap` object exists, set options on it, either when calling
528    `TiledMap.render()` or explicitly, i.e.
529
530      tm.render(opt1 = val1, opt2 = val2, ...)
531
532    or
533
534      tm.opt1 = val1
535      tm.opt2 = val2
536      tm.render()
537
538    Option settings are persistent, i.e. unless a new `TiledMap` object is
539    created the option settings have to be explicitly reset to default
540    values on subsequent calls to `TiledMap.render()`.
541
542    The most important options are the `variables` and `colourmaps`
543    settings.
544
545    `variables` is a dictionary mapping `weavingspace.tileable.Tileable`
546    tile_ids (usually "a", "b", etc.) to variable names in the data. For
547    example,
548
549      tm.variables = dict(zip(["a", "b"], ["population", "income"]))
550
551    `colourmaps` is a dictionary mapping dataset variable names to the
552    matplotlib colourmap to be used for each. For example,
553
554      tm.colourmaps = dict(zip(tm.variables.values(), ["Reds", "Blues"]))
555
556    See [this notebook](https://github.com/DOSull/weaving-space/blob/main/weavingspace/example-tiles-cairo.ipynb)
557    for simple usage.
558    TODO: This more complicated example shows how categorical maps can be
559    created.
560  """
561  # these will be set at instantion by Tiling.get_tiled_map()
562  tiling:Tiling = None
563  """the Tiling with the required tiles"""
564  map:gpd.GeoDataFrame = None
565  """the GeoDataFrame on which this map is based"""
566  variables:dict[str,str] = None 
567  """lookup from tile_id to variable names"""
568  colourmaps:dict[str,Union[str,dict]] = None
569  """lookup from variables to matplotlib cmaps"""
570
571  # the below parameters can be set either before calling self.render()
572  # or passed in as parameters to self.render()
573  # these are solely TiledMap.render() options
574  legend:bool = True
575  """whether or not to show a legend"""
576  legend_zoom:float = 1.0
577  """<1 zooms out from legend to show more context"""
578  legend_dx:float = 0.
579  """x shift of legend relative to the map"""
580  legend_dy:float = 0.
581  """y shift of legend relative to the map"""
582  use_ellipse:bool = False
583  """if True clips legend with an ellipse"""
584  ellipse_magnification:float = 1.0
585  """magnification to apply to clip ellipse"""
586  radial_key:bool = False
587  """if True use radial key even for ordinal/ratio data (normally these will be 
588  shown by concentric tile geometries)"""
589  draft_mode:bool = False
590  """if True plot the map coloured by tile_id"""
591
592  # the parameters below are geopandas.plot options which we intercept to
593  # ensure they are applied appropriately when we plot a GDF
594  scheme:str = "equalinterval"
595  """geopandas scheme to apply"""
596  k:int = 100
597  """geopandas number of classes to apply"""
598  figsize:tuple[float] = (20, 15)
599  """maptlotlib figsize"""
600  dpi:float = 72
601  """dpi for bitmap formats"""
602
603  def render(self, **kwargs) -> Figure:
604    """Renders the current state to a map.
605
606    Note that TiledMap objects will usually be created by calling
607    `Tiling.get_tiled_map()`.
608
609    Args:
610      variables (dict[str,str]): Mapping from tile_id values to
611        variable names. Defaults to None.
612      colourmaps (dict[str,Union[str,dict]]): Mapping from variable
613        names to colour map, either a colour palette as used by
614        geopandas/matplotlib, a fixed colour, or a dictionary mapping
615        categorical data values to colours. Defaults to None.
616      legend (bool): If True a legend will be drawn. Defaults to True.
617      legend_zoom (float): Zoom factor to apply to the legend. Values <1
618        will show more of the tile context. Defaults to 1.0.
619      legend_dx (float): x shift to apply to the legend position.
620        Defaults to 0.0.
621      legend_dy (float): x and y shift to apply to the legend position.
622        Defaults to 0.0.
623      use_ellipse (bool): If True applies an elliptical clip to the
624        legend. Defaults to False.
625      ellipse_magnification (float): Magnification to apply to ellipse
626        clipped legend. Defaults to 1.0.
627      radial_key (bool): If True legend key for TileUnit maps will be
628        based on radially dissecting the tiles. Defaults to False.
629      draft_mode (bool): If True a map of the tiled map coloured by
630        tile_ids (and with no legend) is returned. Defaults to False.
631      scheme (str): passed to geopandas.plot for numeric data. Defaults to
632        "equalinterval".
633      k (int): passed to geopandas.plot for numeric data. Defaults to 100.
634      figsize (tuple[float,floar]): plot dimensions passed to geopandas.
635        plot. Defaults to (20,15).
636      dpi (float): passed to pyplot.plot. Defaults to 72.
637      **kwargs: other settings to pass to pyplot/geopandas.plot.
638
639    Returns:
640      matplotlib.figure.Figure: figure on which map is plotted.
641    """
642    pyplot.rcParams['pdf.fonttype'] = 42
643    pyplot.rcParams['pdf.use14corefonts'] = True
644    matplotlib.rcParams['pdf.fonttype'] = 42
645
646    to_remove = set()  # keep track of kwargs we use to setup TiledMap
647    for k, v in kwargs.items():
648      if k in self.__dict__:
649        self.__dict__[k] = v
650        to_remove.add(k)
651    # remove them so we don't pass them on to pyplot and get errors
652    for k in to_remove:
653      del kwargs[k]
654
655    if self.draft_mode:
656      fig = pyplot.figure(figsize = self.figsize)
657      ax = fig.add_subplot(111)
658      self.map.plot(ax = ax, column = "tile_id", cmap = "tab20",
659              **kwargs)
660      return fig
661
662    if self.legend:
663      # this sizing stuff is rough and ready for now, possibly forever...
664      reg_w, reg_h, *_ = \
665        tiling_utils.get_width_height_left_bottom(self.map.geometry)
666      tile_w, tile_h, *_ = \
667        tiling_utils.get_width_height_left_bottom(
668          self.tiling.tile_unit._get_legend_tiles().rotate(
669            self.tiling.rotation, origin = (0, 0)))
670      sf_w, sf_h = reg_w / tile_w / 3, reg_h / tile_h / 3
671      gskw = {"height_ratios": [sf_h * tile_h, reg_h - sf_h * tile_h],
672              "width_ratios": [reg_w, sf_w * tile_w]}
673
674      fig, axes = pyplot.subplot_mosaic(
675        [["map", "legend"], ["map", "."]],
676        gridspec_kw = gskw, figsize = self.figsize,
677        layout = "constrained", **kwargs)
678    else:
679      fig, axes = pyplot.subplots(
680        1, 1, figsize = self.figsize,
681        layout = "constrained", **kwargs)
682
683    if self.variables is None:
684      # get any floating point columns available
685      default_columns = \
686        self.map.select_dtypes(
687          include = ("float64", "int64")).columns
688      self.variables = dict(zip(self.map.tile_id.unique(),
689                                list(default_columns)))
690      print(f"""No variables specified, picked the first
691            {len(self.variables)} numeric ones available.""")
692    elif isinstance(self.variables, (list, tuple)):
693      self.variables = dict(zip(
694        self.tiling.tile_unit.tiles.tile_id.unique(),
695        self.variables))
696      print(f"""Only a list of variables specified, assigning to
697            available tile_ids.""")
698
699    if self.colourmaps is None:
700      self.colourmaps = {}
701      for var in self.variables.values():
702        if self.map[var].dtype == pd.CategoricalDtype:
703          self.colourmaps[var] = "tab20"
704          print(f"""For categorical data, you should specify colour
705              mapping explicitly.""")
706        else:
707          self.colourmaps[var] = "Reds"
708
709    self._plot_map(axes, **kwargs)
710    return fig
711
712
713  def _plot_map(self, axes:pyplot.Axes, **kwargs) -> None:
714    """Plots map to the supplied axes.
715
716    Args:
717      axes (pyplot.Axes): axes on which maps will be drawn.
718    """
719    bb = self.map.geometry.total_bounds
720    if self.legend:
721      axes["map"].set_axis_off()
722      axes["map"].set_xlim(bb[0], bb[2])
723      axes["map"].set_ylim(bb[1], bb[3])
724      self._plot_subsetted_gdf(axes["map"], self.map, **kwargs)
725      self.plot_legend(ax = axes["legend"], **kwargs)
726      if (self.legend_dx != 0 or self.legend_dx != 0):
727        box = axes["legend"].get_position()
728        box.x0 = box.x0 + self.legend_dx
729        box.x1 = box.x1 + self.legend_dx
730        box.y0 = box.y0 + self.legend_dy
731        box.y1 = box.y1 + self.legend_dy
732        axes["legend"].set_position(box)
733    else:
734      axes.set_axis_off()
735      axes.set_xlim(bb[0], bb[2])
736      axes.set_ylim(bb[1], bb[3])
737      self._plot_subsetted_gdf(axes, self.map, **kwargs)
738    return None
739
740
741  def _plot_subsetted_gdf(self, ax:pyplot.Axes,
742                          gdf:gpd.GeoDataFrame, **kwargs) -> None:
743    """Plots a gpd.GeoDataFrame multiple times based on a subsetting
744    attribute (assumed to be "tile_id").
745
746    NOTE: used to plot both the main map _and_ the legend.
747
748    Args:
749      ax (pyplot.Axes): axes to plot to.
750      gdf (gpd.GeoDataFrame): the GeoDataFrame to plot.
751
752    Raises:
753      Exception: if self.colourmaps cannot be parsed exception is raised.
754    """
755    groups = gdf.groupby("tile_id")
756    for id, var in self.variables.items():
757      subset = groups.get_group(id)
758      # Handle custom color assignments via 'cmaps' parameter.
759      # Result is setting 'cmap' variable used in plot command afterwards.
760      if (isinstance(self.colourmaps[var], dict)):
761        colormap_dict = self.colourmaps[var]
762        data_unique_sorted = sorted(subset[var].unique())
763        cmap = matplotlib.colors.ListedColormap(
764          [colormap_dict[x] for x in data_unique_sorted])
765        subset.plot(ax = ax, column = var, cmap = cmap, **kwargs)
766      else:
767        if (isinstance(self.colourmaps,
768                (str, matplotlib.colors.Colormap,
769                matplotlib.colors.LinearSegmentedColormap,
770                matplotlib.colors.ListedColormap))):
771          cmap = self.colourmaps   # one palette for all ids
772        elif (len(self.colourmaps) == 0):
773          cmap = 'Reds'  # set a default... here, to Brewer's 'Reds'
774        elif (var not in self.colourmaps):
775          cmap = 'Reds'  # no color specified in dict, use default
776        elif (isinstance(self.colourmaps[var],
777                (str, matplotlib.colors.Colormap,
778                matplotlib.colors.LinearSegmentedColormap,
779                matplotlib.colors.ListedColormap))):
780          cmap = self.colourmaps[var]  # specified colors for this var
781        else:
782          raise Exception(f"""Color map for '{var}' is not a known
783                          type, but is {str(type(self.colourmaps[var]))}""")
784
785        subset.plot(ax = ax, column = var, cmap = cmap,
786              scheme = self.scheme, k = self.k, **kwargs)
787
788
789  def to_file(self, fname:str = None) -> None:
790    """Outputs the tiled map to a layered GPKG file.
791
792    Currently delegates to `weavingspace.tiling_utils.write_map_to_layers()`.
793
794    Args:
795      fname (str, optional): Filename to write. Defaults to None.
796    """
797    tiling_utils.write_map_to_layers(self.map, fname)
798    return None
799
800
801  def plot_legend(self, ax: pyplot.Axes = None, **kwargs) -> None:
802    """Plots a legend for this tiled map.
803
804    Args:
805      ax (pyplot.Axes, optional): axes to draw legend. Defaults to None.
806    """
807    # turn off axes (which seems also to make it impossible
808    # to set a background colour)
809    ax.set_axis_off()
810
811    legend_tiles = self.tiling.tile_unit._get_legend_tiles()
812    # this is a bit hacky, but we will apply the rotation to text
813    # annotation so for TileUnits which don't need it, reverse that now
814    if isinstance(self.tiling.tile_unit, TileUnit):
815      legend_tiles.rotation = -self.tiling.rotation
816
817    legend_key = self._get_legend_key_gdf(legend_tiles)
818
819    legend_tiles.geometry = legend_tiles.geometry.rotate(
820      self.tiling.rotation, origin = (0, 0))
821
822    if self.use_ellipse:
823      ellipse = tiling_utils.get_bounding_ellipse(
824        legend_tiles.geometry, mag = self.ellipse_magnification)
825      bb = ellipse.total_bounds
826      c = ellipse.unary_union.centroid
827    else:
828      bb = legend_tiles.geometry.total_bounds
829      c = legend_tiles.geometry.unary_union.centroid
830
831    # apply legend zoom - NOTE that this must be applied even
832    # if self.legend_zoom is not == 1...
833    ax.set_xlim(c.x + (bb[0] - c.x) / self.legend_zoom,
834          c.x + (bb[2] - c.x) / self.legend_zoom)
835    ax.set_ylim(c.y + (bb[1] - c.y) / self.legend_zoom,
836          c.y + (bb[3] - c.y) / self.legend_zoom)
837
838    # plot the legend key tiles (which include the data)
839    self._plot_subsetted_gdf(ax, legend_key, lw = 0, **kwargs)
840
841    for id, tile, rotn in zip(self.variables.keys(),
842                              legend_tiles.geometry,
843                              legend_tiles.rotation):
844      c = tile.centroid
845      ax.annotate(self.variables[id], xy = (c.x, c.y),
846          ha = "center", va = "center", rotation_mode = "anchor",
847          # adjust rotation to favour text reading left to right
848          rotation = (rotn + self.tiling.rotation + 90) % 180 - 90,
849          bbox = {"lw": 0, "fc": "#ffffff40"})
850
851    # now plot background; we include the central tiles, since in
852    # the weave case these may not match the legend tiles
853    context_tiles = self.tiling.tile_unit.get_local_patch(r = 2,
854      include_0 = True).geometry.rotate(self.tiling.rotation, origin = (0, 0))
855    # for reasons escaping all reason... invalid polygons sometimes show up
856    # here I think because of the rotation /shrug... in any case, this
857    # sledgehammer should fix it
858    # context_tiles = gpd.GeoSeries([g.simplify(1e-6)
859    #                                for g in context_tiles.geometry],
860    #                 crs = self.tiling.tile_unit.crs)
861
862    if self.use_ellipse:
863      context_tiles.clip(ellipse, keep_geom_type = False).plot(
864        ax = ax, fc = "#9F9F9F3F", lw = 0.0)
865      tiling_utils.get_tiling_edges(context_tiles.geometry).clip(
866        ellipse, keep_geom_type = True).plot(ax = ax, ec = "#5F5F5F", lw = 1)
867    else:
868      context_tiles.plot(ax = ax, fc = "#9F9F9F3F", ec = "#5F5F5F", lw = 0.0)
869      tiling_utils.get_tiling_edges(context_tiles.geometry).plot(
870        ax = ax, ec = "#5F5F5F", lw = 1)
871
872
873  def _get_legend_key_gdf(self, tiles:gpd.GeoDataFrame) -> gpd.GeoDataFrame:
874    """Returns a GeoDataFrame of tiles dissected and with data assigned 
875    to the slice so a map of them can stand as a legend.
876
877    'Dissection' is handled differently by `WeaveUnit` and `TileUnit`
878    objects and delegated to either `WeaveUnit._get_legend_key_shapes()`
879    or `TileUnit._get_legend_key_shapes()`.
880
881    Args:
882      tiles (gpd.GeoDataFrame): the legend tiles.
883
884    Returns:
885      gpd.GeoDataFrame:  with tile_id, variables and rotation
886        attributes, and geometries of Tileable tiles sliced into a
887        colour ramp or set of nested tiles.
888    """
889    key_tiles = []   # set of tiles to form a colour key (e.g. a ramp)
890    ids = []         # tile_ids applied to the keys
891    unique_ids = []  # list of each tile_id used in order
892    vals = []        # the data assigned to the key tiles
893    rots = []        # rotation of each key tile
894    subsets = self.map.groupby("tile_id")
895    for (id, var), geom, rot in zip(self.variables.items(),
896                 tiles.geometry,
897                 tiles.rotation):
898      subset = subsets.get_group(id)
899      d = subset[var]
900      radial = False
901      # if the data are categorical then it's complicated...
902      if d.dtype == pd.CategoricalDtype:
903        radial = True and self.radial_key
904        # desired order of categorical variable is the
905        # color maps dictionary keys
906        cmap = self.colourmaps[var]
907        num_cats = len(cmap)
908        val_order = dict(zip(cmap.keys(), range(num_cats)))
909        # compile counts of each category
910        freqs = [0] * num_cats
911        for v in list(d):
912          freqs[val_order[v]] += 1
913        # make list of the categories containing appropriate
914        # counts of each in the order needed using a reverse lookup
915        data_vals = list(val_order.keys())
916        data_vals = [data_vals[i] for i, f in enumerate(freqs) if f > 0]
917      else: # any other data is easy!
918        data_vals = sorted(d)
919        freqs = [1] * len(data_vals)
920      key = self.tiling.tile_unit._get_legend_key_shapes(
921        geom, freqs, rot, radial)
922      key_tiles.extend(key)
923      vals.extend(data_vals)
924      n = len(data_vals)
925      ids.extend([id] * n)
926      unique_ids.append(id)
927      rots.extend([rot] * n)
928    # finally make up a data table with all the data in all the
929    # columns (each set of data only gets used in the subset it
930    # applies to). This allows us to reuse the tiling_utils.
931    # plot_subsetted_gdf() function
932    key_data = {}
933    for id in unique_ids:
934      key_data[self.variables[id]] = vals
935
936    key_gdf = gpd.GeoDataFrame(
937      data = key_data | {"tile_id": ids, "rotation": rots},
938      crs = self.map.crs,
939      geometry = gpd.GeoSeries(key_tiles))
940    key_gdf.geometry = key_gdf.rotate(self.tiling.rotation, origin = (0, 0))
941    return key_gdf
942
943
944  def explore(self) -> None:
945    """TODO: add wrapper to make tiled web map via geopandas.explore.
946    """
947    return None