Hexagons, Spillovers, and the Augmentation Stack: Spatial Data Science for Compounding Knowledge

Why Most Enterprise Knowledge Architectures Collapse Three Dimensions Into a String Field — and the Geographic Layer That Fixes It

In the previous two articles I’ve argued for a four-layer Augmentation Stack — structured retrieval grounding, causal reasoning, persistent verifiable memory, and agentic verification loops — and a measurement framework that tells you whether the resulting corpus is actually compounding. Both pieces treated knowledge as if it lived in an undifferentiated semantic space. Documents have entities, entities have relationships, relationships have provenance.

That’s wrong, and it’s wrong in a way that costs organizations enormously.

A huge fraction of organizational knowledge is spatial. Customers live somewhere. Stores are sited somewhere. Supply chains move through somewhere. Patients receive care somewhere. Claims are filed from somewhere. Crops are grown somewhere. Risk events happen somewhere. Drilling permits are issued somewhere. School districts span somewhere. Property values exist somewhere. Election results aggregate somewhere. Most of the high-value decisions enterprises make have a “somewhere” component that is causally relevant, not just descriptive — and yet the typical enterprise knowledge architecture flattens that somewhere into a string field. city: "Kansas City". state: "MO". country: "US". As if geography were a label.

Place is not a label. Place is a structure — points, polygons, networks, hierarchies, neighborhoods, distances, flows, containment, adjacency, scale. The entire field of spatial data science exists because that structure is real, computable, and decision-relevant, and ignoring it loses information that no amount of fancier embeddings will recover. The geospatial foundation model literature of the last two years (Mai et al., 2024; Jakubik et al., 2023; Lu et al., 2024) makes this concrete in a new way: place can now be embedded as richly as language, and combined with language and structured knowledge in unified models. The architectural implication is that any serious knowledge engineering practice in 2026 needs a geographic layer that runs through all four layers of the Augmentation Stack.

This article walks through what that geographic layer looks like in practice. Schema choices, indexing structure, spatial retrieval, spatial causal inference, geospatial foundation models, verification loops, and a worked example pulling them together. Code is in Python. References are at the end.

The Flattening Problem

Walk through the data model of nearly any enterprise knowledge graph or data warehouse and you’ll find geography represented as one of three patterns, all bad.

Pattern one — geography as string. A location field containing free text like "Kansas City, MO" or "headquarters." Retrievable by exact match, occasionally normalized through geocoding, but otherwise inert. No topology, no hierarchy, no neighbors, no distances. Two records in the same neighborhood and two records on opposite sides of the country look identical to the system unless someone explicitly geocodes both and computes a distance.

Pattern two — geography as flat coordinate. A (lat, lon) pair on each record. Better — at least distance is computable — but still no hierarchy, still no topology, still no relationship to administrative geographies (census tracts, school districts, ZIP codes, congressional districts, market areas). Spatial joins require ad-hoc PostGIS queries that are written from scratch every time someone needs them.

Pattern three — geography as denormalized administrative geography. Each record gets enriched with state, county, MSA, ZIP, census tract — but as flat string fields, with no representation of the containment hierarchy, the adjacency graph, or the geometry. The system knows the record is in Jackson County, Missouri, but doesn’t know that Jackson County is adjacent to Cass and Clay, that it sits within the Kansas City MSA, that the MSA crosses two states, or that two records both in “Kansas City” might be in entirely different counties or even different states.

All three patterns have the same root cause: place has been treated as an attribute of the record rather than as a first-class entity in the knowledge model. The fix is not better geocoding. The fix is promoting place to a first-class citizen in the schema, with its own geometry, its own hierarchies, its own topology, and its own relationships — exactly as you’d treat customers, products, or any other domain entity. This is what the geospatial knowledge graph literature has been arguing for years (Janowicz et al., 2022, on KnowWhereGraph; Hu, 2018; Mai et al., 2022) and what the OGC GeoSPARQL standard formalizes for the semantic web. It’s been slow to reach mainstream enterprise practice because most data engineers were trained on relational or document models that have no native concept of spatial topology.

The cost of the flattening, in concrete terms: every spatial decision the organization makes has to be re-derived from scratch every time, because the spatial structure isn’t in the model. Site selection rebuilds the catchment analysis from raw data. Risk scoring re-imports the flood zones. Marketing rebuilds the trade areas. Public health rebuilds the SDOH overlays. The spatial knowledge that should compound across the organization gets lost because the substrate has no place to keep it.

Place as a First-Class Entity

The starting move is to model place explicitly. A spatial knowledge graph treats places as nodes with their own identity, geometry, properties, and relationships. The minimum useful schema:

from dataclasses import dataclass, field
from typing import Optional
from shapely.geometry import shape  # OGC geometry types

