Capstone: A GIS Overlay and Spatial Query Engine

The Project

A GIS system answers a small, recurring set of questions about map layers — parcels, zoning districts, flood plains, roads — that were each built independently and now need to be combined and queried together: "which zoning district is this parcel in?" (overlay), "what's inside this map viewport?" (range query), and "which parcel did the user just click on?" (point location). None of these questions needs a new algorithm. Every one of them is a direct application of a part already covered in this series, and the capstone's actual deliverable is the shared representation that lets all three questions reuse the same underlying structures instead of re-deriving them per query.

Architecture

The engine has exactly four stages, each one a direct reuse of an earlier part:

1. Ingest. Parse each layer's raw polygons into a DCEL (Part 7), normalizing orientation and rejecting or repairing self-intersections using the robust orientation and segment-intersection predicates from Parts 1 and 6.
2. Overlay. To combine two layers, sweep-line-node their boundaries together (Part 6), merge the results into one DCEL, propagate inside/outside labels for each source layer across every face, and select faces per the requested Boolean operation (Part 8) — union for "either layer," intersection for "both," and so on.
3. Index. Build a 2D range tree over each layer's face bounding boxes for window queries (Part 11), and a randomized incremental trapezoidal map over the same DCEL's edges for point-location queries (Part 15).
4. Query. A window query walks the range tree's canonical subtrees; a point-location query walks the trapezoidal map's history DAG to a leaf, then reads that leaf's face label straight off the DCEL.

class SpatialLayer:
    """One GIS layer: a DCEL plus the two indexes built on top of it."""

    def __init__(self, polygons):
        # Part 1's orientation predicate normalizes every ring's winding
        # before Part 7's DCEL builder assigns half-edges and faces.
        self.dcel = build_dcel(normalize_rings(polygons))
        self.range_index = build_range_tree(face_bounding_boxes(self.dcel))
        self.point_index = build_trapezoidal_map(self.dcel)

    @classmethod
    def overlay(cls, layer_a, layer_b, operation):
        """operation in {'union', 'intersection', 'difference', 'xor'}."""
        # Part 6: node both layers' edges at every crossing/touch/overlap.
        noded = sweep_line_node(layer_a.dcel, layer_b.dcel)
        # Part 8: merge into one DCEL, propagate labels, select faces.
        merged_dcel = merge_and_label(noded, layer_a, layer_b)
        kept_faces = select_faces(merged_dcel, operation)
        return cls.from_dcel(trace_boundaries(kept_faces))

    def query_window(self, rect):
        # Part 11: canonical-subtree decomposition of the query rectangle.
        return self.range_index.report(rect)

    def locate(self, point):
        # Part 15: history-DAG walk to the containing trapezoid's face.
        return self.point_index.locate(point).face_label

Why This Is the Capstone, Not Just a Recap

Building each piece in isolation is what Parts 1–15 already did. The capstone's real design problem is keeping all three query paths consistent with a single evolving DCEL: an overlay operation replaces faces, which means the range tree's bounding boxes and the trapezoidal map's history DAG both need to reflect the new topology before the next query runs. Two disciplined choices make this tractable:

• Rebuild indexes, don't patch them, after every overlay. Both the range tree (Part 11) and the trapezoidal map (Part 15) are built for static input; incremental repair of either after an overlay is far more intricate than a O(n log n) rebuild from the merged DCEL's faces, and in a GIS pipeline overlay operations are batch steps, not per-click events.
• Never let two different layers disagree about a boundary that both claim to share. A parcel boundary that coincides with a zoning-district boundary must be noded into the same vertex and edge chain in the merged DCEL (Part 6's shared-endpoint and collinear-overlap policies), or the overlay's face labeling silently produces sliver faces along every shared edge — the single most common visible bug in real GIS overlay tools.

What to Build

A minimal version of this capstone needs only three real GIS layers (say, parcels, zoning, and a flood-plain polygon), one overlay operation, one window query, and one point-location query — enough to prove that all three query types genuinely share the ingest and overlay stages rather than reimplementing geometry per query type. The full system this capstone sketches — multiple simultaneous layers, incremental updates, and a real tile-serving front end — is a natural next step once the shared-representation discipline above is in place.