What Is a Spatial Join? Types and SQL Examples

Q: How do you do a spatial join in Python?

With GeoPandas, call geopandas.sjoin with a predicate such as intersects or within and a how of inner, left, or right. GeoPandas runs on one machine. For tables too large for memory, run the same join in PySpark with WherobotsDB, using ST_Intersects or another predicate in the join condition.

A spatial join combines two tables by where their geometries sit. It pairs rows whose shapes meet a spatial condition, with no shared ID column needed: a point that falls inside a polygon, two parcels that share a boundary, a facility within 2 kilometers of a river. Most data about the physical world lacks a common key across sources, which makes the spatial join the operation that connects building footprints to flood zones, delivery stops to sales territories, and satellite pixels to farm fields.

Key takeaways

A spatial join matches rows on a spatial predicate, such as ST_Intersects, ST_Contains, or ST_DWithin, between two geometry or geography columns.
Both sides need geometry in the same coordinate reference system (CRS), and argument order matters for every predicate except ST_Intersects, ST_Touches, and ST_DWithin.
Spatial joins come as inner, left, and full outer joins. A left join keeps every left row, including those with no match.
The four join types are predicate joins, distance joins, k-nearest neighbor joins, and raster-vector joins.
The cost of a naive spatial join grows with the product of both table sizes, so engines prune candidate pairs before testing exact geometry.

How a spatial join works

Start from a query. This one counts the points of interest inside each Colorado city and town, using open data from the Overture Maps Foundation:

WITH localities AS (
  SELECT id, names.primary AS city, geometry
  FROM wherobots_open_data.overture_maps_foundation.divisions_division_area
  WHERE subtype = 'locality' AND region = 'US-CO'
)
SELECT l.id, l.city, COUNT(p.id) AS places
FROM localities AS l
JOIN wherobots_open_data.overture_maps_foundation.places_place AS p
  ON ST_Intersects(l.geometry, p.geometry)
GROUP BY l.id, l.city

Three things make that JOIN spatial.

A geometry or geography column on each side. A geometry column stores points, lines, or polygons as coordinates on a plane. A geography column stores them on the curved surface of the Earth, so distances come back in meters. A latitude and longitude held in two number columns becomes a geometry only after a constructor such as ST_Point builds it.

The same coordinate reference system. A CRS defines how coordinates map to places on the Earth. EPSG:4326 stores longitude and latitude in degrees. EPSG:2263 stores feet on the New York Long Island State Plane. Comparing geometries from two different systems returns rows, and those rows are wrong. ST_Transform converts one side to match the other.

A spatial predicate. The predicate in the ON clause decides which pairs match. Each left row can match zero, one, or many right rows, so a point on a shared city line can appear twice in the output.

Spatial join predicates

Predicates fall into two groups. Topological predicates test how two shapes relate in space. Distance predicates test how far apart they are.

Predicate	Returns true when	Typical use
`ST_Intersects(A, B)`	A and B have at least one point in common, boundary included	Point in polygon, overlapping polygons, lines that meet
`ST_Contains(A, B)`	Every point of B lies in A, and their interiors share a point	Assigning assets to a territory
`ST_Within(B, A)`	B lies inside A. It is the converse of `ST_Contains`	Filtering records to a boundary
`ST_Covers(A, B)`	No point of B lies outside A, so points on A's boundary count	Containment where an edge point should match
`ST_Touches(A, B)`	A and B share boundary points and no interior points	Finding adjacent parcels
`ST_Crosses(A, B)`	A line passes through a polygon or another line	A pipeline crossing a river
`ST_DWithin(A, B, d)`	A and B lie within distance `d`, where `d` can be a column that varies per row	Proximity with a radius per record

ST_Intersects runs in more production joins than any other predicate. PostGIS defines it as two geometries having any point in common, which makes it the broadest topological test: it matches a point inside a polygon, a point on the polygon's edge, a line that crosses a boundary, and two polygons that overlap. It is also symmetric, so ST_Intersects(A, B) and ST_Intersects(B, A) return the same rows, and no argument order can be wrong.