@dataclass
class Place:
    id: str                          # canonical ID (e.g., GEOID, OSM ID, custom)
    name: str
    geometry: object                 # Shapely geometry: Point, Polygon, MultiPolygon
    place_type: str                  # 'parcel' | 'block' | 'tract' | 'county' | ...
    parent_id: Optional[str] = None  # containment hierarchy
    valid_from: Optional[str] = None # temporal validity start
    valid_to: Optional[str] = None   # temporal validity end (boundaries change)
    properties: dict = field(default_factory=dict)

@dataclass
class SpatialRelation:
    subject_id: str
    predicate: str        # 'contains' | 'within' | 'touches' | 'intersects' |
                          # 'within_distance' | 'connected_via' | ...
    object_id: str
    metric: Optional[float] = None  # e.g., distance in meters, shared border length

This is closer to the OGC Simple Features model and GeoSPARQL relations than to a typical property graph schema, and that’s intentional. The DE-9IM topological relations (contains, within, touches, crosses, overlaps, disjoint, equals, intersects, covers) give you the full algebra of spatial relationships, all of which are computable from geometry and all of which carry knowledge a string field doesn't.

Three structural relationships matter most in practice:

Containment hierarchy. Parcels nest inside blocks, blocks inside tracts, tracts inside counties, counties inside states, states inside nations. Many other hierarchies overlay this — ZIP Code Tabulation Areas, congressional districts, school districts, hospital service areas, designated market areas, watersheds, ecoregions. A useful spatial KG indexes multiple hierarchies and supports queries like “which administrative geographies of any type contain this point?”

Adjacency / topology. Which places share a border with which other places. Adjacency is what lets you compute spatial spillover effects — if a county changes a policy, its neighbors are the natural comparison group, and “neighbor” is a graph relation, not a metric one. Queen contiguity (touching at a vertex or edge) and rook contiguity (touching along an edge) are the standard definitions; libpysal computes both from any polygon dataset.

Network connectivity. Roads, transit lines, power grids, supply chains, watersheds, airline routes. Spatial proximity is not the same as network proximity — two addresses 200 meters apart on opposite sides of a freeway with no crossing are spatially close and network-distant, and the network distance is what determines real accessibility for most decisions.

A spatial entity resolution pass that maps records into this structure is the equivalent of named entity resolution in the linguistic case. You take whatever location reference appears in the source data — a free-text address, a (lat, lon), a place name, a ZIP code — and resolve it to the canonical place IDs that already exist in your spatial KG. The standard tooling is good enough now that this is engineering, not research: geopandas plus shapely for geometry, geopy or a dedicated geocoder (Pelias, Nominatim, Mapbox, Google) for free-text resolution, geopy.distance or pyproj for projections, and libpysal.weights for adjacency graphs.

import geopandas as gpd
from shapely.geometry import Point
import libpysal as ps

class SpatialResolver:
    def __init__(self, place_layers: dict[str, gpd.GeoDataFrame]):
        """
        place_layers: dict like {'tracts': gdf_tracts, 'counties': gdf_counties, ...}
        Each layer must be in the same CRS (use a projected CRS for distance work).
        """
        self.layers = place_layers
        # Pre-build spatial indices for fast point-in-polygon queries
        for name, gdf in self.layers.items():
            gdf.sindex  # builds R-tree

    def resolve_point(self, lon: float, lat: float, crs="EPSG:4326") -> dict:
        """Return all containing places across hierarchies for a single point."""
        pt = gpd.GeoSeries([Point(lon, lat)], crs=crs)
        out = {}
        for name, gdf in self.layers.items():
            pt_proj = pt.to_crs(gdf.crs)
            # spatial join via sindex
            possible = list(gdf.sindex.query(pt_proj.iloc[0], predicate="within"))
            for idx in possible:
                if pt_proj.iloc[0].within(gdf.geometry.iloc[idx]):
                    out[name] = gdf.iloc[idx]["GEOID"]
                    break
        return out

    def adjacency(self, layer_name: str, contiguity: str = "queen"):
        """Build the adjacency weights matrix for a polygon layer."""
        gdf = self.layers[layer_name]
        if contiguity == "queen":
            w = ps.weights.Queen.from_dataframe(gdf, idVariable="GEOID")
        elif contiguity == "rook":
            w = ps.weights.Rook.from_dataframe(gdf, idVariable="GEOID")
        else:
            raise ValueError(contiguity)
        w.transform = "r"  # row-standardize
        return w

The output of resolution is a record enriched with all of its place memberships across hierarchies, ready to be joined to anything else in the corpus that shares any of those place IDs. This is the spatial analog of entity coverage from the measurement framework — and the same metric applies. Spatial entity coverage is the fraction of geo-bearing records whose locations resolve cleanly into the canonical place hierarchy. Below 0.85 and your spatial joins will silently miss large portions of the corpus.

