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.drop(columns = region_vars, inplace = True)
 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.drop(columns = ["joinUID"], inplace = True)
 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    tiled_map.drop(columns = [id_var], inplace = True)
 398    self.region.drop(columns = [id_var], inplace = True)
 399
 400    # if we've retained tiles and want 'clean' edges, then clip
 401    # note that this step is slow: geopandas unary_unions the clip layer
 402    if prioritise_tiles and not ragged_edges:
 403      tiled_map.sindex
 404      tiled_map = tiled_map.clip(self.region_union)
 405      if debug:
 406        print(f"""STEP A7/B3: clip map to region: {perf_counter() - t7:.3f}""")
 407
 408    tm = TiledMap()
 409    tm.tiling = self
 410    tm.map = tiled_map
 411    return tm
 412
 413
 414  def _setup_region_DZID(self) -> str:
 415    """Return guaranteed-unique new attribute name for self.region dataframe.
 416
 417    Avoids a name clash with any existing attribute in the dataframe.
 418
 419    Returns:
 420      str: name of the added attribute.
 421
 422    """
 423    dzid = "DZID"
 424    i = 0
 425    while dzid in self.region.columns:
 426      dzid = "DZID" + str(i)
 427      i = i + 1
 428    self.region[dzid] = list(range(self.region.shape[0]))
 429    return dzid
 430
 431
 432  def make_tiling(self) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
 433    """Tile the region with a tile unit, returning a GeoDataFrame.
 434
 435    Returns:
 436      geopandas.GeoDataFrame: a GeoDataFrame of the region tiled with the
 437        tile unit.
 438
 439    """
 440    # we assume the geometry column is called geometry so make it so...
 441    if self.region.geometry.name != "geometry":
 442      self.region.rename_geometry("geometry", inplace = True)
 443
 444    # chain list of lists of GeoSeries geometries to list of geometries
 445    tiles = itertools.chain(*[
 446      self.tileable.tiles.geometry.translate(p.x, p.y)
 447      for p in self.grid.points])
 448    prototiles = itertools.chain(*[
 449      self.tileable.prototile.geometry.translate(p.x, p.y)
 450      for p in self.grid.points])
 451    # replicate the tile ids
 452    tile_ids = list(self.tileable.tiles.tile_id) * len(self.grid.points)
 453    prototile_ids = list(range(len(self.grid.points)))
 454    tile_prototile_ids = sorted(prototile_ids * self.tileable.tiles.shape[0])
 455    tiles_gs = gpd.GeoSeries(list(tiles))
 456    prototiles_gs = gpd.GeoSeries(list(prototiles))
 457    # assemble and return as GeoDataFrames
 458    tiles_gdf = gpd.GeoDataFrame(
 459      data = {"tile_id": tile_ids, "prototile_id": tile_prototile_ids},
 460      geometry = tiles_gs, crs = self.tileable.crs)
 461    prototiles_gdf = gpd.GeoDataFrame(
 462      data = {"prototile_id": prototile_ids},
 463      geometry = prototiles_gs, crs = self.tileable.crs)
 464    return tiles_gdf, prototiles_gdf
 465
 466
 467  def rotated(self,
 468              rotation:float = 0.0,
 469              ) -> tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]:
 470    """Return the stored tiling rotated.
 471
 472    The stored tiling never changes and
 473    if it was originally made with a Tileable that was rotated it will retain
 474    that rotation. The requested rotation is _additional_ to that baseline
 475    rotation.
 476
 477    Args:
 478      rotation (float, optional): Rotation angle in degrees.
 479        Defaults to None.
 480
 481    Returns:
 482      gpd.GeoDataFrame: Rotated tiling.
 483
 484    """
 485    if self.tiles is None:
 486      self.tiles = self.make_tiling()[0]
 487    if rotation == 0:
 488      return self.tiles, self.prototiles
 489    tiles = gpd.GeoDataFrame(
 490      data = {"tile_id": self.tiles.tile_id,
 491              "prototile_id": self.tiles.tile_id},
 492      crs = self.tiles.crs,
 493      geometry = self.tiles.geometry.rotate(
 494        rotation, origin = self.grid.centre))
 495    prototiles = gpd.GeoDataFrame(
 496      data = {"prototile_id": self.prototiles.prototile_id},
 497      crs = self.prototiles.crs,
 498      geometry = self.prototiles.geometry.rotate(
 499        rotation, origin = self.grid.centre))
 500    self.rotation = rotation
 501    return tiles, prototiles
 502
 503
 504@dataclass(slots=True)
 505class TiledMap:
 506  """Class representing a tiled map.
 507
 508  Should not be accessed directly, but will
 509  be created by calling `Tiling.get_tiled_map()`. After creation the variables
 510  and colourmaps attributes can be set, and then `TiledMap.render()` called to
 511  make a map. Settable attributes are explained in documentation of the
 512  `TiledMap.render()` method.
 513
 514  Examples:
 515    Recommended usage is as follows. First, make a `TiledMap` from a `Tiling`
 516    object:
 517
 518      `tm = tiling.get_tiled_map(...)`
 519
 520    Some options in the `Tiling` constructor affect the map appearance. See
 521    `Tiling` for details.
 522
 523    Once a `TiledMap` object exists, set options on it, either when calling
 524    `TiledMap.render()` or explicitly, i.e.
 525
 526      tm.render(opt1 = val1, opt2 = val2, ...)
 527
 528    or
 529
 530      tm.opt1 = val1
 531      tm.opt2 = val2
 532      tm.render()
 533
 534    Option settings are persistent, i.e. unless a new `TiledMap` object is
 535    created the option settings have to be explicitly reset to new values on
 536    subsequent calls to `TiledMap.render()`.
 537
 538    The most important options are the `vars_map` and `colors_to_use` settings.
 539
 540    `vars_to_map` is a lost of the dataset variable names to match with
 541    `weavingspace.tileable.Tileable` elements with corresponding (ordered)
 542    tile_ids (usually "a", "b", etc.). If you need to control the match, then
 543    also supply `ids_to_map` in matching order. E.g.
 544
 545      tm.ids_to_map = ['d', 'c', 'b', 'a']
 546      tm.vars_to_map = ['x1', 'x2', 'x3', 'x4']
 547
 548    Note that this means that if you really want more than one element in the
 549    tiling to represent the same variable more than once, you can do that.
 550
 551    `colors_to_use` is a parallel list of named matplotlib colormaps,
 552
 553      tm.colors_to_use = ["Reds", "Blues", "Greys", "Purples"]
 554
 555    Similarly, you can specify the classification `schemes_to_use` (such as
 556    'quantiles') and the number of classes `n_classes` in each.
 557
 558    If data are categorical, this is flagged in the `categoricals` list of
 559    booleans, in which case an appropriate colour map should be used. There is
 560    currently no provision for control of which colour in a categorical
 561    colour map is applied to which variable level.
 562
 563    TODO: better control of categorical mapping schemes.
 564
 565  """
 566
 567  # these will be set at instantion by Tiling.get_tiled_map()
 568  tiling:Tiling = None
 569  """the Tiling with the required tiles"""
 570  map:gpd.GeoDataFrame = None
 571  """the GeoDataFrame on which this map is based"""
 572  ids_to_map:list[str] = None
 573  """tile_ids that are to be used to represent data"""
 574  vars_to_map:list[str] = None
 575  """dataset variables that are to be symbolised"""
 576  colors_to_use:list[str|list[str]] = None
 577  """list of matplotlib colormap names."""
 578  categoricals:list[bool] = None
 579  """list specifying if each variable is -- or is to be treated as --
 580  categorical"""
 581  schemes_to_use:list[str|None] = None
 582  """mapclassify schemes to use for each variable."""
 583  n_classes:list[int|None] = None
 584  """number of classes to apply; if set to 0 will be unclassed."""
 585  _colourspecs:dict[str,dict] = None
 586  """dictionary of dictionaries keyed by the items in `ids_to_use` with each
 587  dictionary forming additional kwargs to be supplied to geopandas.plot()."""
 588
 589  # the below parameters can be set either before calling self.render() or
 590  # passed in as parameters to self.render(). These are solely
 591  # `TiledMap.render()` options not geopandas plot options.
 592  legend:bool = True
 593  """whether or not to show a legend"""
 594  legend_zoom:float = 1.0
 595  """<1 zooms out from legend to show more context"""
 596  legend_dx:float = 0.
 597  """x shift of legend relative to the map"""
 598  legend_dy:float = 0.
 599  """y shift of legend relative to the map"""
 600  use_ellipse:bool = False
 601  """if True clips legend with an ellipse"""
 602  ellipse_magnification:float = 1.0
 603  """magnification to apply to clip ellipse"""
 604  radial_key:bool = False
 605  """if True use radial key even for ordinal/ratio data (normally these will be
 606  shown by concentric tile geometries)"""
 607  draft_mode:bool = False
 608  """if True plot the map coloured by tile_id"""
 609
 610  # the parameters below are geopandas.plot options which we intercept to
 611  # ensure they are applied appropriately when we plot a GDF
 612  figsize:tuple[float,float] = (20, 15)
 613  """maptlotlib figsize"""
 614  dpi:float = 72
 615  """dpi for bitmap formats"""
 616
 617  def render(
 618      self,
 619      **kwargs,
 620    ) -> Figure:
 621    """Render the current state to a map.
 622
 623    Note that TiledMap objects will usually be created by calling
 624    `Tiling.get_tiled_map()`.
 625
 626    Args:
 627      ids_to_map (list[str]): tile_ids to be used in the map. Defaults to None.
 628      vars_to_map (list[str]): dataset columns to be mapped. Defaults to None.
 629      colors_to_use (list[str]): list of matplotlib colormaps to be used.
 630        Defaults to None.
 631      categoricals (list[bool]): list of flags indicating if associated variable
 632        should be treated as categorical. Defaults to None.
 633      schemes_to_use (list[str]): list of strings indicating the mapclassify
 634        scheme to use e.g. 'equalinterval' or 'quantiles'. Defaults to None.
 635      n_classes (list[int]): list of ints indicating number of classes to use in
 636        map classification. Defaults to None.
 637      legend (bool): If True a legend will be drawn. Defaults to True.
 638      legend_zoom (float): Zoom factor to apply to the legend. Values <1
 639        will show more of the tile context. Defaults to 1.0.
 640      legend_dx (float): x shift to apply to the legend position in plot area
 641        relative units, i.e. 1.0 is full width of plot. Defaults to 0.0.
 642      legend_dy (float): x and y shift to apply to the legend position in plot
 643        area relative units, i.e. 1.0 is full height of plot. Defaults to 0.0.
 644      use_ellipse (bool): If True applies an elliptical clip to the legend.
 645        Defaults to False.
 646      ellipse_magnification (float): Magnification to apply to ellipse clipped
 647        legend. Defaults to 1.0.
 648      radial_key (bool): If True legend key for TileUnit maps will be based on
 649        radially dissecting the tiles, i.e. pie slices. Defaults to False.
 650      draft_mode (bool): If True a map of the tiled map coloured by tile_ids
 651        (and with no legend) is returned. Defaults to False.
 652      figsize (tuple[float,floar]): plot dimensions passed to geopandas.
 653        plot. Defaults to (20, 15).
 654      dpi (float): passed to pyplot.plot. Defaults to 72.
 655      **kwargs: other settings to pass to pyplot/geopandas.plot.
 656
 657    Returns:
 658      matplotlib.figure.Figure: figure on which map is plotted.
 659
 660    """
 661    plt.rcParams["pdf.fonttype"] = 42
 662    plt.rcParams["pdf.use14corefonts"] = True
 663    matplotlib.rcParams["pdf.fonttype"] = 42
 664
 665    to_remove = set()  # keep track of kwargs we use to setup TiledMap
 666    # kwargs with no corresponding class attribute will be discarded
 667    # because we are using slots, we have to use setattr() here
 668    for k, v in kwargs.items():
 669      if k in self.__slots__:
 670        setattr(self, k, v)
 671        to_remove.add(k)
 672    # remove any them so we don't pass them on to pyplot and get errors
 673    for k in to_remove:
 674      del kwargs[k]
 675
 676    if self.draft_mode:
 677      fig = plt.figure(figsize = self.figsize)
 678      ax = fig.add_subplot(111)
 679      self.map.plot(ax = ax, column = "tile_id", cmap = "tab20", **kwargs)
 680      ax.set_axis_off()
 681      return fig
 682
 683    if self.legend:
 684      # this sizing stuff is rough and ready for now, possibly forever...
 685      reg_w, reg_h, *_ = \
 686        tiling_utils.get_width_height_left_bottom(self.map.geometry)
 687      tile_w, tile_h, *_ = \
 688        tiling_utils.get_width_height_left_bottom(
 689          self.tiling.tileable._get_legend_tiles().rotate(
 690            self.tiling.rotation, origin = (0, 0)))
 691      sf_w, sf_h = reg_w / tile_w / 3, reg_h / tile_h / 3
 692      gskw = {"height_ratios": [sf_h * tile_h, reg_h - sf_h * tile_h],
 693              "width_ratios":  [reg_w, sf_w * tile_w]}
 694
 695      fig, axes = plt.subplot_mosaic([["map", "legend"], ["map", "."]],
 696                                        gridspec_kw = gskw,
 697                                        figsize = self.figsize,
 698                                        layout = "constrained", **kwargs)
 699    else:
 700      fig, axes = plt.subplots(1, 1, figsize = self.figsize,
 701                                  layout = "constrained", **kwargs)
 702
 703    self._plot_map(axes, **kwargs)
 704    return fig
 705
 706
 707  def _plot_map(
 708      self,
 709      ax:plt.Axes,
 710      **kwargs,
 711    ) -> None:
 712    """Plot map to the supplied axes.
 713
 714    Args:
 715      axes (plt.Axes): axes on which maps will be drawn.
 716      kwargs (dict): additional parameters to be passed to plot.
 717
 718    """
 719    self._set_colourspecs()
 720    bb = self.map.geometry.total_bounds
 721    if self.legend:
 722      if (self.legend_dx != 0 or self.legend_dx != 0):
 723        box = ax["legend"].get_position()
 724        box.x0 += self.legend_dx
 725        box.x1 += self.legend_dx
 726        box.y0 += self.legend_dy
 727        box.y1 += self.legend_dy
 728        ax["legend"].set_position(box)
 729      ax["map"].set_axis_off()
 730      ax["map"].set_xlim(bb[0], bb[2])
 731      ax["map"].set_ylim(bb[1], bb[3])
 732      self._plot_subsetted_gdf(ax["map"], self.map, **kwargs)
 733      self.plot_legend(ax = ax["legend"], **kwargs)
 734    else:
 735      ax.set_axis_off()
 736      ax.set_xlim(bb[0], bb[2])
 737      ax.set_ylim(bb[1], bb[3])
 738      self._plot_subsetted_gdf(ax, self.map, **kwargs)
 739
 740
 741  def _plot_subsetted_gdf(
 742      self,
 743      ax:plt.Axes,
 744      gdf:gpd.GeoDataFrame,
 745      grouping_var:str = "tile_id",
 746      **kwargs,
 747    ) -> None:
 748    """Repeatedly plot a gpd.GeoDataFrame based on a subsetting attribute.
 749
 750    NOTE: used to plot both the main map _and_ the legend, which is why
 751    a separate GeoDataframe is supplied and we don't just use self.map.
 752
 753    Args:
 754      ax (pyplot.Axes): axes to plot to.
 755      gdf (gpd.GeoDataFrame): the GeoDataFrame to plot.
 756
 757    Raises:
 758      Exception: if self.colourmaps cannot be parsed exception is raised.
 759
 760    """
 761    groups = gdf.groupby(grouping_var)
 762    for ID, cspec in self._colourspecs.items():
 763      subset = groups.get_group(ID)
 764      n_values = len(subset[cspec["column"]].unique())
 765      if not cspec["categorical"] and n_values == 1:
 766        print(f"""
 767              Only one level in variable {cspec['column']}, replacing requested
 768              colour map with single colour fill.
 769              """)
 770        cspec["color"] = matplotlib.colormaps.get(cspec["cmap"])(0.5)
 771        del cspec["column"]
 772        del cspec["cmap"]
 773        del cspec["scheme"]
 774      elif not cspec["categorical"] and n_values < cspec["k"]:
 775        cspec["k"] = n_values
 776      elif cspec["categorical"]:
 777        del cspec["scheme"]
 778      subset.plot(ax = ax, **cspec, **kwargs)
 779
 780
 781  def to_file(self, fname:str) -> None:
 782    """Output the tiled map to a layered GPKG file.
 783
 784    Currently delegates to `weavingspace.tiling_utils.write_map_to_layers()`.
 785
 786    Args:
 787      fname (str): Filename to write. Defaults to None.
 788
 789    """
 790    tiling_utils.write_map_to_layers(self.map, fname)
 791
 792
 793  def plot_legend(self, ax, **kwargs) -> None:
 794    """Plot a legend for this tiled map.
 795
 796    Args:
 797      ax (pyplot.Axes, optional): axes to draw legend. Defaults to None.
 798
 799    """
 800    ax.set_axis_off()
 801    legend_tiles = self.tiling.tileable._get_legend_tiles()
 802    # this is a bit hacky, but we will apply the rotation to text
 803    # annotation so for TileUnits which don't need it, reverse that now
 804    if isinstance(self.tiling.tileable, TileUnit):
 805      # note that this confuses type hinting because pandas silently assigns 
 806      # a scalar value to a Series
 807      legend_tiles.rotation = -self.tiling.rotation
 808
 809    legend_key = self._get_legend_key_gdf(legend_tiles)
 810    legend_tiles.geometry = legend_tiles.geometry.rotate(
 811      self.tiling.rotation, origin = (0, 0))
 812
 813    if self.use_ellipse:
 814      ellipse = tiling_utils.get_bounding_ellipse(
 815        legend_tiles.geometry, mag = self.ellipse_magnification)
 816      bb = ellipse.total_bounds
 817      c = ellipse.union_all().centroid
 818    else:
 819      bb = legend_tiles.geometry.total_bounds
 820      c = legend_tiles.geometry.union_all().centroid
 821
 822    # apply legend zoom - NOTE that this must be applied even
 823    # if self.legend_zoom is == 1...
 824    ax.set_xlim(c.x + (bb[0] - c.x) / self.legend_zoom,
 825                c.x + (bb[2] - c.x) / self.legend_zoom)
 826    ax.set_ylim(c.y + (bb[1] - c.y) / self.legend_zoom,
 827                c.y + (bb[3] - c.y) / self.legend_zoom)
 828
 829    for cs, tile, rotn in zip(self._colourspecs.values(),
 830                              legend_tiles.geometry,
 831                              legend_tiles.rotation):
 832      c = tile.centroid
 833      ax.annotate(cs["column"], xy = (c.x, c.y),
 834                  ha = "center", va = "center",
 835                  rotation_mode = "anchor",
 836                  # adjust rotation to favour text reading left to right
 837                  rotation = (rotn + self.tiling.rotation + 90) % 180 - 90,
 838                  bbox = {"lw": 0, "fc": "#ffffff60"})
 839
 840    # now plot background; we include the central tiles, since in
 841    # the weave case these may not match the legend tiles
 842    context_tiles = self.tiling.tileable \
 843      .get_local_patch(r = 2, include_0 = True) \
 844      .geometry.rotate(self.tiling.rotation, origin = (0, 0))
 845    if self.use_ellipse:
 846      context_tiles.clip(ellipse, keep_geom_type = False).plot(
 847        ax = ax, fc = "#9F9F9F3F", lw = .35)
 848      tiling_utils.get_tiling_edges(context_tiles.geometry).clip(
 849        ellipse, keep_geom_type = True).plot(
 850          ax = ax, ec = "#5F5F5F", lw = .35)
 851    else:
 852      context_tiles.plot(ax = ax, fc = "#9F9F9F3F",
 853                         ec = "#5F5F5F", lw = .35)
 854      tiling_utils.get_tiling_edges(context_tiles.geometry).plot(
 855        ax = ax, ec = "#5F5F5F", lw = .35)
 856
 857    # plot the legend key tiles (which include the data)
 858    self._plot_subsetted_gdf(ax, legend_key, **kwargs)
 859
 860
 861  def _get_legend_key_gdf(self, tiles:gpd.GeoDataFrame) -> gpd.GeoDataFrame:
 862    """Return tiles dissected and with data assigned for use as a legend.
 863
 864    'Dissection' is handled differently by `WeaveUnit` and `TileUnit`
 865    objects and delegated to either `WeaveUnit._get_legend_key_shapes()`
 866    or `TileUnit._get_legend_key_shapes()`.
 867
 868    Args:
 869      tiles (gpd.GeoDataFrame): the legend tiles.
 870
 871    Returns:
 872      gpd.GeoDataFrame:  with tile_id, variables and rotation
 873        attributes, and geometries of Tileable tiles sliced into a
 874        colour ramp or set of nested tiles.
 875
 876    """
 877    key_tiles = []   # set of tiles to form a colour key (e.g. a ramp)
 878    ids = []         # tile_ids applied to the keys
 879    unique_ids = []  # list of each tile_id used in order
 880    vals = []        # the data assigned to the key tiles
 881    rots = []        # rotation of each key tile
 882    # subsets = self.map.groupby("tile_id")
 883    for (id, cspec), geom, rot in zip(self._colourspecs.items(),
 884                                      tiles.geometry, 
 885                                      tiles.rotation):
 886      d = list(self.map.loc[self.map.tile_id == id][cspec["column"]])
 887      # if the data are categorical then it's complicated...
 888      # if cs["categorical"]:
 889      #   radial = True and self.radial_key
 890      #   # desired order of categorical variable is the
 891      #   # color maps dictionary keys
 892      #   num_cats = len(cmap)
 893      #   val_order = dict(zip(cmap.keys(), range(num_cats)))
 894      #   # compile counts of each category
 895      #   freqs = [0] * num_cats
 896      #   for v in list(d):
 897      #     freqs[val_order[v]] += 1
 898      #   # make list of the categories containing appropriate
 899      #   # counts of each in the order needed using a reverse lookup
 900      #   data_vals = list(val_order.keys())
 901      #   data_vals = [data_vals[i] for i, f in enumerate(freqs) if f > 0]
 902      # else: # any other data is easy!
 903      #   data_vals = sorted(d)
 904      #   freqs = [1] * len(data_vals)
 905      data_vals = sorted(d)
 906      key = self.tiling.tileable._get_legend_key_shapes(          # type: ignore
 907        geom, [1] * len(data_vals), rot, False)
 908      key_tiles.extend(key)
 909      vals.extend(data_vals)
 910      n = len(data_vals)
 911      ids.extend([id] * n)
 912      unique_ids.append(id)
 913      rots.extend([rot] * n)
 914    # finally make up a data table with all the data in all the columns. This
 915    # allows us to reuse the tiling_utils.plot_subsetted_gdf() function. To be
 916    # clear: all the data from all variables are added in all columns but when
 917    # sent for plotting only the subset associated with each tile_id will get
 918    # plotted. It's wasteful of space... but note that the same is true of the
 919    # original data - each tile_id has data for all the variables even if it's
 920    # not being used to plot them: tables gonna table!
 921    key_data = {}
 922    for ID in unique_ids:
 923      key_data[self.vars_to_map[self.ids_to_map.index(ID)]] = vals
 924    key_gdf = gpd.GeoDataFrame(
 925      data = key_data | {"tile_id": ids, "rotation": rots},
 926      crs = self.map.crs,
 927      geometry = gpd.GeoSeries(key_tiles))
 928    key_gdf.geometry = key_gdf.rotate(self.tiling.rotation, origin = (0, 0))
 929    return key_gdf
 930
 931
 932  def explore(self) -> None:
 933    """TODO: add wrapper to make tiled web map via geopandas.explore.
 934    """
 935    return None
 936
 937
 938  def _set_colourspecs(self) -> None:
 939    """Set _colourspecs dictionary based on instance member variables.
 940
 941    Set up is to the extent it is possible to follow user requested variables.
 942    Each requested `ids_to_map` item keys a dictionary in `_colourspecs` which
 943    contains the `column`, `cmap`, `scheme`, `categorical`, and `k` parameters
 944    to be passed on for use by the `geopandas.GeoDataFrame.plot()` calls in the
 945    _plot_subsetted_gdf()` method.
 946
 947    This is the place to make 'smart' adjustments to how user requests for map
 948    styling are handled.
 949    """
 950    numeric_columns = list(self.map.select_dtypes(
 951      include = ("float", "int")).columns)
 952    # note that some numeric columns can be considered categorical
 953    categorical_columns = list(self.map.select_dtypes(
 954      include = ("category", "int")).columns)
 955    try:
 956      if isinstance(self.ids_to_map, str):
 957        # wrap a single string in a list - this would be an unusual request...
 958        if self.ids_to_map in list(self.map.tile_id):
 959          print("""You have only requested a single attribute to map.
 960                That's fine, but perhaps not what you intended?""")
 961          self.ids_to_map = [self.ids_to_map]
 962        else:
 963          raise KeyError(
 964            """You have requested a single non-existent attribute to map!""")
 965      elif self.ids_to_map is None or not isinstance(self.ids_to_map, Iterable):
 966        # default to using all of them in order
 967        print("""No tile ids provided: setting all of them!""")
 968        self.ids_to_map = sorted(list(set(self.map.tile_id)))
 969
 970      if self.vars_to_map is None or not isinstance(self.vars_to_map, Iterable):
 971        self.vars_to_map = []
 972        if len(numeric_columns) == 0:
 973          # if there are none then we can't do it
 974          raise IndexError("""Attempting to set default variables, but
 975                          there are no numeric columns in the data!""")
 976        if len(numeric_columns) < len(self.ids_to_map):
 977          # if there are fewer available than we need then repeat some
 978          print("""Fewer numeric columns in the data than elements in the 
 979                tile unit. Reusing as many as needed to make up the numbers""")
 980          reps = len(self.ids_to_map) // len(numeric_columns) + 1
 981          self.vars_to_map = (numeric_columns * reps)[:len(self.ids_to_map)]
 982        elif len(numeric_columns) > len(self.ids_to_map):
 983          # if there are more than we need let the user know, but trim list
 984          print("""Note that you have supplied more variables to map than 
 985                there are distinct elements in the tile unit. Ignoring the
 986                extras.""")
 987          self.vars_to_map = numeric_columns[:len(self.ids_to_map)]
 988        else:
 989          self.vars_to_map = numeric_columns
 990      # print(f"{self.vars_to_map=}")
 991
 992      if self.categoricals is None or not isinstance(self.categoricals, Iterable):
 993        # provide a set of defaults
 994        self.categoricals = [col not in numeric_columns for col in self.vars_to_map]
 995      # print(f"{self.categoricals=}")
 996      
 997      if isinstance(self.schemes_to_use, str):
 998        self.schemes_to_use = [self.schemes_to_use] * len(self.ids_to_map)
 999      elif self.schemes_to_use is None or not isinstance(self.schemes_to_use, Iterable):