ST_Covers exists because of a quirk in ST_Contains. Under the PostGIS definition, a polygon does not contain a point that sits exactly on its boundary, because the interiors share no point. ST_Covers only asks that no point of B lie outside A, so the boundary point matches.

Inner, left, and full outer spatial joins

The join type decides what happens to rows that find no match.

Inner join keeps only matched pairs. In the query above, a city with no points of interest drops out of the result.
Left join keeps every row of the left table. Unmatched rows carry nulls in the right table's columns, and COUNT(p.id) returns 0 for them.
Full outer join keeps unmatched rows from both tables: cities with no places, and places outside every city.

On Overture data, the inner join above returns Colorado localities that hold at least one place. Switching JOIN to LEFT JOIN returns all 304 localities, including the 11 with no places at all. Those 11 rows are often the finding: towns with no mapped businesses, territories with no customers, zones with no sensors. The filter on localities sits inside the WITH clause, before the join, so it never removes the unmatched rows.

A full outer join needs both sides filtered before the join, and it needs an ID from each side in the output, because an unmatched place has no city to group by:

-- Boulder area
WITH localities AS (
  SELECT id, geometry
  FROM wherobots_open_data.overture_maps_foundation.divisions_division_area
  WHERE subtype = 'locality' AND region = 'US-CO'
    AND bbox.xmin BETWEEN -105.35 AND -105.0 AND bbox.ymin BETWEEN 39.9 AND 40.1
),
places AS (
  SELECT id, geometry
  FROM wherobots_open_data.overture_maps_foundation.places_place
  WHERE bbox.xmin BETWEEN -105.3 AND -105.2 AND bbox.ymin BETWEEN 39.98 AND 40.06
)
SELECT l.id AS locality_id, p.id AS place_id
FROM localities l
FULL OUTER JOIN places p ON ST_Intersects(l.geometry, p.geometry)

Rows with a null place_id are localities with no places. Rows with a null locality_id are places outside every locality. Around Boulder that comes to 11,914 matched pairs, 6 localities with no places, and 321 places outside every locality.

A left join followed by WHERE p.id IS NULL is an anti join. It returns the rows with no spatial match, such as buildings outside every flood zone polygon, and it is the correct way to ask "what lies outside." ST_Disjoint in a join condition matches almost every pair and produces an enormous result.

In WherobotsDB, predicate joins, distance joins, and raster-vector joins all accept LEFT JOIN and FULL OUTER JOIN. The ST_KNN join is an inner join on the query table: a query row with no neighbor is left out.

The four types of spatial join

Predicate join

A predicate join, also called a range join, pairs rows whose geometries satisfy a topological predicate. Point in polygon is the standard case.

Parameters: the two geometry columns, in the order the predicate expects.
Example: ON ST_Contains(zones.geometry, assets.geometry) assigns each asset to the zone that contains it.

Distance join

A distance join pairs rows whose geometries lie within a distance of each other.

Parameters: two geometry columns and a distance. The distance can be a constant or a column, so each row can carry its own radius.
Units: for geometry, the distance follows the CRS. On EPSG:4326 data, ST_DWithin(a, b, 2) means 2 degrees, a span of hundreds of kilometers.
Meters: transform both sides to a projected CRS in meters, or pass useSpheroid = true, which measures between the centroids of the two geometries. On the geography type, ST_DWithin measures the minimum great-circle distance between the shapes, in meters.
Example: ON ST_DWithin(stations.geom_m, homes.geom_m, stations.walk_radius_m) matches homes within each station's own walking radius.

Nearest-neighbor join

A nearest-neighbor join answers a different question: for each row on one side, which k rows on the other side are closest? It needs no radius and no containment.

ST_KNN(R, S, k, use_sphere, search_radius) takes five arguments:

R: the query geometry, one row per question, such as each customer.
S: the object geometry to search, such as every store.
k: the number of neighbors to return for each query row.
use_sphere: the distance model. false uses planar distance in the data's coordinate units and measures the true minimum distance between shapes, boundary to boundary. true uses great-circle distance in meters and first reduces each line or polygon to its centroid, which shifts results for large or elongated shapes.
search_radius (optional): the farthest distance to search, in the units of the chosen model.

Planar distance on EPSG:4326 data is measured in degrees, and a degree of longitude shrinks toward the poles, so the ranking can differ from ground distance. For nearest features in meters, transform both sides to a projected CRS for the area first. ST_AKNN takes the same arguments and returns approximate neighbors faster on large tables.

Raster-vector join

A raster stores a value in every cell of a grid: elevation, reflectance, land cover class. A raster-vector join pairs raster tiles with the vector geometries they overlap, using RS_Intersects, RS_Contains, or RS_Within, and then summarizes the pixels inside each geometry with zonal statistics such as RS_ZonalStats.

Parameters: a raster column, a geometry column, and for the summary a band number and a statistic such as mean, sum, or count.
Example: mean elevation per town from a digital elevation model, or the share of each field covered by healthy vegetation from a Sentinel-2 NDVI raster.

Spatial join vs. attribute join

	Attribute join	Spatial join
Join condition	Equal values in a shared column	A spatial predicate between two geometries
Requirement	A common key in both tables	Geometry in both tables, in one CRS
Example	Policies to customers on `customer_id`	Policies to flood zones by location
Result size	Bounded by key matches	One row per matching pair, which can exceed both inputs

Attribute joins and spatial joins often run in the same pipeline. A spatial join assigns each building a parcel ID, and an attribute join on that ID then brings in the assessor's records.

How to run a spatial join in Wherobots

The examples below run in a Wherobots notebook or SQL session against the open data catalog. Python cells start a Sedona session first.

from sedona.spark import SedonaContext

config = SedonaContext.builder().getOrCreate()
sedona = SedonaContext.create(config)

Point-in-polygon join with a left join

Count Overture places in every Colorado locality, keeping localities with none.

places_per_city = sedona.sql("""
WITH localities AS (
  SELECT id, names.primary AS city, geometry
  FROM wherobots_open_data.overture_maps_foundation.divisions_division_area
  WHERE subtype = 'locality' AND region = 'US-CO'
),
places AS (
  -- Prefilter places to Colorado's bounding box before the join
  SELECT id, geometry
  FROM wherobots_open_data.overture_maps_foundation.places_place
  WHERE bbox.xmin BETWEEN -109.06 AND -102.04 AND bbox.ymin BETWEEN 36.99 AND 41.0
)
SELECT l.city, COUNT(p.id) AS place_count
FROM localities l
LEFT JOIN places p ON ST_Intersects(l.geometry, p.geometry)
GROUP BY l.id, l.city
ORDER BY place_count DESC
""")
places_per_city.show(10)

Nearest-neighbor join in meters

Find the nearest hospital to every school in a box around central Denver, which takes in parts of neighboring cities. Both sides move to UTM zone 13N (EPSG:26913), so planar distances come back in meters.

WITH schools AS (
  SELECT id, names.primary AS school,
         ST_Transform(geometry, 'EPSG:4326', 'EPSG:26913') AS geom_m
  FROM wherobots_open_data.overture_maps_foundation.places_place
  WHERE taxonomy.primary IN ('elementary_school', 'middle_school', 'high_school')
    AND bbox.xmin BETWEEN -105.11 AND -104.6 AND bbox.ymin BETWEEN 39.61 AND 39.91
),
hospitals AS (
  -- Search a wider area so schools near the edge still find their true nearest hospital
  SELECT id, names.primary AS hospital,
         ST_Transform(geometry, 'EPSG:4326', 'EPSG:26913') AS geom_m
  FROM wherobots_open_data.overture_maps_foundation.places_place
  WHERE taxonomy.primary = 'hospital'
    AND bbox.xmin BETWEEN -105.3 AND -104.4 AND bbox.ymin BETWEEN 39.4 AND 40.1
)
SELECT s.school, h.hospital, ROUND(ST_Distance(s.geom_m, h.geom_m)) AS distance_m
FROM schools s
JOIN hospitals h ON ST_KNN(s.geom_m, h.geom_m, 1, false)
ORDER BY distance_m DESC