Spatial Indexing at Scale

The naive way to do spatial joins — point-in-polygon, nearest-neighbor, distance buffers — scales badly. R-trees help up to a few hundred million features. Beyond that, and especially in distributed compute environments, you need a hierarchical spatial indexing scheme that lets you partition the world into cells that can be sharded, joined, and aggregated without re-computing geometry every time.

The two production-grade systems worth knowing are H3 (Uber’s hexagonal hierarchical spatial index) and S2 (Google’s). Both tile the globe at multiple resolutions, both give every cell a canonical 64-bit integer ID, and both support fast neighbor lookup, containment, and resolution change. They differ in cell shape — H3 uses hexagons, S2 uses spherical quadrilaterals — and in the precise hierarchical structure. Hexagons have a single, isotropic neighbor distance (six neighbors, all equidistant) which makes them generally preferable for analytics; S2’s quad structure is preferable when you need clean nesting at every level.

The practical advantages of moving spatial work onto an H3 grid are large. Spatial joins become JOIN ON h3_cell_id — a string equality, fast at any scale. Aggregations become GROUP BY h3_cell_id. Neighborhoods become h3.grid_disk(cell, k) — pure integer math, no geometry computation. Visualizations become hex maps, which avoid the visual lying that goes with arbitrary administrative polygons. Most importantly, an H3 grid is a stable, deterministic substrate on which you can co-register many heterogeneous datasets that arrive with different native geometries — points, polygons, rasters, networks — and have them all index into the same hex IDs at the same resolution.

import h3
import pandas as pd
import geopandas as gpd

def to_h3(gdf: gpd.GeoDataFrame, resolution: int = 9) -> pd.DataFrame:
    """
    Map a GeoDataFrame to H3 cells at the given resolution.
    Resolution 9 ≈ 0.1 km² per cell — appropriate for neighborhood-level analytics.
    See https://h3geo.org/docs/core-library/restable for a full resolution table.
    """
    gdf_wgs = gdf.to_crs("EPSG:4326")
    rows = []
    for idx, row in gdf_wgs.iterrows():
        geom = row.geometry
        if geom.geom_type == "Point":
            cells = [h3.latlng_to_cell(geom.y, geom.x, resolution)]
        elif geom.geom_type in ("Polygon", "MultiPolygon"):
            cells = list(h3.geo_to_cells(geom.__geo_interface__, resolution))
        else:
            continue
        for c in cells:
            r = row.drop("geometry").to_dict()
            r["h3_cell"] = c
            rows.append(r)
    return pd.DataFrame(rows)

def neighbors_within_k(cell: str, k: int = 2) -> list[str]:
    """All cells within k rings of the given cell (including the cell itself)."""
    return list(h3.grid_disk(cell, k))

Resolution choice matters. H3 resolution 7 is roughly 5 km² per hex (county-sized analytics), 8 is roughly 0.7 km² (neighborhood-scale), 9 is 0.1 km² (block-scale), 10 is 15,000 m² (parcel-scale). The right resolution depends on the modifiable areal unit problem, which I’ll come to in the spatial causal inference section — TL;DR, the choice of analytical unit is itself a substantive modeling decision that should be made deliberately, not by accident of the data’s native geometry.

The DuckDB spatial extension (Holub et al., 2024, in the DuckDB project) and BigQuery’s GIS support both have native H3 functions now, which means you can do spatial work at warehouse scale without leaving SQL. PostGIS remains the most feature-complete spatial database; the trade-off is operational footprint.

Spatial-Aware Retrieval

Standard RAG retrieval is dense-vector similarity over text chunks, sometimes hybridized with BM25. It is geographically blind. Ask “what’s the average household income in the neighborhoods around our Liberty store?” and a vanilla retriever will surface chunks that talk about Liberty stores and chunks that talk about household income, with no understanding that around is a spatial operator that should constrain retrieval to chunks geographically near a specific point. The result is text that sounds plausible and is regularly wrong.

Spatial-aware retrieval treats geographic constraints as first-class filters and rankers. The pattern that works in practice has three components.

Geographic intent extraction. Parse the query for spatial operators (“near,” “within,” “across,” “between,” “from X to Y”) and reference points (place names, addresses, coordinates), and extract them into a structured spatial filter. An LLM with a structured output schema does this well; the output is a dict with geographic_constraint, reference_geometry, and distance_or_topology fields.

Spatial pre-filter. Apply the constraint as a hard filter on the candidate corpus before semantic ranking. If the query specifies “within 5 km of the Liberty store,” the retrieval pool shrinks to chunks whose source documents have a place ID resolving inside that 5 km buffer. If the query specifies “in the Kansas City MSA,” the pool shrinks to chunks tied to places contained within the MSA’s geometry. The H3 grid makes this operationally cheap — pre-compute every chunk’s H3 cells, and the spatial filter becomes a set intersection.

