A set of functions for replacing a collection of polygons with hexagons ...

hex_my_map_up.py

import geopandas as gpd
from shapely import geometry
from shapely.affinity import translate, scale
import numpy as np
from matplotlib import pyplot as plt
from shapely.ops import cascaded_union
from scipy.optimize import linear_sum_assignment

# Define the hexagonal tile as a shapely polygon
def gen_hex(h=1, off=(0,0), orig=(0,0)):
    """
    off is the indes of this tile (col, row) (with the odd numbered cols shifted down by half a tile)
    orig is the starting location of the array
    """
    # The overall height of the polygon
    w = 2*h/3**0.5
    # the basic shape used this as a guide: https://www.redblobgames.com/grids/hexagons/
    pts = [
        [0, h/2],
        [w/4, h],
        [3*w/4, h],
        [w, h/2],
        [3*w/4, 0],
        [w/4, 0],
    ]
    # make the shapely pol
    pol = geometry.Polygon(pts)
    pol = translate(pol, 0, -h/2)
    # move the bottom left pol to the bottom left of the area we want to cover
    pol = translate(pol, orig[0], orig[1])
    # depending on what part of the array we are in we offset
    # NB the horozontal offset depends on the col we are in
    pol = translate(pol, 3/4*w*off[0], h*(off[1] - (off[0]%2)/2))
    return pol

# Generate an array of hexagonal tiles in a geopandas DF approsimating the shape of
# the geometry in a shapely polygon
def gen_grid(outline, divs):
    """
    outline is the polygon we wish to approximate
    divs is the number of polygons we want our array to have
    """
    # what area are we trying to cover?
    bounds = outline.bounds
    # let's start with an array 2 hexagons high
    h = (bounds[3] - bounds[1])/2
    # this is to initialize the binary search below
    d_h = h/2
    # binary search for the right hex size to get the desired number of divs
    for i in range(100):
        # width of a tile
        w = 2*h/3**0.5
        # number of tiles in the arrar (with extras due to the offset in the alternating
        # columns and because, better too many than too few
        hig = int((bounds[3] - bounds[1])/h//1+2)
        wid = int((bounds[2] - bounds[0])/w*4/3//1+2)
        # an index list for our hexs
        ind = [(i,j) for i in range(wid) for j in range(hig)]
        # generate the GDF
        grid = gpd.GeoDataFrame({
            'i':[p[0] for p in ind],
            'j':[p[1] for p in ind],
            'geometry':[gen_hex(h=h, off=p, orig=(bounds[0], bounds[1])) for p in ind],
        })
        # clip the grid to the tiles that overlap the map (by centroid)
        grid = grid.loc[grid.geometry.centroid.intersects(outline)]
        #print(f'Grid size: {grid.shape[0]}')
        # Now do a binary search for the hex height that produces the desired
        # number of divs, thank god for the intermediate value theorem
        if grid.shape[0] == divs:
            # when we have the right number of hex's stop looking
            break
        elif grid.shape[0] < divs:
            h = h - d_h
            d_h = d_h/2
            continue
        else:
            h = h + d_h
            d_h = d_h/2
            continue
    return grid.reset_index(drop=True)

def map_to_hex_grid(map_gdf):
    """
    take a GDF map_gdf and produce a hex tile for each row
    associate the tiles to the rows so that neighbours stay close
    NB the algorithm for associating the tiles is probably not 
    the right one ... seems to work mostly
    """
    # how many hexs do we need?
    n_areas = map_gdf.shape[0]
    # union of all the geom in the GDF
    outline = cascaded_union(map_gdf.geometry.buffer(10).simplify(10))
    # generate the required grid
    grid = gen_grid(outline, divs=n_areas)
    # calculate the distance^8 from all GDF row geoms to the hexs
    # NB distance^8 because it seems to work by trial and error ... this could do with improvement
    # NB2 vectorized AF FTW
    x1 = np.array(map_gdf.geometry.centroid.x).repeat(n_areas).reshape((n_areas,n_areas))
    y1 = np.array(map_gdf.geometry.centroid.y).repeat(n_areas).reshape((n_areas,n_areas))
    x2 = np.array(grid.geometry.centroid.x).repeat(n_areas).reshape((n_areas,n_areas)).T
    y2 = np.array(grid.geometry.centroid.y).repeat(n_areas).reshape((n_areas,n_areas)).T
    dists = np.power(np.square(x1-x2) + np.square(y1-y2), 4)
    # allocate the tiles to the GDF geoms to minimize the sum of the dists^8
    row_ind, col_ind = linear_sum_assignment(dists)
    # for the merge below
    grid.loc[:, 'key'] = row_ind
    map_gdf.loc[:, 'key'] = col_ind
    # attach the hexs to the GDF
    grid_assoc = grid.merge(
        map_gdf,
        on='key',
        suffixes=('', '_map'),
    )
    # time to clean 'ouse!
    grid_assoc.drop(columns=['i', 'j', 'key'], inplace=True)
    return gpd.GeoDataFrame(grid_assoc)

# Just some plotting to see what it did
def inspect(map_gdf, grid_assoc, label_col, arrows=True, size=12):
    fl = 'grid_assoc'
    plt.close(fl)
    fig, ax = plt.subplots(figsize=(size,size))
    fig.set_label(fl)

    ax_ = grid_assoc.plot(ax=ax, facecolor="none", edgecolor="black")
    grid_assoc.apply(
        lambda x: ax_.annotate(
            text=x[label_col],
            xy=x.geometry.centroid.coords[0],
            ha='center'
        ),axis=1);

    if arrows:
        for r in grid_assoc.iterrows():
            r = r[1]
            ax.plot(
                (r.geometry.centroid.x, r.geometry_map.centroid.x),
                (r.geometry.centroid.y, r.geometry_map.centroid.y), 'r'
            )
    for r in grid_assoc.iterrows():
        r = r[1]
        ax.plot(
            (r.geometry_map.centroid.x),
            (r.geometry_map.centroid.y), 'k.'
        )
    map_gdf.plot(alpha=0.5, ax=ax, column=label_col)