Fixing topological issues, overlaps, and gaps

Precinct shapefiles are often created by stitching together collections of precinct geometries sourced from different counties or different years. As a result, the shapefile often has gaps or overlaps between precincts where the different sources disagree about the boundaries. (And by “often,” we mean “for almost every shapefile that isn’t produced by the U.S. Census Burueau.”) As we saw in the examples above, these issues can pose problems when moving data between shapefiles.

Even when working with a single shapefile, gaps and overlaps may cause problems if you are interested in working with the adjacency graph of the precincts. This adjacency information is especially important when studying redistricting, because districts are almost always expected to be contiguous.

Before doing anything else, it is wise to understand the current status of a shapefile with regard to topological issues. maup provides a doctor function to diagnose gaps, overlaps, and invalid geometries. If a shapefile has none of these issues, maup.doctor returns a value of True; otherwise it returns False along with a brief summary of the problems that it found.

The blocks shapefile, like most shapefiles from the Census Bureau, is clean:

>>> maup.doctor(blocks)
True

The old precincts shapefile, however, has some minor issues:

>>> maup.doctor(old_precincts)
There are 2 overlaps.
There are 3 holes.
False

As of version 2.0.0, maup provides two repair functions with a variety of options for fixing these issues:

  1. quick_repair is the new name for the autorepair function from version 1.x (and autorepair still works as a synonym). This function makes fairly simplistic repairs to gaps and overlaps:

    • Any polygon \(Q\) created by the overlapping intersection of two geometries \(P_1\) and \(P_2\) is removed from both polygons and reassigned to the one with which it shares the greatest perimeter.

    • Any polygon \(Q\) representing a gap between geometries \(P_1,\ldots, P_n\) is assigned to the one with which it shares the greatest perimeter.

    This function is probably sufficient when gaps and overlaps are all very small in area relative to the areas of the geometries, AND when the repaired file will only be used for operations like aggregating and prorating data. But it should NOT be relied upon when it is important for the repaired file to accurately represent adjacency relations between neighboring geometries, such as when a precinct shapefile is used as a basis for creating districting plans with contiguous districts.

    For instance, when a gap adjoins many geometries (which happens frequently along county boundaries in precinct shapefiles!), whichever geometry the gap is adjoined to becomes “adjacent” to all the other geometries adjoining the gap, which can lead to the creation of discontiguous districts in plans based on the repaired shapefile.

  2. smart_repair is a more sophisticated repair function designed to reproduce the “true” adjacency relations between geometries as accurately as possible. In the case of gaps that adjoin several geometries, this is accomplished by an algorithm that divides the gap into pieces to be assigned to different geometries instead of assigning the entire gap to a single geometry.

    In addition to repairing gaps and overlaps, smart_repair includes two optional features:

    • In many cases, the shapefile geometries are intended to nest cleanly into some larger units; e.g., in many states, precincts should nest cleanly into counties. smart_repair allows the user to optionally specify a second shapefile—e.g., a shapefile of county boundaries within a state—and then performs the repair process so that the repaired geometries nest cleanly into the units in the second shapefile.

    • Whether as a result of inaccurate boundaries in the original map or as an artifact of the repair algorithm, it may happen that some units share boundaries with very short perimeter but should actually be considered “queen adjacent”—i.e., intersecting at only a single point—rather than “rook adjacent”—i.e., intersecting along a boundary of positive length. smart_repair includes an optional step in which all rook adjacencies of length below a user-specified parameter are converted to queen adjacencies.

smart_repair can accept either a GeoSeries or GeoDataFrame as input, and the output type will be the same as the input type. The input must be projected to a non-geographic coordinate reference system (CRS)—i.e., not lat/long coordinates—in order to have sufficient precision for the repair. One option is to reproject a GeoDataFrame called gdf to a suitable UTM (Universal Transverse Mercator) projection via

gdf = gdf.to_crs(gdf.estimate_utm_crs())

At a minimum, all overlaps will be repaired in the output. Optional arguments include: * snapped (default value True): If True, all polygon vertices are snapped to a grid of size no more than \(10^{-10}\) times the maximum of width/height of the entire shapefile extent. HIGHLY RECOMMENDED to avoid topological exceptions due to rounding errors. * fill_gaps (default value True): If True, all simply connected gaps with area less than fill_gaps_threshold times the largest area of all geometries adjoining the gap are filled. Default threshold is \(0.1\); setting fill_gaps_threshold = None will fill all simply connected gaps. * nest_within_regions (default value None): If nest_within_regions is a secondary GeoSeries or GeoDataFrame of region boundaries (e.g., counties within a state) then the repair will be performed so that repaired geometries nest cleanly into the region boundaries; specifically, each repaired geometry will be contained in the region with which the original geometry has the largest area of intersection. Note that the CRS for the region GeoSeries/GeoDataFrame must be the same as that for the primary input. * min_rook_length (default value None): If min_rook_length is given a numerical value, all rook adjacencies with length below this value will be replaced with queen adjacencies. Note that this is an absolute value and not a relative value, so make sure that the value provided is in the correct units with respect to the input GeoSeries/GeoDataFrame’s CRS.