Geographically-weighted re-ranking. For chunks that pass the hard filter, the final ranker can blend semantic similarity with geographic centrality: chunks closer to the reference geometry get a boost. The functional form that’s worked well for me is a multiplicative reranker with a spatial decay term:

score(c, q)  =  semantic_sim(c, q)  ·  spatial_decay( d(c, q) )

where:
   semantic_sim  =  cross-encoder or cosine similarity score
   d(c, q)       =  geographic distance from chunk c to query reference q
                    (or 0 if c is inside the reference geometry)
   spatial_decay = exp( - d / σ ),  σ chosen for the query type
                   (e.g., σ = 1 km for hyperlocal, 50 km for regional)

This is the Tobler-flavored extension of standard RAG ranking and it dramatically improves hyperlocal query quality. The downside is that you need clean geographic metadata on every retrievable chunk — which is exactly the schema discipline the spatial KG section was about. The work compounds across applications, though, and the same metadata supports site selection, risk modeling, and a dozen other spatial workloads.

Spatial Causal Inference: The Spillover Problem

The measurement article walked through standard causal inference methods — difference-in-differences, BSTS, synthetic control, instrumental variables. Every one of those methods has a spatial extension, and ignoring the spatial dimension is one of the most common ways causal claims fail when scrutinized.

Three issues matter most.

Spatial autocorrelation. Tobler’s First Law of Geography: everything is related to everything else, but near things are more related than distant things. Empirically this means observations close in space have correlated outcomes — neighboring counties have correlated unemployment, adjacent neighborhoods have correlated home prices, contiguous farms have correlated yields. Standard regression models assume independent residuals, and when the residuals are spatially autocorrelated, the standard errors are wrong, often by a lot. The fix is spatial econometric models — spatial lag (SAR), spatial error (SEM), and spatial Durbin (SDM) — that explicitly model the dependency. Anselin (1988) and LeSage & Pace (2009) are the canonical references; the spreg package in PySAL is the production-grade Python implementation.

The diagnostic is Moran’s I (Moran, 1950), which tests whether residuals show spatial structure beyond what would be expected by chance:

from libpysal.weights import Queen
from esda.moran import Moran
import geopandas as gpd

def diagnose_spatial_autocorrelation(gdf: gpd.GeoDataFrame,
                                      residual_col: str) -> dict:
    """Moran's I on a residual series; tells you whether OLS is invalid."""
    w = Queen.from_dataframe(gdf)
    w.transform = "r"
    moran = Moran(gdf[residual_col].values, w, permutations=999)
    return {
        "I": moran.I,
        "expected_I": moran.EI,
        "p_value_sim": moran.p_sim,    # permutation-based, robust to non-normality
        "z_score": moran.z_sim,
    }

A Moran’s I significantly different from its expected value (typically very close to zero) is the empirical signature of spatial dependence in your residuals. If your causal lift estimate from the previous article was computed without checking for it, the standard error is suspect.

Spatial spillovers. A treatment in one location often affects outcomes in neighboring locations. A new highway in one county changes traffic patterns in adjacent counties. A minimum wage increase in one state changes labor markets in border counties of neighboring states. A Medicaid expansion in one state affects hospital utilization across the state line. The classic difference-in-differences design assumes the control group is unaffected by treatment of the treated group — a SUTVA assumption. Spatial spillovers violate this directly. The remediation is spatial DiD with spillover terms (Verbitsky-Savitz & Raudenbush, 2012; Delgado & Florax, 2015) or the more recent spatial-econometric approaches (Reich et al., 2021; Akbari et al., 2023). When the spillovers are large relative to the direct effect, ignoring them produces biased estimates of both.

Modifiable areal unit problem (MAUP). The same underlying spatial data, aggregated to different units (tracts vs. counties vs. ZIPs vs. arbitrary hex grids), can produce different and even contradictory statistical results (Openshaw, 1984; Fotheringham & Wong, 1991). This is not a quirk; it’s a fundamental property of spatial data and one of the strongest arguments for working in stable hierarchical grids like H3 rather than administrative boundaries when the analysis allows it. When the analysis requires administrative units (because the policy intervention or business unit is administrative), the MAUP should be addressed by sensitivity analysis across multiple aggregation schemes — if the conclusion holds across census tracts, ZCTAs, and an H3 r8 grid, you have a more defensible claim than if it only holds at one aggregation.

A spatial DiD specification looks like:

Y_it  =  α_i  +  λ_t  +  β · D_it  +  ρ · (W · D)_it  +  X_it · γ  +  ε_it

