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