What is a spatial join in GeoPandas?

A spatial join matches rows from two geographic datasets by their location relationship, such as points falling within polygons, rather than by a shared key column.

How do we avoid wrong results in a spatial join?

A basic safeguard is to confirm both layers use a compatible coordinate reference system, inspect unmatched records, and test points near polygon boundaries before relying on the join.

Is the article's spatial analysis code publicly reproducible?

No. As of July 21, 2026, the linked public repository does not include the promised implementation or spatial validation script. The code on this page is instructional, and the empirical claims are article-reported rather than publicly reproduced.

GeoPandas Spatial Joins: An Illustrative Workflow

The foundation of many spatial analyses is a spatial join: merging point data, such as provider locations, with polygon boundaries, such as census tracts, using gpd.sjoin(). Once the data are joined, we can build neighborhood weights and test for spatial patterns. But the workflow can still produce plausible-looking output when coordinate reference systems are incompatible or geometries are missing. Spatial dependence can also violate the independence assumptions of downstream statistical tests.

This walkthrough covers loading boundaries, computing centroids, constructing spatial weights, testing for autocorrelation with Moran's I, and identifying local clusters via LISA. It uses an article-reported example from a food security mobility study of a California county, described here as containing 108 census tracts. No public data-and-code package currently reproduces that tract count or the diagnostics. The goal is to show where checks belong and why defaults require scrutiny.

The PySAL ecosystem offers esda for exploratory spatial statistics and libpysal.weights for neighborhood definitions. GeoPandas provides the geometry operations. The corrected workflow below states each assumption explicitly while using maintained implementations for CRS transformations, centroids, weights, and Moran inference. That matters because a hand-built distance matrix in longitude and latitude measures degrees, not ground distance, and a hand-copied Moran variance formula can combine incompatible null distributions without raising an error.

Use checks that report the intermediate inputs and outputs needed to detect these failures. The table below describes the instructional stack used in this article, not a publicly released project environment.

Tool Stack: GeoPandas, libpysal, and esda

Instructional tool stack; no public implementation package is currently available.
Component	Tool	Purpose
Geometries	Census TIGER/Line shapefiles (2023)	Tract boundary polygons
Spatial operations	`geopandas`	CRS alignment, point-in-polygon matching, projected centroids
Spatial weights	`libpysal.weights`	k-nearest-neighbor construction and row standardization
Spatial statistics	`esda`	Permutation-based global and local Moran diagnostics

Step 1: Load the Census Tract Shapefile

This teaching example uses the 2023 California TIGER/Line tract file. The statewide and county counts vary by vintage, so the run should record the source archive and verify the selected count. GeoPandas retains the file CRS, which gives us an explicit object to inspect before any join or distance calculation.

import geopandas as gpd
import numpy as np

tracts_all = gpd.read_file("tl_2023_06_tract.shp")
tracts = tracts_all.loc[tracts_all["COUNTYFP"] == "085"].copy()

if tracts.crs is None:
    raise ValueError("The tract layer has no declared CRS")
if tracts.geometry.isna().any() or tracts.geometry.is_empty.any():
    raise ValueError("The tract layer contains missing or empty geometries")

print(tracts.crs)
print(f"Loaded {len(tracts)} tracts for the target county")
# The article reports 108. Verify that count for the selected vintage.

# A point layer must declare its own CRS, then be transformed to the
# polygon CRS before a point-in-polygon join.
points = gpd.read_file("provider_locations.geojson")
if points.crs is None:
    raise ValueError("The point layer has no declared CRS")
points_for_join = points.to_crs(tracts.crs)
joined = gpd.sjoin(
    points_for_join,
    tracts[["GEOID", "geometry"]],
    how="left",
    predicate="within",
)

Step 2: Extract Tract Centroids from Polygon Geometries

With tract boundaries loaded, we need centroids for a distance-based neighborhood definition. Computing centroids or Euclidean distances in EPSG:4269 or EPSG:4326 would operate in angular degrees. Instead, estimate a local projected CRS, transform the polygons, and calculate centroids in metric coordinates. GeoPandas documents estimate_utm_crs() for this purpose. A project spanning several UTM zones should choose and justify a projection explicitly rather than rely on this convenience method.

Stop the workflow when any tract has a missing or empty geometry. The guard in Step 1 raises an error rather than skipping rows, so the analyst must repair or explicitly exclude each affected record before centroids and weights are constructed. The article reports that all 108 county tracts had valid geometries and that 3 tracts dropped from a separate statewide test. Neither run is publicly reproduced, so these counts illustrate why the guard matters rather than establish verified project results.

analysis_crs = tracts.geometry.estimate_utm_crs()
if analysis_crs is None:
    raise ValueError("Choose an appropriate projected analysis CRS explicitly")

tracts_metric = tracts.to_crs(analysis_crs).reset_index(drop=True)
centroids = tracts_metric.geometry.centroid
coords = np.column_stack([centroids.x, centroids.y])
geoids = tracts_metric["GEOID"].to_numpy()