where:
   D_it       =  treatment indicator for unit i at time t
   W          =  spatial weights matrix (e.g., row-standardized queen contiguity)
   (W · D)_it =  fraction of i's neighbors that are treated at time t
   β          =  direct treatment effect
   ρ          =  spillover effect on neighbors

The β + ρ partition is what tells you the total treatment effect, which is what most policy questions actually care about. Reporting only β is the spatial-econometric analog of cherry-picking a subgroup analysis. If you’re doing causal lift estimates for an organizational intervention that has any spatial structure — a region-by-region rollout, a store-by-store deployment — you should be reporting both.

Geospatial Foundation Models: The New Layer

The most consequential development in spatial data science over the last two years is the arrival of geospatial foundation models that embed satellite imagery, location, and spatial context as richly as language models embed text. This is the layer that lets you bring continuous-space, raster-derived, remotely-sensed knowledge into the same architecture as the discrete-space, document-derived, semantic knowledge that the rest of the Augmentation Stack handles.

The four checkpoints worth knowing as of this writing:

Prithvi (NASA-IBM, 2023; Jakubik et al.). A masked-autoencoder-pretrained ViT trained on Harmonized Landsat-Sentinel imagery, designed for downstream geospatial tasks like flood mapping, burn-scar detection, crop type classification, and land use. Open-weights via Hugging Face. The 100M-parameter version is the practical default; a 600M version is also available. Strong on multispectral and temporal signals.

SatMAE (Cong et al., 2022, NeurIPS). Earlier MAE-based pretraining on satellite imagery, with explicit handling of multi-spectral bands and temporal sequences. Still a useful comparison baseline and good for cases where the Prithvi pretraining domain doesn’t overlap your application.

Clay Foundation Model (2024, open source). A more recent, fully-open geospatial foundation model with a strong embedding API and good documentation for downstream fine-tuning. The Clay project’s positioning as a community-governed open model makes it a defensible choice for organizations with open-source compliance constraints.

SpectralGPT (Hong et al., 2024, IEEE TPAMI). Spectral-band-aware foundation model handling hyperspectral and multispectral data more rigorously than natural-image-pretrained backbones.

Plus the survey-level work that gives the broader landscape: Mai et al. (2024) in ACM Computing Surveys on foundation models for GeoAI, Lu et al. (2024) on AI foundation models for Earth observation, and Bountos et al. (2025) on vision foundation models for remote sensing.

The architectural pattern that brings these into the Augmentation Stack: use a geospatial foundation model to embed every cell of your H3 grid (or every parcel, tract, hex, whatever your spatial unit is) using whatever imagery and spatial context is available — recent satellite imagery, historical imagery, weather, terrain, road network density, points of interest. The output is a dense vector per spatial unit that captures place character in the same way a language embedding captures text character. These embeddings are queryable by similarity (“what other neighborhoods in our region are most similar to this target neighborhood?”), composable with semantic embeddings (“find document chunks about retail performance that are anchored in places similar to our top-performing stores”), and stable enough to be cached and reused across many downstream applications.

# Sketch of using a geospatial foundation model for place embeddings.
# The exact API differs by checkpoint; this shows the pattern with Clay.
import torch
from clay_model import ClayMAE  # illustrative; check current Clay package

def embed_places(
    place_grid: gpd.GeoDataFrame,        # one row per spatial unit (e.g., H3 cell polygon)
    imagery_loader,                       # function: place -> tensor of imagery patches
    model: ClayMAE,
    device: str = "cuda",
) -> dict[str, torch.Tensor]:
    """Compute a foundation-model embedding for each place."""
    model = model.to(device).eval()
    embeddings = {}
    with torch.no_grad():
        for _, row in place_grid.iterrows():
            pid = row["place_id"]
            patch = imagery_loader(row).to(device)   # [B, C, H, W]
            emb = model.encode(patch).mean(dim=0)    # pool to a single vector
            embeddings[pid] = emb.cpu()
    return embeddings

The integration into the broader Augmentation Stack is the interesting part. Layer 1 retrieval gets richer because you can now find geographically-relevant places by similarity, not just containment or distance. Layer 2 causal modeling gets richer because place-character embeddings are valid covariates for adjustment, capturing dimensions of place that no list of administrative attributes ever fully captures. Layer 3 memory persists per-place embeddings as a stable substrate that updates over time as new imagery and context arrive — a place’s embedding trajectory tells you whether it’s changing. Layer 4 verification gets richer because remote sensing provides ground truth for many spatial claims that previously had to be taken on faith from a single data source.

Spatial Verification Loops