On Overture data the query pairs 596 schools in that box with a hospital. To report on one city, filter the schools with ST_Intersects against its boundary first. The average school sits 1,589 meters from its nearest hospital, and the farthest sits 9,915 meters away.

Raster-vector join with zonal statistics

Compute the mean elevation of every Colorado city and town from the Copernicus 30 m DEM. A town can span several raster tiles, so the query sums pixel values and pixel counts per tile and divides once per town, which weights each tile by its share of the town.

WITH towns AS (
  SELECT id, names.primary AS town, geometry
  FROM wherobots_open_data.overture_maps_foundation.divisions_division_area
  WHERE subtype = 'locality' AND region = 'US-CO'
),
dem AS (
  -- Keep only elevation tiles that touch Colorado
  SELECT rast
  FROM wherobots_open_data.copernicus_dem.glo_30m
  WHERE ST_Intersects(footprint, ST_GeomFromText(
    'POLYGON((-109.1 36.9, -102 36.9, -102 41.1, -109.1 41.1, -109.1 36.9))'))
),
per_tile AS (
  SELECT t.id, t.town,
         RS_ZonalStats(d.rast, t.geometry, 1, 'sum') AS elev_sum,
         RS_ZonalStats(d.rast, t.geometry, 1, 'count') AS cells
  FROM towns t
  JOIN dem d ON RS_Intersects(d.rast, t.geometry)
)
SELECT town, ROUND(SUM(elev_sum) / SUM(cells)) AS mean_elevation_m
FROM per_tile
GROUP BY id, town
ORDER BY mean_elevation_m DESC

Alma comes out at a mean of 3,185 meters and Leadville at 3,088 meters. The same pattern summarizes canopy height, land cover, or any other raster inside parcels, fields, or flood zones.

Why spatial joins get expensive at scale

An attribute join on a key can hash one side and probe it with the other, one pass over each table. A spatial predicate has no hash. Tested naively, every left row is compared with every right row: 10 million buildings against 50,000 flood zone polygons is 500 billion geometry tests, and each test walks coordinate arrays that can hold thousands of vertices for a detailed coastline or county boundary.

Every spatial engine attacks this the same way at a high level. It compares cheap bounding boxes first, discards pairs whose boxes cannot touch, and runs the exact predicate only on the candidates that remain. The engines differ in how well that pruning scales across many machines, how they handle skew, where a few dense cities hold most of the points, and how much tuning they demand from the user. A join that works on a sample in a single-node database such as PostGIS can stall on a full country.

Spatial joins in WherobotsDB

WherobotsDB is the spatial engine in Wherobots, built by the original creators of Apache Sedona and compatible with Sedona's APIs. It plans spatial joins automatically: the docs describe spatial indexing as managed for you, without manual tuning. The same SQL runs on a sample and on a national table.

Vector and raster in one query. Predicate, distance, KNN, and raster-vector joins run in the same SQL, so a pipeline can join buildings to parcels, parcels to flood zones, and flood zones to an elevation raster without moving data between tools.
Join types beyond inner. Predicate, distance, and raster-vector joins take left and full outer joins, and ST_KNN and ST_AKNN add exact and approximate nearest-neighbor joins.
Measured speed. On SpatialBench at scale factor 1000, the current WherobotsDB runs spatial queries and joins up to 2.5 times faster than the previous generation, with a mean speedup of 1.9 times.

Spatial join use cases

