weavingspace.tile_map

   1"""Classes for tiling maps.
   2
   3`weavingspace.tile_map.Tiling` and
   4`weavingspace.tile_map.TiledMap` are exposed in the  public API and
   5respectively enable creation of a tiling and plotting of the tiling as a
   6multivariate map.
   7"""
   8
   9from __future__ import annotations
  10
  11import copy
  12import itertools
  13import typing
  14from collections.abc import Iterable
  15from dataclasses import dataclass
  16from time import perf_counter
  17
  18import geopandas as gpd
  19import matplotlib.colors
  20import matplotlib.pyplot as plt
  21import numpy as np
  22import pandas as pd
  23import shapely.affinity as affine
  24import shapely.geometry as geom
  25import shapely.ops
  26
  27from weavingspace import Tileable, TileUnit, tiling_utils
  28
  29if typing.TYPE_CHECKING:
  30  from matplotlib.figure import Figure
  31
  32CMAPS_SEQUENTIAL = list(itertools.chain(*[[x, x + "_r"] for x in
  33  ["Greys", "Purples", "Blues", "Greens", "Oranges", "Reds",
  34   "viridis", "plasma", "inferno", "magma", "cividis",
  35   "YlOrBr", "YlOrRd", "OrRd", "PuRd", "RdPu", "BuPu",
  36   "GnBu", "PuBu", "YlGnBu", "PuBuGn", "BuGn", "YlGn",
  37   "binary", "gist_yarg", "gist_gray", "gray", "bone",
  38   "pink", "spring", "summer", "autumn", "winter", "cool",
  39   "Wistia", "hot", "afmhot", "gist_heat", "copper"]]))
  40
  41CMAPS_DIVERGING = list(itertools.chain(*[[x, x + "_r"] for x in
  42  ["PiYG", "PRGn", "BrBG", "PuOr", "RdGy", "RdBu", "RdYlBu",
  43   "RdYlGn", "Spectral", "coolwarm", "bwr", "seismic"]]))
  44# these ones not yet in MPL 3.8.4: "berlin", "managua", "vanimo"
  45
  46CMAPS_CATEGORICAL = list(itertools.chain(*[[x, x + "_r"] for x in
  47  ["Pastel1", "Pastel2", "Paired", "Accent", "Dark2",
  48   "Set1", "Set2", "Set3", "tab10", "tab20", "tab20b",
  49   "tab20c"]]))
  50
  51@dataclass(slots=True, init=False)
  52class _TileGrid:
  53  """A class to represent the translation centres of a tiling.
  54
  55  We store the grid as a GeoSeries of Point objects to make it simple to plot
  56  in map views if required.
  57
  58  Implementation relies on transforming the translation vectors into a square
  59  space where the tile spacing is unit squares, then transforming this back
  60  into the original map space. Some member variables of the class are in
  61  the transformed grid generation space, some in the map space. We store the
  62  translation vectors as shapely Points, and the extent of the grid in grid
  63  space as a shapely Polygon so we can more easily visualise if needed.
  64  """
  65
  66  tile_unit:Tileable
  67  """the base tile in map space."""
  68  oriented_rect_to_tile:geom.Polygon
  69  """oriented rectangular region in map space that is to be tiled."""
  70  to_grid_space:tuple[float,...]
  71  """the forward transformation from map space to tiling grid generation space,
  72  stored as a shapely.affinity transform tuple of 6 floats."""
  73  to_map_space:tuple[float,...]
  74  """the inverse transform from tiling grid space to map space."""
  75  extent_in_grid_space:geom.Polygon
  76  """geometry of the circular extent of the tiling transformed into grid
  77  generation space."""
  78  centre:geom.Point
  79  """centre point of the grid in map space - this is required for later
  80  rotations of a generated tiling in map space"""
  81  points:gpd.GeoSeries
  82  """geom.Points recording translation vectors of the tiling in map space."""
  83
  84  def __init__(
  85      self,
  86      tile_unit:Tileable,
  87      to_tile:gpd.GeoSeries,
  88      at_centroids:bool = False) -> None:
  89    self.tile_unit = tile_unit
  90    self.oriented_rect_to_tile = self._get_rect_to_tile(to_tile)
  91    self.to_map_space, self.to_grid_space = self._get_transforms()
  92    self._set_centre_in_map_space()
  93    self._set_extent_in_grid_space()
  94    if at_centroids:
  95      self.points = to_tile.representative_point()
  96    else:
  97      self.points = self._get_grid()
  98    self.points.crs = self.tile_unit.crs
  99
 100
 101  def _get_rect_to_tile(
 102      self,
 103      region_to_tile:gpd.GeoSeries) -> geom.Polygon:
 104    """Generate an oriented rectangle that encompasses the region to be tiled.
 105
 106    This will be rotated to generate a circular are that will actually be tiled
 107    so that if new orientations of the tiled pattern are requested they can be
 108    generated quickly without much recalculation.
 109
 110    Args:
 111      region_to_tile (gpd.GeoSeries): the region to be tiled.
 112
 113    Returns:
 114      geom.Polygon: an oriented rectangle encompassing the region to be tiled
 115        with some buffering to avoid cutting off tiling elements.
 116
 117    """
 118    # buffer the region by an amount dictated by the size of the tile unit
 119    bb = self.tile_unit.tiles.total_bounds
 120    diagonal = np.hypot(bb[2] - bb[0], bb[3] - bb[1])
 121    return (region_to_tile.union_all().buffer(diagonal)
 122            .minimum_rotated_rectangle)
 123
 124
 125  def _get_transforms(self) -> tuple[tuple[float,...],tuple[float,...]]:
 126    """Return forward and inverse transforms from map space to grid space.
 127
 128    In grid generation space the translation vectors are (1, 0) and (0, 1)
 129    so we simply form a matrix from two translation vectors to get the
 130    transform from grid space to map space, and invert it to get transform from
 131    map space to grid space.
 132
 133    See pages 18-22 in Kaplan CS, 2009. Introductory Tiling Theory for Computer
 134    Graphics (Morgan & Claypool) for the logic of this approach.
 135
 136    The results are returned as shapely affine transform tuples.
 137
 138    Returns:
 139      tuple[tuple[float,...],tuple[float,...]]: shapely affine transform tuples
 140        (a, b, d, e, dx, dy).
 141
 142    """
 143    v = self.tile_unit.get_vectors()
 144    # unpack the first two vectors and funnel them into a 2x2 matrix
 145    grid_to_map = np.array([[v[0][0], v[1][0]], [v[0][1], v[1][1]]])
 146    map_to_grid = np.linalg.inv(grid_to_map)
 147    return (self._np_to_shapely_transform(grid_to_map),
 148            self._np_to_shapely_transform(map_to_grid))
 149
 150
 151  def _set_centre_in_map_space(self) -> None:
 152    """Set the tiling centre in map space."""
 153    self.centre = self.oriented_rect_to_tile.centroid
 154
 155
 156  def _set_extent_in_grid_space(self) -> None:
 157    """Set extent of the grid in grid generation space."""
 158    corner = geom.Point(self.oriented_rect_to_tile.exterior.coords[0])
 159    radius = self.centre.distance(corner)
 160    self.extent_in_grid_space = \
 161      affine.affine_transform(self.centre.buffer(radius), self.to_grid_space)
 162
 163
 164  def _get_grid(self) -> gpd.GeoSeries:
 165    """Generate the grid transformed into map space.
 166
 167    Obtain dimensions of the transformed region, then set down a uniform grid.
 168    Grid generation is greatly accelerated by using the numpy meshgrid method
 169    which we can do because we're working in a grid-generation space where
 170    tile units have been transformed to the unit square.
 171
 172    Returns:
 173      gpd.GeoSeries: the grid as a collection of geom.Points.
 174
 175    """
 176    bb = self.extent_in_grid_space.bounds
 177    _w, _h, _l, _b = bb[2] - bb[0], bb[3] - bb[1], bb[0], bb[1]
 178    w = int(np.ceil(_w))
 179    h = int(np.ceil(_h))
 180    l = _l - (w - _w) / 2
 181    b = _b - (h - _h) / 2
 182    xs, ys = np.array(np.meshgrid(np.arange(w) + l,
 183                                  np.arange(h) + b)).reshape((2, w * h))
 184    pts = [geom.Point(x, y) for x, y in zip(xs, ys, strict = True)]
 185    return (gpd.GeoSeries(
 186      [p for p in pts if p.within(self.extent_in_grid_space)])
 187        .affine_transform(self.to_map_space))
 188
 189
 190  def _np_to_shapely_transform(
 191      self,
 192      array:np.ndarray) -> tuple[float,...]:
 193    """Convert a numpy affine transform matrix to shapely format.
 194
 195        [[a b c]
 196         [d e f] --> (a b d e c f)
 197         [g h i]]
 198
 199    There is no translation in transforms in this use-case so c == f == 0
 200
 201    Args:
 202      array (np.ndarray): numpy affine transform array to convert.
 203
 204    Returns:
 205      tuple[float,...]: shapely affine transform tuple
 206
 207    """
 208    return (*list(array[:2, :2].flatten()), 0, 0)
 209
 210
 211@dataclass(slots=True, init=False)
 212class Tiling:
 213  """Class that applies a `Tileable` object to a region to be mapped.
 214
 215  The result of the tiling procedure is stored in the `tiles` variable and
 216  covers a region sufficient that the tiling can be rotated to any desired
 217  angle. Rotation can be requested when the render method is called.
 218  """
 219
 220  tileable:Tileable
 221  """Tileable on which the tiling is based."""
 222  region:gpd.GeoDataFrame
 223  """the region to be tiled."""
 224  region_union:geom.Polygon
 225  """a single polygon of all the areas in the region to be tiled"""
 226  grid:_TileGrid
 227  """the grid which will be used to apply the tiling."""
 228  tiles:gpd.GeoDataFrame
 229  """the tiles after tiling has been carried out."""
 230  prototiles:gpd.GeoDataFrame
 231  """the prototiles after tiling has been carried out."""
 232  rotation:float
 233  """additional rotation applied to the tiling beyond any that might have
 234  been 'baked in' to the Tileable."""
 235
 236  def __init__(
 237      self,
 238      tileable:Tileable,
 239      region:gpd.GeoDataFrame,
 240      as_icons:bool = False,
 241    ) -> None:
 242    """Construct a tiling by polygons extending beyond supplied region.
 243
 244    The tiling is extended sufficiently to allow for its application at any
 245    rotation.
 246
 247    Args:
 248      tileable (Tileable): the TileUnit or WeaveUnit to use.
 249      region (gpd.GeoDataFrame): the region to be tiled.
 250      as_icons (bool, optional): if True prototiles will only be placed at the
 251        region's zone centroids, one per zone. Defaults to False.
 252
 253    """
 254    self.tileable = tileable
 255    self.rotation = 0
 256    self.region = region
 257    self.region.sindex # this probably speeds up overlay
 258    self.region_union = self.region.geometry.union_all()
 259    self.grid = _TileGrid(
 260      self.tileable,
 261      self.region.geometry if as_icons else gpd.GeoSeries([self.region_union]),
 262      as_icons)
 263    self.tiles, self.prototiles = self.make_tiling()
 264    self.tiles.sindex # again this probably speeds up overlay
 265
 266
 267  def get_tiled_map(
 268      self,
 269      rotation:float = 0.,
 270      join_on_prototiles:bool = False,
 271      retain_tileables:bool = False,
 272      prioritise_tiles:bool = True,
 273      ragged_edges:bool = True,
 274      use_centroid_lookup_approximation:bool = False,
 275      debug:bool = False,
 276    ) -> TiledMap:
 277    """Return a `TiledMap` filling a region at the requested rotation.
 278
 279    HERE BE DRAGONS! This function took a lot of trial and error to get right,
 280    so modify with CAUTION!
 281
 282    The `proritise_tiles = True` option means that the tiling will not break
 283    up the tiles in `TileUnit`s at the boundaries between areas in the mapped
 284    region, but will instead ensure that tiles remain complete, picking up
 285    their data from the region zone which they overlap the most.
 286
 287    The exact order in which operations are performed affects performance. For
 288    example, the final clipping to self.region when ragged_edges = False is
 289    _much_ slower if it is carried out before the dissolving of tiles into the
 290    region zones. So... again... modify CAREFULLY!
 291
 292    Args:
 293      rotation (float, optional): Optional rotation to apply. Defaults to 0.
 294      join_on_prototiles (bool, optional): if True data from the region dataset
 295        are joined to tiles based on the prototile to which they belong. If
 296        False the join is based on the tiles in relation to the region areas.
 297        For weave-based tilings False is probably to be preferred. Defaults to
 298        False.
 299      retain_tileables (bool, optional): if True complete tileable units will
 300        be retained. If False tile unit elements that do not overlap the map
 301        area will be discarded.
 302      prioritise_tiles (bool, optional): if True tiles will not be broken at
 303        boundaries in the region dataset. Defaults to True.
 304      ragged_edges (bool, optional): if True tiles at the edge of the region
 305        will not be cut by the region extent - ignored if prioritise_tiles is
 306        False when edges will always be clipped to the region extent. Defaults
 307        to True.
 308      use_centroid_lookup_approximation (bool, optional): if True use tile
 309        centroids for lookup of region data - ignored if prioritise_tiles is
 310        False when it is irrelevant. Defaults to False.
 311      debug (bool, optional): if True prints timing messages. Defaults to
 312        False.
 313
 314    Returns:
 315      TiledMap: a TiledMap of the region with attributes attached to tiles.
 316
 317    """
 318    if debug:
 319      t1 = perf_counter()
 320
 321    id_var = self._setup_region_DZID()
 322    if join_on_prototiles:
 323      if rotation == 0:
 324        tiled_map, join_layer = self.tiles, self.prototiles
 325      else:
 326        tiled_map, join_layer = self.rotated(rotation)
 327      tiled_map["joinUID"] = self.tiles["prototile_id"]
 328    else:
 329      tiled_map = self.tiles if rotation == 0 else self.rotated(rotation)[0]
 330      tiled_map["joinUID"] = self.tiles["tile_id"]
 331      join_layer = tiled_map
 332    join_layer["joinUID"] = list(range(join_layer.shape[0]))
 333
 334    # compile a list of the variable names we are NOT going to change
 335    # i.e. everything except the geometry and the id_var
 336    region_vars = [column for column in self.region.columns
 337                   if "geom" not in column and column != id_var]
 338
 339    if debug:
 340      t2 = perf_counter()
 341      print(f"STEP 1: prep data (rotation if requested): {t2 - t1:.3f}")
 342
 343    if prioritise_tiles:
 344      # maintain tile continuity across zone boundaries
 345      # so we have to do more work than a simple overlay
 346      if use_centroid_lookup_approximation:
 347        t5 = perf_counter()
 348        tile_pts = copy.deepcopy(join_layer)
 349        tile_pts.geometry = tile_pts.centroid
 350        lookup = tile_pts.sjoin(
 351          self.region, how = "inner")[["joinUID", id_var]]
 352      else:
 353        # determine areas of overlapping tiles and drop the data we join the
 354        # data back later, so dropping makes that easier overlaying in region.
 355        # overlay(tiles) seems to be faster??
 356        # TODO: also... this part is performance-critical, think about fixes --
 357        # possibly including the above centroid-based approximation
 358        overlaps = self.region.overlay(join_layer, make_valid = False)
 359        if debug:
 360          t3 = perf_counter()
 361          print(f"STEP A2: overlay zones with tiling: {t3 - t2:.3f}")
 362        overlaps["area"] = overlaps.geometry.area
 363        if debug:
 364          t4 = perf_counter()
 365          print(f"STEP A3: calculate areas: {t4 - t3:.3f}")
 366        overlaps = overlaps.drop(columns = region_vars)
 367        if debug:
 368          t5 = perf_counter()
 369          print(f"STEP A4: drop columns prior to join: {t5 - t4:.3f}")
 370        # make a lookup by largest area tile to region id
 371        lookup = overlaps \
 372          .iloc[overlaps.groupby("joinUID")["area"] \
 373          .agg(pd.Series.idxmax)][["joinUID", id_var]]
 374      # now join the lookup and from there the region data
 375      if debug:
 376        t6 = perf_counter()
 377        print(f"STEP A5: build lookup for join: {t6 - t5:.3f}")
 378      tiled_map = tiled_map \
 379        .merge(lookup, on = "joinUID") \
 380        .merge(self.region.drop(columns = ["geometry"]), on = id_var)
 381      if debug:
 382        t7 = perf_counter()
 383        print(f"STEP A6: perform lookup join: {t7 - t6:.3f}")
 384      tiled_map = tiled_map.drop(columns = ["joinUID"])
 385
 386    else:  
 387      # here it's a simple overlay
 388      tiled_map = self.region.overlay(tiled_map)
 389      t7 = perf_counter()
 390      if debug:
 391        print(f"STEP B2: overlay tiling with zones: {t7 - t2:.3f}")
 392
 393    if not retain_tileables:
 394      tiled_map = tiled_map.loc[
 395        shapely.intersects(self.region_union, np.array(tiled_map.geometry)), :]
 396    
 397    # inplace changes considered unsafe, BUT not dropping id_var in this
 398    # causes it to persist in the tiled_map and region dataframes!
 399    tiled_map.drop(columns = [id_var], inplace = True)
 400    self.region.drop(columns = [id_var], inplace = True)
 401
 402    # if we've retained tiles and want 'clean' edges, then clip
 403    # note that this step is slow: geopandas unary_unions the clip layer
 404    if prioritise_tiles and not ragged_edges:
 405      tiled_map.sindex
 406      tiled_map = tiled_map.clip(self.region_union)
 407      if debug:
 408        print(f"""STEP A7/B3: clip map to region: {perf_counter() - t7:.3f}""")
 409
 410    tm = TiledMap()
 411    tm.tiling = self
 412    tm.map = tiled_map
 413    return tm
 414
 415
 416  def _setup_region_DZID(self) -> str:
 417    """Return guaranteed-unique new attribute name for self.region dataframe.
 418
 419    Avoids a name clash with any existing attribute in the dataframe.
 420
 421    Returns:
 422      str: name of the added attribute.
 423
 424    """
 425    dzid = "DZID"
 426    i = 0
 427    while dzid in self.region.columns:
 428      dzid = "DZID" + str(i)
 429      i = i + 1
 430    self.region[dzid] = list(range(self.region.shape[0]))
 431    return dzid
 432
 433
 434  def make_tiling(self) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
 435    """Tile the region with a tile unit, returning a GeoDataFrame.
 436
 437    Returns:
 438      geopandas.GeoDataFrame: a GeoDataFrame of the region tiled with the
 439        tile unit.
 440
 441    """
 442    # we assume the geometry column is called geometry so make it so...
 443    if self.region.geometry.name != "geometry":
 444      self.region = self.region.rename_geometry("geometry")
 445
 446    # chain list of lists of GeoSeries geometries to list of geometries
 447    tiles = itertools.chain(*[
 448      self.tileable.tiles.geometry.translate(p.x, p.y)
 449      for p in self.grid.points])
 450    prototiles = itertools.chain(*[
 451      self.tileable.prototile.geometry.translate(p.x, p.y)
 452      for p in self.grid.points])
 453    # replicate the tile ids
 454    tile_ids = list(self.tileable.tiles.tile_id) * len(self.grid.points)
 455    prototile_ids = list(range(len(self.grid.points)))
 456    tile_prototile_ids = sorted(prototile_ids * self.tileable.tiles.shape[0])
 457    tiles_gs = gpd.GeoSeries(list(tiles))
 458    prototiles_gs = gpd.GeoSeries(list(prototiles))
 459    # assemble and return as GeoDataFrames
 460    tiles_gdf = gpd.GeoDataFrame(
 461      data = {"tile_id": tile_ids, "prototile_id": tile_prototile_ids},
 462      geometry = tiles_gs, crs = self.tileable.crs)
 463    prototiles_gdf = gpd.GeoDataFrame(
 464      data = {"prototile_id": prototile_ids},
 465      geometry = prototiles_gs, crs = self.tileable.crs)
 466    return tiles_gdf, prototiles_gdf
 467
 468
 469  def rotated(self,
 470              rotation:float = 0.0,
 471              ) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
 472    """Return the stored tiling rotated.
 473
 474    The stored tiling never changes and
 475    if it was originally made with a Tileable that was rotated it will retain
 476    that rotation. The requested rotation is _additional_ to that baseline
 477    rotation.
 478
 479    Args:
 480      rotation (float, optional): Rotation angle in degrees.
 481        Defaults to None.
 482
 483    Returns:
 484      gpd.GeoDataFrame: Rotated tiling.
 485
 486    """
 487    if self.tiles is None:
 488      self.tiles = self.make_tiling()[0]
 489    if rotation == 0:
 490      return self.tiles, self.prototiles
 491    tiles = gpd.GeoDataFrame(
 492      data = {"tile_id": self.tiles.tile_id,
 493              "prototile_id": self.tiles.tile_id},
 494      crs = self.tiles.crs,
 495      geometry = self.tiles.geometry.rotate(
 496        rotation, origin = self.grid.centre))
 497    prototiles = gpd.GeoDataFrame(
 498      data = {"prototile_id": self.prototiles.prototile_id},
 499      crs = self.prototiles.crs,
 500      geometry = self.prototiles.geometry.rotate(
 501        rotation, origin = self.grid.centre))
 502    self.rotation = rotation
 503    return tiles, prototiles
 504
 505
 506@dataclass(slots=True)
 507class TiledMap:
 508  """Class representing a tiled map.
 509
 510  Should not be accessed directly, but will
 511  be created by calling `Tiling.get_tiled_map()`. After creation the variables
 512  and colourmaps attributes can be set, and then `TiledMap.render()` called to
 513  make a map. Settable attributes are explained in documentation of the
 514  `TiledMap.render()` method.
 515
 516  Examples:
 517    Recommended usage is as follows. First, make a `TiledMap` from a `Tiling`
 518    object:
 519
 520      `tm = tiling.get_tiled_map(...)`
 521
 522    Some options in the `Tiling` constructor affect the map appearance. See
 523    `Tiling` for details.
 524
 525    Once a `TiledMap` object exists, set options on it, either when calling
 526    `TiledMap.render()` or explicitly, i.e.
 527
 528      tm.render(opt1 = val1, opt2 = val2, ...)
 529
 530    or
 531
 532      tm.opt1 = val1
 533      tm.opt2 = val2
 534      tm.render()
 535
 536    Option settings are persistent, i.e. unless a new `TiledMap` object is
 537    created the option settings have to be explicitly reset to new values on
 538    subsequent calls to `TiledMap.render()`.
 539
 540    The most important options are the `vars_map` and `colors_to_use` settings.
 541
 542    `vars_to_map` is a lost of the dataset variable names to match with
 543    `weavingspace.tileable.Tileable` elements with corresponding (ordered)
 544    tile_ids (usually "a", "b", etc.). If you need to control the match, then
 545    also supply `ids_to_map` in matching order. E.g.
 546
 547      tm.ids_to_map = ['d', 'c', 'b', 'a']
 548      tm.vars_to_map = ['x1', 'x2', 'x3', 'x4']
 549
 550    Note that this means that if you really want more than one element in the
 551    tiling to represent the same variable more than once, you can do that.
 552
 553    `colors_to_use` is a parallel list of named matplotlib colormaps,
 554
 555      tm.colors_to_use = ["Reds", "Blues", "Greys", "Purples"]
 556
 557    Similarly, you can specify the classification `schemes_to_use` (such as
 558    'quantiles') and the number of classes `n_classes` in each.
 559
 560    If data are categorical, this is flagged in the `categoricals` list of
 561    booleans, in which case an appropriate colour map should be used. There is
 562    currently no provision for control of which colour in a categorical
 563    colour map is applied to which variable level.
 564
 565    TODO: better control of categorical mapping schemes.
 566
 567  """
 568
 569  # these will be set at instantion by Tiling.get_tiled_map()
 570  tiling:Tiling = None
 571  """the Tiling with the required tiles"""
 572  map:gpd.GeoDataFrame = None
 573  """the GeoDataFrame on which this map is based"""
 574  ids_to_map:list[str] = None
 575  """tile_ids that are to be used to represent data"""
 576  vars_to_map:list[str] = None
 577  """dataset variables that are to be symbolised"""
 578  colors_to_use:list[str|list[str]] = None
 579  """list of matplotlib colormap names."""
 580  categoricals:list[bool] = None
 581  """list specifying if each variable is -- or is to be treated as --
 582  categorical"""
 583  schemes_to_use:list[str|None] = None
 584  """mapclassify schemes to use for each variable."""
 585  n_classes:list[int|None] = None
 586  """number of classes to apply; if set to 0 will be unclassed."""
 587  _colourspecs:dict[str,dict] = None
 588  """dictionary of dictionaries keyed by the items in `ids_to_use` with each
 589  dictionary forming additional kwargs to be supplied to geopandas.plot()."""
 590  add_buffer:bool = False
 591  """if True then include a buffer of all tiles as a background"""
 592  buffer_colour:str = "grey"
 593  """colour of any include buffer layer"""
 594
 595  # the below parameters can be set either before calling self.render() or
 596  # passed in as parameters to self.render(). These are solely
 597  # `TiledMap.render()` options not geopandas plot options.
 598  legend:bool = True
 599  """whether or not to show a legend"""
 600  legend_zoom:float = 1.0
 601  """<1 zooms out from legend to show more context"""
 602  legend_dx:float = 0.
 603  """x shift of legend relative to the map"""
 604  legend_dy:float = 0.
 605  """y shift of legend relative to the map"""
 606  use_ellipse:bool = False
 607  """if True clips legend with an ellipse"""
 608  ellipse_magnification:float = 1.0
 609  """magnification to apply to clip ellipse"""
 610  radial_key:bool = False
 611  """if True use radial key even for ordinal/ratio data (normally these will be
 612  shown by concentric tile geometries)"""
 613  draft_mode:bool = False
 614  """if True plot the map coloured by tile_id"""
 615
 616  # the parameters below are geopandas.plot options which we intercept to
 617  # ensure they are applied appropriately when we plot a GDF
 618  figsize:tuple[float,float] = (20, 15)
 619  """maptlotlib figsize"""
 620  dpi:float = 72
 621  """dpi for bitmap formats"""
 622
 623  def render(
 624      self,
 625      **kwargs,
 626    ) -> Figure:
 627    """Render the current state to a map.
 628
 629    Note that TiledMap objects will usually be created by calling
 630    `Tiling.get_tiled_map()`.
 631
 632    Args:
 633      ids_to_map (list[str]): tile_ids to be used in the map. Defaults to None.
 634      vars_to_map (list[str]): dataset columns to be mapped. Defaults to None.
 635      colors_to_use (list[str]): list of matplotlib colormaps to be used.
 636        Defaults to None.
 637      categoricals (list[bool]): list of flags indicating if associated variable
 638        should be treated as categorical. Defaults to None.
 639      schemes_to_use (list[str]): list of strings indicating the mapclassify
 640        scheme to use e.g. 'equalinterval' or 'quantiles'. Defaults to None.
 641      n_classes (list[int]): list of ints indicating number of classes to use in
 642        map classification. Defaults to None.
 643      legend (bool): If True a legend will be drawn. Defaults to True.
 644      legend_zoom (float): Zoom factor to apply to the legend. Values <1
 645        will show more of the tile context. Defaults to 1.0.
 646      legend_dx (float): x shift to apply to the legend position in plot area
 647        relative units, i.e. 1.0 is full width of plot. Defaults to 0.0.
 648      legend_dy (float): x and y shift to apply to the legend position in plot
 649        area relative units, i.e. 1.0 is full height of plot. Defaults to 0.0.
 650      use_ellipse (bool): If True applies an elliptical clip to the legend.
 651        Defaults to False.
 652      ellipse_magnification (float): Magnification to apply to ellipse clipped
 653        legend. Defaults to 1.0.
 654      radial_key (bool): If True legend key for TileUnit maps will be based on
 655        radially dissecting the tiles, i.e. pie slices. Defaults to False.
 656      draft_mode (bool): If True a map of the tiled map coloured by tile_ids
 657        (and with no legend) is returned. Defaults to False.
 658      figsize (tuple[float,floar]): plot dimensions passed to geopandas.
 659        plot. Defaults to (20, 15).
 660      dpi (float): passed to pyplot.plot. Defaults to 72.
 661      **kwargs: other settings to pass to pyplot/geopandas.plot.
 662
 663    Returns:
 664      matplotlib.figure.Figure: figure on which map is plotted.
 665
 666    """
 667    plt.rcParams["pdf.fonttype"] = 42
 668    plt.rcParams["pdf.use14corefonts"] = True
 669    matplotlib.rcParams["pdf.fonttype"] = 42
 670
 671    to_remove = set()  # keep track of kwargs we use to setup TiledMap
 672    # kwargs with no corresponding class attribute will be discarded
 673    # because we are using slots, we have to use setattr() here
 674    for k, v in kwargs.items():
 675      if k in self.__slots__:
 676        setattr(self, k, v)
 677        to_remove.add(k)
 678    # remove any them so we don't pass them on to pyplot and get errors
 679    for k in to_remove:
 680      del kwargs[k]
 681
 682    if self.draft_mode:
 683      fig = plt.figure(figsize = self.figsize)
 684      ax = fig.add_subplot(111)
 685      self.map.plot(ax = ax, column = "tile_id", cmap = "tab20", **kwargs)
 686      ax.set_axis_off()
 687      return fig
 688
 689    if self.legend:
 690      # this sizing stuff is rough and ready for now, possibly forever...
 691      reg_w, reg_h, *_ = \
 692        tiling_utils.get_width_height_left_bottom(self.map.geometry)
 693      tile_w, tile_h, *_ = \
 694        tiling_utils.get_width_height_left_bottom(
 695          self.tiling.tileable._get_legend_tiles().rotate(
 696            self.tiling.rotation, origin = (0, 0)))
 697      sf_w, sf_h = reg_w / tile_w / 3, reg_h / tile_h / 3
 698      gskw = {"height_ratios": [sf_h * tile_h, reg_h - sf_h * tile_h],
 699              "width_ratios":  [reg_w, sf_w * tile_w]}
 700
 701      fig, axes = plt.subplot_mosaic([["map", "legend"], ["map", "."]],
 702                                        gridspec_kw = gskw,
 703                                        figsize = self.figsize,
 704                                        layout = "constrained", **kwargs)
 705    else:
 706      fig, axes = plt.subplots(1, 1, figsize = self.figsize,
 707                               layout = "constrained", **kwargs)
 708
 709    self._plot_map(axes, **kwargs)
 710    return fig
 711
 712
 713  def _plot_map(
 714      self,
 715      ax:plt.Axes,
 716      **kwargs,
 717    ) -> None:
 718    """Plot map to the supplied axes.
 719
 720    Args:
 721      axes (plt.Axes): axes on which maps will be drawn.
 722      kwargs (dict): additional parameters to be passed to plot.
 723
 724    """
 725    self._set_colourspecs()
 726    bb = self.map.geometry.total_bounds
 727    if self.legend:
 728      if self.add_buffer:
 729        self._plot_buffer(ax["map"], **kwargs)
 730      if (self.legend_dx != 0 or self.legend_dx != 0):
 731        box = ax["legend"].get_position()
 732        box.x0 += self.legend_dx
 733        box.x1 += self.legend_dx
 734        box.y0 += self.legend_dy
 735        box.y1 += self.legend_dy
 736        ax["legend"].set_position(box)
 737      ax["map"].set_axis_off()
 738      ax["map"].set_xlim(bb[0], bb[2])
 739      ax["map"].set_ylim(bb[1], bb[3])
 740      self._plot_subsetted_gdf(ax["map"], self.map, **kwargs)
 741      self.plot_legend(ax = ax["legend"], **kwargs)
 742    else:
 743      if self.add_buffer:
 744        self._plot_buffer(ax, **kwargs)
 745      ax.set_axis_off()
 746      ax.set_xlim(bb[0], bb[2])
 747      ax.set_ylim(bb[1], bb[3])
 748      self._plot_subsetted_gdf(ax, self.map, **kwargs)
 749
 750
 751  def _plot_subsetted_gdf(
 752      self,
 753      ax:plt.Axes,
 754      gdf:gpd.GeoDataFrame,
 755      grouping_var:str = "tile_id",
 756      **kwargs,
 757    ) -> None:
 758    """Repeatedly plot a gpd.GeoDataFrame based on a subsetting attribute.
 759
 760    NOTE: used to plot both the main map _and_ the legend, which is why
 761    a separate GeoDataframe is supplied and we don't just use self.map.
 762
 763    Args:
 764      ax (pyplot.Axes): axes to plot to.
 765      gdf (gpd.GeoDataFrame): the GeoDataFrame to plot.
 766
 767    Raises:
 768      Exception: if self.colourmaps cannot be parsed exception is raised.
 769
 770    """
 771    groups = gdf.groupby(grouping_var)
 772    for ID, cspec in self._colourspecs.items():
 773      subset = groups.get_group(ID)
 774      n_values = len(subset[cspec["column"]].unique())
 775      if not cspec["categorical"] and n_values == 1:
 776        print(f"""
 777              Only one level in variable {cspec['column']}, replacing requested
 778              colour map with single colour fill.
 779              """)
 780        cspec["color"] = matplotlib.colormaps.get(cspec["cmap"])(0.5)
 781        del cspec["column"]
 782        del cspec["cmap"]
 783        del cspec["scheme"]
 784      elif not cspec["categorical"] and n_values < cspec["k"]:
 785        cspec["k"] = n_values
 786      elif cspec["categorical"]:
 787        del cspec["scheme"]
 788      subset.plot(ax = ax, **cspec, **kwargs)
 789
 790
 791  def to_file(self, fname:str) -> None:
 792    """Output the tiled map to a layered GPKG file.
 793
 794    Currently delegates to `weavingspace.tiling_utils.write_map_to_layers()`.
 795
 796    Args:
 797      fname (str): Filename to write. Defaults to None.
 798
 799    """
 800    tiling_utils.write_map_to_layers(self.map, fname)
 801
 802
 803  def _plot_buffer(self, ax) -> None:
 804    buffer = self.map.geometry \
 805      .buffer(10, cap_style = "square", join_style = "mitre", resolution = 1) \
 806      .union_all()
 807    gdf = gpd.GeoDataFrame(
 808      geometry = gpd.GeoSeries([buffer]), crs = self.map.crs)
 809    gdf.plot(ax = ax, fc = self.buffer_colour)
 810
 811
 812  def plot_legend(self, ax, **kwargs) -> None:
 813    """Plot a legend for this tiled map.
 814
 815    Args:
 816      ax (pyplot.Axes, optional): axes to draw legend. Defaults to None.
 817
 818    """
 819    ax.set_axis_off()
 820    legend_tiles = self.tiling.tileable._get_legend_tiles()
 821    # this is a bit hacky, but we will apply the rotation to text
 822    # annotation so for TileUnits which don't need it, reverse that now
 823    if isinstance(self.tiling.tileable, TileUnit):
 824      # note that this confuses type hinting because pandas silently assigns 
 825      # a scalar value to a Series
 826      legend_tiles.rotation = -self.tiling.rotation
 827
 828    legend_key = self._get_legend_key_gdf(legend_tiles)
 829    legend_tiles.geometry = legend_tiles.geometry.rotate(
 830      self.tiling.rotation, origin = (0, 0))
 831
 832    if self.use_ellipse:
 833      ellipse = tiling_utils.get_bounding_ellipse(
 834        legend_tiles.geometry, mag = self.ellipse_magnification)
 835      bb = ellipse.total_bounds
 836      c = ellipse.union_all().centroid
 837    else:
 838      bb = legend_tiles.geometry.total_bounds
 839      c = legend_tiles.geometry.union_all().centroid
 840
 841    # apply legend zoom - NOTE that this must be applied even
 842    # if self.legend_zoom is == 1...
 843    ax.set_xlim(c.x + (bb[0] - c.x) / self.legend_zoom,
 844                c.x + (bb[2] - c.x) / self.legend_zoom)
 845    ax.set_ylim(c.y + (bb[1] - c.y) / self.legend_zoom,
 846                c.y + (bb[3] - c.y) / self.legend_zoom)
 847
 848    for cs, tile, rotn in zip(self._colourspecs.values(),
 849                              legend_tiles.geometry,
 850                              legend_tiles.rotation):
 851      c = tile.centroid
 852      ax.annotate(cs["column"], xy = (c.x, c.y),
 853                  ha = "center", va = "center",
 854                  rotation_mode = "anchor",
 855                  # adjust rotation to favour text reading left to right
 856                  rotation = (rotn + self.tiling.rotation + 90) % 180 - 90,
 857                  bbox = {"lw": 0, "fc": "#ffffff60"})
 858
 859    # now plot background; we include the central tiles, since in
 860    # the weave case these may not match the legend tiles
 861    context_tiles = self.tiling.tileable \
 862      .get_local_patch(r = 2, include_0 = True) \
 863      .geometry.rotate(self.tiling.rotation, origin = (0, 0))
 864    if self.use_ellipse:
 865      context_tiles.clip(ellipse, keep_geom_type = False).plot(
 866        ax = ax, fc = "#9F9F9F3F", lw = .35)
 867      tiling_utils.get_tiling_edges(context_tiles.geometry).clip(
 868        ellipse, keep_geom_type = True).plot(
 869          ax = ax, ec = "#5F5F5F", lw = .35)
 870    else:
 871      context_tiles.plot(ax = ax, fc = "#9F9F9F3F",
 872                         ec = "#5F5F5F", lw = .35)
 873      tiling_utils.get_tiling_edges(context_tiles.geometry).plot(
 874        ax = ax, ec = "#5F5F5F", lw = .35)
 875
 876    # plot the legend key tiles (which include the data)
 877    self._plot_subsetted_gdf(ax, legend_key, **kwargs)
 878
 879
 880  def _get_legend_key_gdf(self, tiles:gpd.GeoDataFrame) -> gpd.GeoDataFrame:
 881    """Return tiles dissected and with data assigned for use as a legend.
 882
 883    'Dissection' is handled differently by `WeaveUnit` and `TileUnit`
 884    objects and delegated to either `WeaveUnit._get_legend_key_shapes()`
 885    or `TileUnit._get_legend_key_shapes()`.
 886
 887    Args:
 888      tiles (gpd.GeoDataFrame): the legend tiles.
 889
 890    Returns:
 891      gpd.GeoDataFrame:  with tile_id, variables and rotation
 892        attributes, and geometries of Tileable tiles sliced into a
 893        colour ramp or set of nested tiles.
 894
 895    """
 896    key_tiles = []   # set of tiles to form a colour key (e.g. a ramp)
 897    ids = []         # tile_ids applied to the keys
 898    unique_ids = []  # list of each tile_id used in order
 899    vals = []        # the data assigned to the key tiles
 900    rots = []        # rotation of each key tile
 901    # subsets = self.map.groupby("tile_id")
 902    for (id, cspec), geom, rot in zip(self._colourspecs.items(),
 903                                      tiles.geometry, 
 904                                      tiles.rotation):
 905      d = list(self.map.loc[self.map.tile_id == id][cspec["column"]])
 906      # if the data are categorical then it's complicated...
 907      # if cs["categorical"]:
 908      #   radial = True and self.radial_key
 909      #   # desired order of categorical variable is the
 910      #   # color maps dictionary keys
 911      #   num_cats = len(cmap)
 912      #   val_order = dict(zip(cmap.keys(), range(num_cats)))
 913      #   # compile counts of each category
 914      #   freqs = [0] * num_cats
 915      #   for v in list(d):
 916      #     freqs[val_order[v]] += 1
 917      #   # make list of the categories containing appropriate
 918      #   # counts of each in the order needed using a reverse lookup
 919      #   data_vals = list(val_order.keys())
 920      #   data_vals = [data_vals[i] for i, f in enumerate(freqs) if f > 0]
 921      # else: # any other data is easy!
 922      #   data_vals = sorted(d)
 923      #   freqs = [1] * len(data_vals)
 924      data_vals = sorted(d)
 925      key = self.tiling.tileable._get_legend_key_shapes(          # type: ignore
 926        geom, [1] * len(data_vals), rot, False)
 927      key_tiles.extend(key)
 928      vals.extend(data_vals)
 929      n = len(data_vals)
 930      ids.extend([id] * n)
 931      unique_ids.append(id)
 932      rots.extend([rot] * n)
 933    # finally make up a data table with all the data in all the columns. This
 934    # allows us to reuse the tiling_utils.plot_subsetted_gdf() function. To be
 935    # clear: all the data from all variables are added in all columns but when
 936    # sent for plotting only the subset associated with each tile_id will get
 937    # plotted. It's wasteful of space... but note that the same is true of the
 938    # original data - each tile_id has data for all the variables even if it's
 939    # not being used to plot them: tables gonna table!
 940    key_data = {}
 941    for ID in unique_ids:
 942      key_data[self.vars_to_map[self.ids_to_map.index(ID)]] = vals
 943    key_gdf = gpd.GeoDataFrame(
 944      data = key_data | {"tile_id": ids, "rotation": rots},
 945      crs = self.map.crs,
 946      geometry = gpd.GeoSeries(key_tiles))
 947    key_gdf.geometry = key_gdf.rotate(self.tiling.rotation, origin = (0, 0))
 948    return key_gdf
 949
 950
 951  def explore(self) -> None:
 952    """TODO: add wrapper to make tiled web map via geopandas.explore.
 953    """
 954    return None
 955
 956
 957  def _set_colourspecs(self) -> None:
 958    """Set _colourspecs dictionary based on instance member variables.
 959
 960    Set up is to the extent it is possible to follow user requested variables.
 961    Each requested `ids_to_map` item keys a dictionary in `_colourspecs` which
 962    contains the `column`, `cmap`, `scheme`, `categorical`, and `k` parameters
 963    to be passed on for use by the `geopandas.GeoDataFrame.plot()` calls in the
 964    _plot_subsetted_gdf()` method.
 965
 966    This is the place to make 'smart' adjustments to how user requests for map
 967    styling are handled.
 968    """
 969    numeric_columns = list(self.map.select_dtypes(
 970      include = ("float", "int")).columns)
 971    # note that some numeric columns can be considered categorical
 972    categorical_columns = list(self.map.select_dtypes(
 973      include = ("category", "int")).columns)
 974    try:
 975      if isinstance(self.ids_to_map, str):
 976        # wrap a single string in a list - this would be an unusual request...
 977        if self.ids_to_map in list(self.map.tile_id):
 978          print("""You have only requested a single attribute to map.
 979                That's fine, but perhaps not what you intended?""")
 980          self.ids_to_map = [self.ids_to_map]
 981        else:
 982          raise KeyError(
 983            """You have requested a single non-existent attribute to map!""")
 984      elif self.ids_to_map is None or not isinstance(self.ids_to_map, Iterable):
 985        # default to using all of them in order
 986        print("""No tile ids provided: setting all of them!""")
 987        self.ids_to_map = sorted(list(set(self.map.tile_id)))
 988
 989      if self.vars_to_map is None or not isinstance(self.vars_to_map, Iterable):
 990        self.vars_to_map = []
 991        if len(numeric_columns) == 0:
 992          # if there are none then we can't do it
 993          raise IndexError("""Attempting to set default variables, but
 994                          there are no numeric columns in the data!""")
 995        if len(numeric_columns) < len(self.ids_to_map):
 996          # if there are fewer available than we need then repeat some
 997          print("""Fewer numeric columns in the data than elements in the 
 998                tile unit. Reusing as many as needed to make up the numbers""")
 999          reps = len(self.ids_to_map) // len(numeric_columns) + 1