The fourth layer of the Augmentation Stack — agentic verification — has a direct spatial extension. Many claims an AI system might produce about the world are spatially verifiable in ways that purely textual claims aren’t. “There’s a new building at this location” can be checked against recent imagery. “This farm shifted from corn to soybeans” is visible in NDVI time series. “This neighborhood gained a transit stop” shows up in OSM diffs. “This region experienced flooding” is observable in Sentinel-1 SAR backscatter changes.

The pattern is: when an agent produces a claim that has a spatial referent, route the claim to a spatial verifier that checks it against the appropriate remote-sensing or geographic data source. The verifier doesn’t need to be elaborate — a small set of well-chosen checks (recent imagery comparison, OSM lookup, NDVI delta, building-footprint database query) covers a surprising fraction of the spatial claims that come up in business workloads. What matters is that the verification loop is real and the validated outputs are what get retained.

This is the same FunSearch-style propose-execute-verify-retain pattern from the architecture article, applied in the spatial domain. The difference from the textual case is that the verifier is grounded in physical, sensed reality, which is the strongest possible kind of fitness function.

A Worked Example: Healthcare Access in Practice

Let me make this concrete with a sketch of how the layers come together for one application.

Suppose the question is: which neighborhoods in our service area have the worst combination of healthcare access barriers and demographic vulnerability, and how does that map to the locations of our outreach programs?

A naive answer geocodes the program addresses, joins to ZIP-level demographic data, sorts by some composite score, and produces a list. It’s wrong in five ways: ZIP is the wrong unit (ZCTAs are administrative postal areas, not coherent communities); the demographic data is correlational, not causal; the access measure is straight-line distance, not network distance; there’s no validation against ground truth; and the corpus the answer comes from doesn’t compound — next quarter’s analysis rebuilds the whole thing from scratch.

The spatial-aware version:

Layer 1 (spatial KG + indexing). Resolve every program address to a canonical place ID. Build a KG of clinics, transit stops, road network, demographic units (census block groups), and administrative geographies (county, MSA, HSA, congressional district), all indexed to H3 r9 and linked through containment, adjacency, and network connectivity edges. Pull SDOH data from AHRQ, demographics from ACS, healthcare facility data from HRSA and CMS, transit from GTFS feeds. Every record is co-registered to the H3 grid.

Layer 2 (spatial causal modeling). Build a spatial autoregressive model of healthcare utilization as a function of access (network distance to nearest provider of each type), demographic composition, insurance status, and a place-character embedding from Prithvi or Clay. Estimate spillover effects of clinic openings on neighboring tracts. Use modern DiD methods (Callaway-Sant’Anna) to estimate the causal effect of past program rollouts on utilization, with explicit spatial spillover terms. Run sensitivity analyses across multiple aggregation schemes (block group, tract, ZCTA, H3 r8) to bound the MAUP risk.

Layer 3 (spatio-temporal memory). Persist every analytic decision — the units chosen, the model specification, the estimated effects, the neighborhoods identified — with full provenance back to source data, methodology, and model version. Update the place-character embeddings quarterly as new imagery arrives. Track which program locations were recommended by past analyses and what the realized outcomes were.

Layer 4 (spatial verification). When the system claims a neighborhood has a particular access barrier, validate it against recent imagery (is there construction blocking the road on this OSM segment?), against transit feeds (is the bus route still active?), against facility data (is the clinic still open and accepting new patients?). Where the system produces a recommendation, run a placebo test — does the model assign similar priority scores to placebo neighborhoods that are demographically similar but have no real access barriers? Calibrate the recommendation confidences against historical hit rates.

The output is not a static list. It’s a continuously-updated, causally-grounded, spatially-rigorous, verified ranking with explicit confidence intervals and known sensitivities. A board-level chart shows the compounding metric: every quarter, the corpus of place-knowledge grows, the causal estimates tighten as more outcomes accumulate, the verifier disagreement rate falls as resolution improves, and the program targeting decisions get measurably better.

This is what happens when geography is a first-class citizen of the knowledge architecture rather than a string field.

The Spatial Stack in Tools

A defensible spatial knowledge engineering stack in 2026, from substrate to deployment:

Geometry and projections. Shapely 2.0 (Gillies et al.), GeoPandas (Jordahl et al.), pyproj (Snow et al.), Fiona / GDAL.

Spatial indexing. H3 (Uber, h3-py), S2 (Google, s2sphere), DuckDB spatial extension, PostGIS.

Spatial econometrics. PySAL ecosystem — libpysal (weights), esda (Moran's I, LISA, Geary's C), spreg (spatial regression), mgwr (geographically weighted regression), pointpats (point pattern analysis), segregation (spatial segregation indices).

Spatial machine learning. scikit-learn for baselines; pytorch-geometric and pytorch-geometric-temporal for spatio-temporal GNNs; torchgeo and terratorch for remote-sensing model fine-tuning.