Insurance and catastrophe risk. A join between policy locations and flood risk polygons returns exposed policies and insured value per zone, the first step of a catastrophe model.
Property and land records. Building footprints join to parcel data to give each building its assessor attributes, and parcels join to hazard layers to rank exposure.
Site selection. Candidate sites join to transmission lines, substations, and flood zones by distance and containment to screen data center sites.
Agriculture. Field boundaries join to Sentinel-2 imagery to give each field a vegetation index time series across a season.
Logistics. A KNN join assigns each delivery stop to its nearest depot, and a distance join flags stops outside every depot's service radius.
Environmental compliance. Operating sites join to protected areas and waterways to flag every overlap in one query.

Common spatial join mistakes

Mismatched coordinate reference systems. Two geometries in different CRSs still join and return rows, and the rows are wrong. Transform one side with ST_Transform first.
Degrees where meters were intended. ST_Distance and ST_DWithin on EPSG:4326 data measure degrees. Transform to a projected CRS in meters, or use the geography type.
Centroids where shapes were intended. Spheroid distance and ST_KNN with use_sphere = true measure between centroids for lines and polygons. A large polygon's centroid can sit kilometers from its nearest edge.
Reversed arguments. ST_Contains(zone, asset) asks whether the zone holds the asset. ST_Contains(asset, zone) asks the reverse and returns almost nothing for points.
Boundary sets that skip places. In Overture, Denver is a consolidated city and county, so it appears as subtype = 'county', and a query for subtype = 'locality' leaves it out. Check that the polygon layer covers every place the question names.
Inner joins that hide gaps. An inner join drops every row without a match. Use a left join when the unmatched rows matter.
Invalid geometry. Self-intersecting polygons return wrong answers or errors in the middle of a join. ST_IsValid finds them and ST_MakeValid repairs them.

Raster vs. vector data: the two data models a spatial join connects
Digital elevation model: a raster often summarized per polygon with zonal statistics
Parcel data: land records joined to buildings and hazards by location
Flood risk data: hazard zones joined to assets in insurance and planning
Geospatial AI: models whose outputs become geometries ready to join

Frequently asked questions

A spatial join combines two tables by a spatial relationship between their geometries, such as one shape containing, touching, or lying within a set distance of another. Each output row pairs a record from the left table with a record from the right table that meets the condition. A spatial join needs no shared ID column, because location does the matching.

A store point that falls inside a city boundary, a pipeline that crosses a river, two parcels that share a boundary, and a school within 500 meters of a highway are all spatial relationships. In SQL each one is a predicate: ST_Contains, ST_Crosses, ST_Touches, and ST_DWithin.

Use a spatial join when two datasets describe the same places but share no key. Common cases are assigning points to the polygon they fall in, counting features per zone, finding the nearest facility to each customer, and summarizing raster values such as elevation or land cover inside each polygon.

In ArcGIS Pro, run the Spatial Join tool in the Analysis toolbox. Set the target features and the join features, choose a match option such as Intersect or Within a distance, and choose Join one to one or Join one to many. The tool writes the target features and the joined attributes to a new feature class.

An attribute join matches rows on equal values in a shared column, such as a parcel ID or a ZIP code. A spatial join matches rows on a spatial predicate between two geometry columns. Use an attribute join when a reliable key exists, and a spatial join when the only link between the datasets is location.

Yes. An inner spatial join keeps only pairs that match. A left spatial join keeps every row of the left table and fills the right columns with nulls where nothing matches, and a full outer join keeps unmatched rows from both sides. In WherobotsDB, predicate joins, distance joins, and raster-vector joins take LEFT JOIN and FULL OUTER JOIN in SQL. The ST_KNN nearest-neighbor join is an inner join.

With GeoPandas, call geopandas.sjoin with a predicate such as intersects or within and a how of inner, left, or right. GeoPandas runs on one machine. For tables too large for memory, run the same join in PySpark with WherobotsDB, using ST_Intersects or another predicate in the join condition.

A KNN join returns, for each row in one table, the k closest rows in another table. In WherobotsDB, ST_KNN(queries, objects, k, use_sphere) takes the query geometry, the object geometry, the number of neighbors, and a distance model, with an optional search radius. ST_AKNN returns approximate neighbors faster on large tables.

What Is a Spatial Join? Predicates, Types, and SQL