1000          self.vars_to_map = (numeric_columns * reps)[:len(self.ids_to_map)]
1001        elif len(numeric_columns) > len(self.ids_to_map):
1002          # if there are more than we need let the user know, but trim list
1003          print("""Note that you have supplied more variables to map than 
1004                there are distinct elements in the tile unit. Ignoring the
1005                extras.""")
1006          self.vars_to_map = numeric_columns[:len(self.ids_to_map)]
1007        else:
1008          self.vars_to_map = numeric_columns
1009      # print(f"{self.vars_to_map=}")
1010
1011      if self.categoricals is None or not isinstance(self.categoricals, Iterable):
1012        # provide a set of defaults
1013        self.categoricals = [col not in numeric_columns for col in self.vars_to_map]
1014      # print(f"{self.categoricals=}")
1015      
1016      if isinstance(self.schemes_to_use, str):
1017        self.schemes_to_use = [self.schemes_to_use] * len(self.ids_to_map)
1018      elif self.schemes_to_use is None or not isinstance(self.schemes_to_use, Iterable):
1019        # provide a set of defaults
1020        self.schemes_to_use = [None if cat else "EqualInterval"
1021                              for cat in self.categoricals]
1022      # print(f"{self.schemes_to_use=}")
1023      
1024      if isinstance(self.colors_to_use, str):
1025        self.colors_to_use = [self.colors_to_use] * len(self.ids_to_map)
1026      elif self.colors_to_use is None or not isinstance(self.colors_to_use, Iterable):
1027        # provide starter defaults
1028        print(f"""No colour maps provided! Setting some defaults.""")
1029        self.colors_to_use = ["Set1" if cat else "Reds" for cat in self.categoricals]
1030      for i, (col, cat) in enumerate(zip(self.colors_to_use, self.categoricals)):
1031        if cat and col not in CMAPS_CATEGORICAL:
1032          self.colors_to_use[i] = CMAPS_CATEGORICAL[i]
1033        # we'll allow diverging schemes for now...
1034        elif not cat and col not in CMAPS_SEQUENTIAL and col not in CMAPS_DIVERGING:
1035          self.colors_to_use[i] = CMAPS_SEQUENTIAL[i]
1036      # print(f"{self.colors_to_use=}")
1037
1038      if isinstance(self.n_classes, int):
1039        if self.n_classes == 0:
1040          self.n_classes = [255] * len(self.ids_to_map)
1041        else:
1042          self.n_classes = [self.n_classes] * len(self.ids_to_map)
1043      elif self.n_classes is None or not isinstance(self.n_classes, Iterable):
1044        # provide a set of defaults
1045        self.n_classes = [None if cat else 100 for cat in self.categoricals] 
1046      # print(f"{self.n_classes=}")
1047
1048    except IndexError as e:
1049      e.add_note("""One or more of the supplied lists of mapping settings is 
1050                 an inappropriate length""")
1051      raise
1052
1053    self._colourspecs = {
1054      ID: {"column": v,
1055           "cmap": c,
1056           "categorical": cat,
1057           "scheme": s,
1058           "k": k}
1059      for ID, v, c, cat, s, k 
1060      in zip(self.ids_to_map,
1061             self.vars_to_map,
1062             self.colors_to_use,
1063             self.categoricals,
1064             self.schemes_to_use,
1065             self.n_classes, strict = False)}