Geospatial foundation models. Prithvi via Hugging Face; Clay Foundation Model from the Clay project; SatMAE from the original repo; SpectralGPT for spectral-aware tasks.

Raster and imagery handling. rasterio, xarray, rioxarray, dask for large-scale raster, stackstac and pystac-client for STAC catalogs, planetary-computer for Microsoft's hosted catalog.

Network analysis. osmnx for OpenStreetMap-derived street networks, networkx for general graph algorithms, r5py and valhalla / OpenTripPlanner for transit and multimodal routing.

Geocoding. Pelias (open source, self-hostable), Nominatim (OSM-backed), commercial APIs (Mapbox, Google, Esri) when accuracy and coverage matter more than cost.

Visualization. lonboard for high-performance interactive maps in notebooks, kepler.gl for exploratory analytics, pydeck for production deck.gl integration, geemap for Earth Engine integration, folium for simple Leaflet output.

Spatial causal inference. PySAL’s spreg for SAR/SEM/SDM models, differences and did for modern DiD, custom code for explicit spatial spillover specifications. The R ecosystem (spdep, spatialreg, geex) is more mature for cutting-edge spatial-causal methods; Python is catching up.

Pipeline glue. dbt for transformation, prefect or dagster for orchestration, delta or iceberg for storage, dvc or lakefs for data versioning. None of these are spatial-specific, but the spatial work has to live inside whatever data engineering substrate the rest of the organization uses.

Where to Start

If you’re bringing spatial knowledge engineering into an organization that doesn’t have it, the order of operations:

Start with the schema move. Promote place to a first-class entity. Build the canonical place hierarchy from public data (TIGER/Line for the US administrative geography is free and authoritative; OSM for finer-grained features; STAC catalogs for imagery). This is a few weeks of work and produces a substrate every other spatial workload will reuse for years.

Add H3 indexing as a co-registration layer across all existing geo-bearing data. This is mostly an ETL pass that adds an h3_cell column at one or two resolutions to every relevant table. It pays off the first time someone asks a spatial question that previously required a custom GIS analyst engagement.

Layer in spatial entity resolution. The metric to track is spatial entity coverage — what fraction of geo-bearing records resolve cleanly into the canonical hierarchy. Drive it above 0.85 for canonical entities.

Add spatial diagnostics to your existing causal models. Moran’s I on the residuals is a one-line addition that often reveals that existing causal claims are quietly violating independence assumptions. Where it matters, escalate to spatial econometric models.

Bring in geospatial foundation model embeddings as features and as similarity primitives. Prithvi or Clay weights on Hugging Face, applied per H3 cell, produce embeddings that drop into existing ML pipelines as additional covariates with substantial lift in many use cases (Mai et al., 2024, surveys the empirical evidence).

Close the loop with spatial verification. Pick three or four spatial claim types your system produces routinely, build a verifier for each, and start retaining only validated outputs into the long-term corpus.

The full stack is a multi-quarter build. The first phase — schema, indexing, entity resolution — pays for itself within two quarters in any organization with non-trivial spatial workloads.

The Through-Line

The argument across this three-article series has been a single one. Knowledge in organizations compounds when the substrate is structurally sound, the reasoning over it is causally honest, the memory of decisions is persistent and verifiable, and the verification loop is grounded in external reality. The Augmentation Stack is the architecture; the measurement framework is the instrumentation; and the spatial layer is the dimension that turns a flat knowledge architecture into one that respects how the world is actually shaped.

The world is three-dimensional, time-varying, networked, and hierarchical at multiple scales. Knowledge architectures that flatten any of those dimensions lose information in ways that cost real money. The spatial flattening is the most common, the most costly in spatially-driven industries (insurance, healthcare, retail, logistics, real estate, agriculture, energy, public sector), and the most fixable with current open-source tooling.

Geography is knowledge. Place is structure. The substrate evolves under selection. The selection layer is the architecture. The architecture has a geographic dimension, and the organizations that build it that way are the ones whose knowledge corpus actually reflects the world it’s supposed to model.

That’s the job.

References

Geospatial foundation models and GeoAI

Bountos, N. et al. (2025). Vision Foundation Models in Remote Sensing: A Survey. IEEE Geoscience and Remote Sensing Magazine.
Cong, Y. et al. (2022). SatMAE: Pre-training Transformers for Temporal and Multi-Spectral Satellite Imagery. NeurIPS.
Hong, D. et al. (2024). SpectralGPT: Spectral Foundation Model. IEEE TPAMI.
Jakubik, J. et al. (2023). Foundation Models for Generalist Geospatial Artificial Intelligence (Prithvi). arXiv.
Lu, S. et al. (2024). AI Foundation Models for Earth Observation. IEEE Geoscience and Remote Sensing Magazine.
Mai, G. et al. (2022). Towards a Foundation Model for Geospatial Artificial Intelligence. SIGSPATIAL.
Mai, G. et al. (2024). On the Opportunities and Challenges of Foundation Models for Geospatial Artificial Intelligence. ACM Computing Surveys.
Reed, C.J. et al. (2023). Scale-MAE: A Scale-Aware Masked Autoencoder for Multiscale Geospatial Representation Learning. ICCV.