print(analysis_crs)
print(f"Computed projected centroids for {len(coords)} tracts")

Step 3: Define Which Tracts Are Neighbors (k-NN Weights)

Spatial autocorrelation statistics require a weights matrix that defines which tracts are "neighbors." There are several ways to define this: contiguity (shared borders), distance bands, or k-nearest-neighbors.

Why k=8? This teaching example uses eight nearest neighbors as a starting specification, not a universal default. For the article-reported 108-tract case, 8 divided by 108 is roughly 7%, which produces a sparse weights matrix. An applied analysis should justify the choice and test whether conclusions change under alternative neighborhood definitions.

For contiguity weights, libpysal can construct shared-boundary graphs directly from the GeoDataFrame with Queen.from_dataframe(tracts_metric, use_index=True) or Rook.from_dataframe(tracts_metric, use_index=True). Compare those choices with k-nearest-neighbor and fixed-distance definitions when the substantive mechanism does not select one in advance.

The code below also row-standardizes the weights. This makes every non-isolated row sum to one, so each observation's spatial lag is expressed on a comparable scale. It does not by itself establish that the chosen neighbor definition is appropriate.

from libpysal.weights import KNN

n = len(coords)
k = 8
if n <= k:
    raise ValueError("k must be smaller than the number of tracts")

# coords are in the projected analysis CRS, so Euclidean distances
# are measured in that CRS's linear units rather than angular degrees.
w = KNN.from_array(coords, k=k)
w.transform = "R"  # row-standardize non-isolated rows

Step 4: Test Whether SNAP Participation Clusters Spatially

Is SNAP participation spatially clustered, or consistent with spatial randomness? Global Moran's I summarizes whether similar values tend to occur near one another under a specified weights matrix. Positive values generally indicate clustering and negative values generally indicate dispersion, but the attainable range and null distribution depend on the data and weights.

from esda import Moran, Moran_Rate

snap_rates = np.asarray(snap_rates, dtype=float)
if len(snap_rates) != n or not np.isfinite(snap_rates).all():
    raise ValueError("Rates must align one-to-one with the weights")

np.random.seed(20260721)
global_moran = Moran(snap_rates, w, permutations=9_999)
print(f"Moran's I: {global_moran.I:.4f}")
print(f"Permutation pseudo-p: {global_moran.p_sim:.4f}")

# When the outcome is an event rate with unequal denominators, retain
# the event counts and populations and assess the rate-specific statistic.
# rate_moran = Moran_Rate(events, population, w, permutations=9_999)

This permutation procedure holds the tract locations and weights fixed and randomly reassigns the observed values across locations. PySAL documents p_sim as a one-sided extreme-tail pseudo-p-value under that spatial-randomness reference distribution; the code reports that value without converting it to a two-sided probability. It is not a test of a causal process, and its interpretation depends on the pre-specified weights and permutation convention. For unstable rates with unequal exposure denominators, PySAL also provides esda.Moran_Rate; the choice between a raw rate and an adjusted rate statistic should be made from the measurement process, not after seeing significance.

The article reports a Moran's I above 0.3 for SNAP participation and uses Moran's I = 0.32 with p < 0.001 as its concrete example. No public inputs, weights matrix, output, or test log currently reproduces those values. If a verified run produced that combination under a valid inference procedure, it would be evidence of positive spatial autocorrelation under the specified weights. It would not by itself establish causality, and a model that assumes independent errors would need further spatial diagnostics.

The article uses Moran's I = 0.32 and z = 3.1 as an illustration, then reports that a development cross-check found an off-by-one error that shifted p-values by about 0.02. That debugging account is not backed by a public test log. The corrected teaching code therefore does not reproduce the handwritten variance calculation. PySAL separately documents analytical normality and randomization approximations and permutation inference; those are different reference distributions and must not be blended into one unnamed p-value.

Step 5: Identify Local Clusters via LISA

A statistically supported global Moran's I can indicate that clustering exists, but it does not locate the clusters. Local Indicators of Spatial Association (LISA) decompose the pattern into observation-level contributions. A common quadrant classification distinguishes High-High, Low-Low, High-Low, and Low-High patterns. The code below combines quadrants with unadjusted permutation screens; a final local-inference plan must also specify multiplicity control.

from esda import Moran_Local

local = Moran_Local(
    snap_rates,
    w,
    permutations=9_999,
    seed=20260721,
)
quadrant_name = {1: "HH", 2: "LH", 3: "LL", 4: "HL"}

# These are unadjusted permutation screens, not final discoveries.
cluster_labels = np.array(["NS"] * n, dtype=object)
for i, (quadrant, p_value) in enumerate(zip(local.q, local.p_sim)):
    if p_value < 0.05:
        cluster_labels[i] = quadrant_name[quadrant]

Moran_Local.p_sim is also a one-sided conditional-randomization screen. The seed belongs in the constructor because PySAL's numba implementation does not reliably honor an external NumPy seed. The article reports a handful of High-Low tracts in otherwise low-SNAP areas and interprets them as possible localized pockets of food insecurity. No public cluster table, map, local p-values, or multiple-testing procedure currently reproduces that account. Treat it as an article-reported example of a hypothesis to investigate, not a verified empirical finding.

