Working With Geometries

In the course of working with legislative redistricting data, it is inevitable that we will have to work with files that contain geometries. Most often, these geometries come in the form of shapefiles, which, while nice in theory, can be a bit of a pain to work with in practice. For now, we will focus on the basics of working with geometries, but the interested reader is encouraged to explore our partnered library, maup, which is specifically designed to help fix tricky geometry problems. Specifically, in the event that you are working with a shapefile and run in to an error of the flavour:

UserWarning: Found overlaps among the given polygons.
Indices of overlaps: {(887, 892), (893, 915), (892, 914), (887, 893)}

or

UserWarning: Found islands (degree-0 nodes). Indices of islands: {2552, 3107}
"Found islands (degree-0 nodes). Indices of islands: {}".format(islands)

then you should consider consulting the maup documentation to see if it can help you out.

Loading and Running a Plan


For this example, we will make use of a Minnesota GeoJSON file that contains the geometries of the state’s precincts (you will need to unzip the above folder to get to the file – it’s a bit large). We will follow a similar workflow to what we already covered in the ReCom section, but with an eye towards some of the conveniences afforded by GeographicPartition objects. As always, we’ll start with the imports:

import matplotlib.pyplot as plt
from gerrychain import (
    Partition,
    Graph,
    MarkovChain,
    updaters,
    constraints,
    accept,
    GeographicPartition,
)
from gerrychain.proposals import build_recom_proposal_fn
from gerrychain.tree import bipartition_tree
from gerrychain.constraints import contiguous
import pandas

Matplotlib is building the font cache; this may take a moment.

And now we load the graph from the GeoJSON file

import zipfile

with zipfile.ZipFile("MN.zip") as z:
    z.extractall()

graph = Graph.from_file("MN_precincts.geojson")

as well as create our chain, with its initial partition and updaters

recom_chain = MarkovChain(
    total_steps=20,
    rng=42,
)

recom_chain.initial_partition = GeographicPartition(
    graph,
    assignment="CONGDIST",
)

recom_chain.add_updaters(
    {
        "population": updaters.Tally("TOTPOP", alias="population"),
        "cut_edges": updaters.cut_edges,
        "perimeter": updaters.perimeter,
        "area": updaters.Tally("area", alias="area"),
    }
)

The observant reader will notice that we have added two new updaters, perimeter, and area, [1] and we are now using the GeographicPartition class instead of the Partition class. The GeographicPartition class is a subclass of the Partition class that allows us the capability of working with geometries throughout our Markov chain, and the perimeter and area updaters are examples of such a geometric updater that was previously unavailable to us. These updaters are necessary for monitoring things like geometric compactness and area via metrics such as the Polsby-Popper test. [2]

And now it is time for one of the first conveniences of the GeographicPartition class: we can plot our map and see the initial partition!

recom_chain.initial_partition.plot()

<Axes: >
../../_images/a120e8b2eec8a975596c8d14fea90ff136764fd9a1ff4b7bdc1fc288e71645e3.png

Of course, this isn’t very pretty, so let’s pass it some additional arguments to things a bit nicer:

fig, ax = plt.subplots(figsize=(8, 8))
ax.set_yticks([])
ax.set_xticks([])
ax.set_title("Initial Partition in MN")
recom_chain.initial_partition.plot(ax=ax, cmap="tab20c")

<Axes: title={'center': 'Initial Partition in MN'}>
../../_images/42297f3b17e39758c328e04de8dc2abeeba7c1cf4d6b529f6da87090f623853b.png

Under the hood, the plot method is using the geodataframe.plot method from geopandas to plot the geometries, and all of this is built on top of matplotlib, so most of the standard methods for modifying a matplotlib plot will work here as well.

The chain already knows where it starts and what to track, so all that is left is to tell it how to propose a plan, which plans are valid, and when to accept one:

ideal_population = sum(recom_chain.initial_partition["population"].values()) / len(
    recom_chain.initial_partition
)

recom_chain.proposal_fn = build_recom_proposal_fn(
    pop_col="TOTPOP",
    pop_target=ideal_population,
    epsilon=0.01,
)
recom_chain.add_constraint(contiguous)
recom_chain.acceptance_fn = accept.always_accept

The next cell builds an interactive viewer for watching the chain. Its Back and Forward buttons run entirely in the browser, so they also work in the rendered documentation.

from io import BytesIO

from matplotlib.animation import ArtistAnimation
from PIL import Image
from IPython.display import HTML

frames = []
district_data = []

for i, partition in enumerate(recom_chain):
    for district_name in partition["perimeter"]:
        district_data.append(
            (
                i,
                district_name,
                partition["population"][district_name],
                partition["perimeter"][district_name],
                partition["area"][district_name],
            )
        )

    with BytesIO() as buffer:
        fig, ax = plt.subplots(figsize=(10, 10))
        partition.plot(ax=ax, cmap="tab20")
        ax.set_xticks([])
        ax.set_yticks([])
        fig.savefig(buffer, format="png", bbox_inches="tight", pad_inches=0)
        frames.append(Image.open(buffer).copy())
        plt.close(fig)

df = pandas.DataFrame(
    district_data,
    columns=["step", "district_name", "population", "perimeter", "area"],
)

fig, ax = plt.subplots(figsize=(8, 8))
ax.axis("off")
fig.subplots_adjust(left=0, right=1, bottom=0, top=1)
artists = [[ax.imshow(frame, animated=True)] for frame in frames]
animation = ArtistAnimation(fig, artists, interval=500)
plt.close(fig)
HTML(animation.to_jshtml())