Geospatial knowledge graphs and semantic GIS

Hu, Y. (2018). Geo-text Data and Data-Driven Geospatial Semantics. Geography Compass.
Janowicz, K. et al. (2022). Know, Know Where, KnowWhereGraph: A Densely Connected, Cross-Domain Knowledge Graph and Geo-Enrichment Service. AI Magazine.
OGC GeoSPARQL — A Geographic Query Language for RDF Data, OGC Standard 11–052r4.

LLMs and geographic reasoning

Bhandari, P. et al. (2023). Are Large Language Models Geospatially Knowledgeable? SIGSPATIAL.
Kuckreja, K. et al. (2024). GeoChat: Grounded Large Vision-Language Model for Remote Sensing. CVPR.
Mooney, P. et al. (2023). ChatGPT and OpenStreetMap. Transactions in GIS.
Roberts, J. et al. (2023). GPT4GEO: How a Language Model Sees the World’s Geography. arXiv.

Spatial econometrics and statistics

Anselin, L. (1988). Spatial Econometrics: Methods and Models. Kluwer.
Anselin, L. (2010). Thirty Years of Spatial Econometrics. Papers in Regional Science.
Fotheringham, A.S. & Wong, D. (1991). The Modifiable Areal Unit Problem in Multivariate Statistical Analysis. Environment and Planning A.
LeSage, J. & Pace, R.K. (2009). Introduction to Spatial Econometrics. CRC Press.
Moran, P.A.P. (1950). Notes on Continuous Stochastic Phenomena. Biometrika.
Openshaw, S. (1984). The Modifiable Areal Unit Problem. CATMOG 38.
Tobler, W. (1970). A Computer Movie Simulating Urban Growth in the Detroit Region. Economic Geography.

Spatial causal inference

Akbari, K. et al. (2023). Spatial Causal Inference in the Presence of Unmeasured Confounding. arXiv.
Delgado, M.S. & Florax, R.J.G.M. (2015). Difference-in-Differences Techniques for Spatial Data: Local Autocorrelation and Spatial Interaction. Economics Letters.
Papadogeorgou, G. et al. (2022). Causal Inference with Spatio-Temporal Data: Estimating the Effects of Airstrikes on Insurgent Violence in Iraq. Journal of the Royal Statistical Society B.
Reich, B.J. et al. (2021). A Review of Spatial Causal Inference Methods for Environmental and Epidemiological Applications. International Statistical Review.
Verbitsky-Savitz, N. & Raudenbush, S.W. (2012). Causal Inference Under Interference in Spatial Settings. Statistics, Politics, and Policy.

Spatial indexing

Brodsky, I. (2018). H3: Uber’s Hexagonal Hierarchical Spatial Index. Uber Engineering.
Sahr, K. et al. (2003). Geodesic Discrete Global Grid Systems. Cartography and Geographic Information Science.

Remote sensing and machine learning

Rolf, E. et al. (2021). A Generalizable and Accessible Approach to Machine Learning with Global Satellite Imagery (MOSAIKS). Nature Communications.
Stewart, A.J. et al. (2022). TorchGeo: Deep Learning with Geospatial Data. SIGSPATIAL.

Modern DiD (carried from prior article for spatial extensions)

Callaway, B. & Sant’Anna, P. (2021). Difference-in-Differences with Multiple Time Periods. Journal of Econometrics.
Goodman-Bacon, A. (2021). Difference-in-Differences with Variation in Treatment Timing. Journal of Econometrics.
Roth, J. et al. (2023). What’s Trending in Difference-in-Differences? Journal of Econometrics.

Tools and software

Gillies, S. et al. Shapely Documentation. https://shapely.readthedocs.io
Jordahl, K. et al. GeoPandas Documentation. https://geopandas.org
PySAL Project. https://pysal.org
Holub, M. et al. (2024). DuckDB Spatial Extension. https://duckdb.org/docs/extensions/spatial

Brian Curry is an applied AI researcher and the founder of Vector1 Research, an independent lab focused on knowledge engineering, causal inference, agent evaluation, and the architecture of compounding organizational intelligence. He is the author of open-source frameworks including MeaningFlow (semantic content modeling), Papilon (causal inference and complex systems simulation), and Memory-Node Encapsulation (autobiographical memory architecture for agentic systems). He writes regularly on applied AI, knowledge architecture, and decision science.

Connect on Medium · GitHub · LinkedIn