1000        # provide a set of defaults
1001        self.schemes_to_use = [None if cat else "EqualInterval"
1002                              for cat in self.categoricals]
1003      # print(f"{self.schemes_to_use=}")
1004      
1005      if isinstance(self.colors_to_use, str):
1006        self.colors_to_use = [self.colors_to_use] * len(self.ids_to_map)
1007      elif self.colors_to_use is None or not isinstance(self.colors_to_use, Iterable):
1008        # provide starter defaults
1009        print(f"""No colour maps provided! Setting some defaults.""")
1010        self.colors_to_use = ["Set1" if cat else "Reds" for cat in self.categoricals]
1011      for i, (col, cat) in enumerate(zip(self.colors_to_use, self.categoricals)):
1012        if cat and col not in CMAPS_CATEGORICAL:
1013          self.colors_to_use[i] = CMAPS_CATEGORICAL[i]
1014        # we'll allow diverging schemes for now...
1015        elif not cat and col not in CMAPS_SEQUENTIAL and col not in CMAPS_DIVERGING:
1016          self.colors_to_use[i] = CMAPS_SEQUENTIAL[i]
1017      # print(f"{self.colors_to_use=}")
1018
1019      if isinstance(self.n_classes, int):
1020        if self.n_classes == 0:
1021          self.n_classes = [255] * len(self.ids_to_map)
1022        else:
1023          self.n_classes = [self.n_classes] * len(self.ids_to_map)
1024      elif self.n_classes is None or not isinstance(self.n_classes, Iterable):
1025        # provide a set of defaults
1026        self.n_classes = [None if cat else 100 for cat in self.categoricals] 
1027      # print(f"{self.n_classes=}")
1028
1029    except IndexError as e:
1030      e.add_note("""One or more of the supplied lists of mapping settings is 
1031                 an inappropriate length""")
1032      raise
1033
1034    self._colourspecs = {
1035      ID: {"column": v,
1036           "cmap": c,
1037           "categorical": cat,
1038           "scheme": s,
1039           "k": k}
1040      for ID, v, c, cat, s, k 
1041      in zip(self.ids_to_map,
1042             self.vars_to_map,
1043             self.colors_to_use,
1044             self.categoricals,
1045             self.schemes_to_use,
1046             self.n_classes, strict = False)}