Step 6: Validate Against Known Urban/Rural Patterns

Do the spatial clusters align with external knowledge about the county's geography? An unexpected pattern should prompt checks of inputs, coordinate systems, classifications, and model choices. It can also be a genuine result, so disagreement with prior expectations is a diagnostic signal rather than proof that something went wrong.

# Quartile classification for quick visual validation
quartiles = np.percentile(snap_rates, [25, 50, 75])
classifications = np.digitize(snap_rates, quartiles)

# Cross-tabulate with screened local labels
for q in range(4):
    mask = classifications == q
    hh_count = np.sum(cluster_labels[mask] == "HH")
    print(f"Quartile {q+1}: {hh_count} hot-spot tracts out of {mask.sum()}")

This cross-tab is an internal bookkeeping check, not external validation: HH and LL labels already depend on whether a tract is above or below the mean. Investigate impossible IDs, missing rows, and unexpected label counts, then compare mapped patterns with an independent source or known fixture. Agreement with prior expectations alone cannot validate the analysis.

The snippets identify the required validation points, but they do not reproduce the article's empirical results. A release-ready analysis would publish the input vintage, environment, full script, weights specification, tests, and saved outputs together.

What Can Go Wrong

Several failure modes here produce plausible-looking output rather than errors:

Incompatible coordinate reference systems can produce empty or incorrect joins. The example TIGER/Line file uses NAD83 (EPSG:4269), while many point datasets use WGS84 (EPSG:4326). Transform the point layer to the polygon CRS before joining, then use a justified projected CRS for centroids and Euclidean distance.
An unnamed Moran null produces an uninterpretable p-value. Normality approximations, randomization approximations, and permutations are distinct procedures. State which values are exchangeable, which features are fixed, how the weights were chosen, and whether the reported permutation probability is one- or two-sided.
Large shapefiles can require substantial memory. Filter by geography before retaining complex geometries when the analysis does not require the full source file, and measure memory use in the actual environment.

Given these pitfalls, most of which produce failures with no error rather than crashes, when does this workflow make sense?

When to Use This Approach

The workflow keeps the design decisions explicit while delegating numerical spatial statistics to maintained libraries. Whether it fits depends on the project.

Good fit:

Moderate-scale teaching or exploratory analyses where spatial dependence must be diagnosed before modeling
Projects where CRS, neighborhood, and inference choices can be logged and reviewed explicitly
Exploratory work where understanding the mechanics of spatial statistics, rather than the output alone, is the goal

Less suitable:

Large-scale production workflows that require a versioned environment, batch tests, memory benchmarks, and operational monitoring beyond these standalone snippets
Analyses whose neighborhood mechanism depends on a routed network, travel time, or another relation that polygon contiguity and centroid distance do not represent
Projects that require a complete released implementation and frozen inputs rather than an article-only teaching workflow

Limitations

Using maintained spatial operations reduces implementation risk, but the analyst still owns the geographic and inferential assumptions.

Centroid neighborhood definition. Even projected geometric centroids can fall outside elongated or non-convex polygons. Contiguity or another substantively justified distance definition may be preferable.
Fixed k for all tracts. Eight neighbors can span a much longer distance in sparse rural areas than in dense urban areas. That changing geographic reach may or may not match the mechanism; compare it with contiguity and fixed-distance alternatives.
Permutation exchangeability. The permutation test assumes the observed values are exchangeable across the fixed locations under the spatial-randomness null. Spatially varying variance or rate precision can make that reference distribution a poor description of the data-generating process.
Boundary truncation. A county-only k-NN graph still assigns eight neighbors, but it excludes potentially closer tracts across the county line. Build weights on the substantively relevant surrounding geography before subsetting results when cross-boundary interaction matters.
Single time period. This workflow analyzes a cross-sectional snapshot. Temporal dynamics in SNAP participation or store openings/closures require panel-data extensions.

Code Availability

No public reproduction package is currently available. As of July 21, 2026, the linked foodsecurity_mobility repository contains only a one-line README. It does not include the previously claimed scripts/phase5_spatial_validation.py, data preprocessing, visualization scripts, environment specification, or saved outputs. The snippets on this page are instructional examples and should not be described as the full project implementation.

References and Notes

Source note. TIGER/Line shapefiles are produced by the U.S. Census Bureau and available at census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html. The code example names the 2023 California tract vintage, but the exact source archive and checksum used for the article have not been published.
Inference note. PySAL's official esda.Moran documentation distinguishes normality, randomization, and permutation outputs. Its esda.Moran_Local documentation also requires the reproducibility seed inside the constructor. This page uses the documented one-sided permutation output and states its spatial-randomness null.
Count note. This article reports 108 tracts for the study county. That count is not currently reproduced by a public source manifest or run output and may differ across TIGER/Line vintages.