Spatial Join: which zone is this point in, really?

Project at a glance

• Track: Geospatial
• Difficulty: Beginner · 3/10
• Time: ~1.0 h for Part A · ~2.5 h more for the optional Part B
• You’ll use: geopandas, shapely, scipy, pandas, numpy, pytest
• 💻 Local setup needed (no in-browser ▶ Run yet) — a one-time, free Python install, then you just run files. Start Here walks you through it.
• New to this? This is a beginner project — you can’t break anything. Never coded? You can still run and read every step — the writing is only the optional Part B. This one needs a one-time free setup on your computer; Start Here walks you through it.

The data in this project is synthetic — geometry layers generated locally by data/make_data.py, with no external basemap. That is deliberate: grading a spatial join needs the known-correct answer (which zone contains each point, which service areas cover it), which real map data doesn’t come labelled with — and it keeps the project download-free. No signup, no GPU.

Before you start

New to geospatial joins? Start here.

A spatial join links two layers of geometry by their spatial relationship instead of a shared key. The most common one: given points (deliveries, customers, incidents) and polygons (zones, districts, service areas), tag each point with the polygon that contains it — point-in-polygon.

The tempting shortcut is to assign each point to the polygon whose centre is nearest. It’s fast and it’s wrong: a zone’s centre doesn’t define its area. A point deep inside a big, wide zone can be closer to the centre of a small zone next door — so “nearest centre” mislabels it. Only true containment (does this polygon’s outline enclose this point?) is correct.

A tiny example — a wide zone Z1 and a small zone Z3 sitting above it:

Z3 centre •           a point here (just inside Z1, near the top edge)
Z1 ............X....   is closer to Z3's centre than to Z1's centre,
Z1 centre •            so "nearest centre" says Z3 — but it's inside Z1.

What you’ll learn: do a real point-in-polygon spatial join with GeoPandas, see how far the nearest-centre shortcut drifts, then handle the harder case where areas overlap and a point can belong to several at once. Section 3 lists exactly what you’ll be able to do.

New to Python too? See Start Here — it starts from zero.

1. The Brief

You’ve got a pile of delivery coordinates and a set of operating zones, and every downstream report — deliveries per zone, cost per district, which depot serves what — depends on tagging each delivery with the right zone. You’ll do it the correct way (containment), show how badly the nearest-centre shortcut distorts the per-zone counts, then extend the join to overlapping service areas where a point can be covered by more than one facility.

Spatial joins are the backbone of geospatial analytics — every choropleth, catchment analysis, and “points within region” query is one. Getting the predicate right (contains vs within-distance vs intersects) is the whole skill, and it’s identical whether the geometry is delivery zones, sales territories, or flood plains.

2. Difficulty & time

Difficulty 3 / 10. One core operation (the spatial join) with one genuine judgment — use true geometry, not a proximity proxy — and a fully specified path in Part A. Calibrated against the ledger (G10): a peer of silent-row-loss, idempotent-loads, and duplicate-invoices (all 3) — single-mechanic builds where the lesson is a specific data-handling judgment (here: containment ≠ nearest-centre), not a pipeline or a model. It sits below maverick-spend (4), which layers a conditional rule and a second detection type. It stays above the 1–2 band because the naive shortcut looks right and silently biases every per-zone number, and Part B needs a real change of predicate (containment → within-radius, one → many).

Time is separate from difficulty: Part A is about an hour; Part B is 2–3 hours. The scripts run in under a second.

3. What you’ll be able to do after

• Do a point-in-polygon spatial join with GeoPandas (sjoin, predicate="within").
• Explain why nearest-centre is wrong — and quantify how much it distorts per-zone counts.
• Change the join predicate — from containment to within-a-distance, and handle a point matching several areas.
• Report coverage — for overlapping service areas, the full set of areas covering each point.

The finished result

By the end of Part A you’ll have the true per-zone counts from a real point-in-polygon join, and the proof that the nearest-centre shortcut quietly biases them:

containment (correct):   Z1 594   Z2 354   Z4 278   Z5 148   Z3 126
nearest-centre shortcut agrees on only 88.4% of points -- 174 misassigned,
   biased toward the small central zones and against the big/irregular ones

The shortcut looks fine and runs fast — but ~12% of points land in the wrong zone, and every per-zone number built on it inherits that bias. Containment is the whole point.

4. Prereqs & time box

You can write a function and use a pandas DataFrame. No geospatial or GIS background needed — Part A introduces GeoPandas and Shapely. Coordinates are plain planar (x, y in km), so no map-projection knowledge is required.

Where to write and run code

Run this project on your computer. It uses GeoPandas and Shapely (the geospatial libraries); the in-browser ▶ Run option isn’t confirmed for GeoPandas yet, so this is a local-only project for now.

The one-time setup — installing Python, downloading the code from the course site (no GitHub, no git), opening a terminal, creating a virtual environment, installing the requirements, and running Part A step by step — is covered once in Start Here. Come back here once python -c "import geopandas" runs without error.

Time box

• Part A — Guided Build: 60–90 minutes. Hard cap.
• Part B — Your Turn: 2–4 hours. Hard cap.

5. Setup & data

The data is synthetic and generated locally — see the banner above and SOURCES.md for full provenance. One command creates both parts’ layers:

python data/make_data.py

That writes four CSVs into data/ (git-ignored; regenerate anytime, byte-for-byte identical because the seed is fixed):

• Part A — zones_a.csv (5 zones: zone_id, geometry as WKT polygons) and points_a.csv (1,500 points: point_id, x, y).
• Part B — facilities_b.csv (25 facilities: facility_id, x, y, radius) and demand_b.csv (2,000 points: point_id, x, y). Service circles overlap.

Coordinates are a projected planar grid (km) — treat them as ordinary x/y geometry.

6. Part A — Guided Build

You’ll join the 1,500 points to their zones. Run the whole thing, or (better) one # %% cell at a time — see Start Here (this project runs on your computer, not in the browser):

python part_a/join_points.py

The full, commented script is in part_a/join_points.py. The spine:

Step 1 — Load the layers. Build GeoDataFrames from the WKT polygons and the x/y points.

Checkpoint: 5 zones of very different areas (Z1 the wide strip is biggest, Z3 the small square smallest), 1,500 points.

Step 2 — The shortcut: nearest zone centre. Assign each point to the closest centroid.

Checkpoint: it runs and looks plausible — but a centre doesn’t define an area.

Step 3 — The spatial join. gpd.sjoin(..., predicate="within").

Checkpoint: true per-zone counts — Z1 ~594, Z2 ~354, Z4 ~278, Z5 ~148, Z3 ~126.

Step 4 — How wrong was the shortcut?

Checkpoint: nearest-centre agrees with containment on only ~0.884 of points (~174 misassigned), biased toward small central zones and against the big/irregular ones. Writes part_a/assigned_a.csv.

Why this and not that

• Why containment, not nearest centre. A polygon’s centroid is a single point; its area is the whole outline. Proximity to the centre and being inside the outline are different questions, and they diverge most for large or irregular zones — exactly where the errors cluster.
• Why the error matters. ~12% misassignment isn’t noise: it systematically inflates small zones and shrinks big ones. Every per-zone metric built on it — counts, rates, choropleths, cost allocation — inherits that bias.
• Why GeoPandas sjoin. It runs the containment test with a spatial index, so it’s both correct and fast, and it generalises to other predicates (intersects, contains, dwithin) — including the one Part B needs.

Functions you’ll meet (and where to read more)

Function	What it does here	Docs
gpd.GeoDataFrame(..., geometry=...)	Attach geometry to a table	GeoDataFrame
shapely.wkt.loads(s)	Parse a WKT polygon string	shapely.wkt
gpd.points_from_xy(x, y)	Build point geometry from coordinates	points_from_xy
gpd.sjoin(a, b, predicate="within")	Point-in-polygon spatial join	sjoin
geom.centroid	A polygon’s centre (for the naive baseline)	centroid

If something looks off

• Some Part B points have an empty facilities set — is that a bug? No. ~96 demand points sit outside every service circle, so “covered by no facility” is the correct answer for them (an empty set, not an error). If every point has at least one, you’re snapping to the nearest instead of testing the radius.
• ModuleNotFoundError: No module named 'geopandas'. The install step was skipped — run pip install -r requirements.txt inside your activated virtual environment (see Start Here).
• A checkpoint number doesn’t match. The data is seeded, so the counts are exact. Re-run python data/make_data.py (you may have skipped or edited it), then re-run Part A.

7. Part B — Your Turn

Part B is the optional “your turn” half — a guided fill-in-the-blank, not a blank page. You get a working starter (part_b/starter.py) that already handles the easy version; your job is to finish the # TODO it marks out. More open than Part A’s step-by-step, but you are never staring at an empty file. If you finished Part A and want to prove the skill, this is where it happens. New to Python? It’s completely fine to stop after Part A and come back later.

Same core skill, a different predicate. Each facility has a circular service area (radius), and the circles overlap — so a demand point can be covered by zero, one, or several facilities. For every point in demand_b.csv, report the full set of facilities that cover it. The starter assigns each point to its nearest facility only and scores ~0.18 — it can’t represent multi-coverage, and here most points sit in two or more overlapping areas.

The transfer: Part A’s join was “which one polygon contains me” (each point in exactly one zone); Part B is “which of the overlapping circles am I within” (each point in many). The predicate moves from containment-of-one to within-radius-of-many.

Acceptance criteria (checkable)

Your solution writes into part_b/out/:

• coverage.csv — columns point_id, n_covering, facilities (the covering facility ids joined by |, empty if none).
• report.json — with keys demand_points, achieved_metric (you can leave achieved_metric null; the grader scores you).

Then:

python part_b/starter.py
pytest part_b/test_solution.py -q

If pytest reports everything as skipped, it just means the output files don’t exist yet — run python part_b/starter.py first (once your solution fills in the TODOs), then re-run.

You’re done when the tests pass: ≥ 90% of points have the exact correct covering set (the reference reaches 1.00). Nearest-facility-only scores ~0.18 — reporting all facilities within the radius is what clears the floor.

Constraints — what must be true of your solution, not how to get there

• Report every facility whose service area covers the point, not just the closest.
• Use the radius from facilities_b.csv (a point is covered iff its distance ≤ radius).
• A point covered by no facility gets an empty set (that’s a valid answer, not an error).
• Emit one row per demand point.

Hints

Hint 1 — the right predicate

Coverage is a within-distance test, not nearest. For each demand point you need all facilities with distance ≤ radius — a range query, not an argmin.

Hint 2 — do it fast

scipy.spatial.cKDTree(facility_xy).query_ball_point(point_xy, r=radius) returns, for each point, the indices of all facilities within the radius. (GeoPandas sjoin with predicate="dwithin" or buffered polygons works too.)

Hint 3 — build the set (still not the code)

Map those indices back to facility_id, sort them, and join with |. Points with an empty result are genuinely uncovered — leave their facilities blank.

8. Self-check

You don’t need the answer key. You’re done when all of these hold:

• pytest reports set accuracy ≥ 0.90 (aim for 1.00 — the distance test is exact).
• Your n_covering has plenty of 2s and 3+s, not all 0/1 — if it’s capped at 1, you’re still doing nearest-only.
• Re-running from python data/make_data.py reproduces the same numbers — the data is seeded.

And the tell-tale: a handful of points (~96) are covered by no facility. If every point has at least one, you’re snapping to the nearest instead of testing the radius.

9. Stretch

• Buffered-polygon join. Rebuild Part B as a GeoPandas sjoin using each facility’s circle as a buffered polygon (geom.buffer(radius)), and confirm it matches the KDTree result.
• Nearest within a limit. Report, per point, the single closest facility but only if within the radius — and compare that assignment to the full-coverage view.
• (Genuinely hard) Real coordinates. Switch to lat/lon and use great-circle distance (or project to a metric CRS with to_crs) so “within 22 km” is correct on the globe, not on a flat grid.

10. Ship it

Put this in a portfolio repo and it counts. In that repo’s README:

• State the finding in one sentence an ops lead gets: “nearest-centre mis-tagged 12% of deliveries; a proper point-in-polygon join fixed the per-zone counts.”
• Show the naive-vs-containment count table and the overlapping-coverage result side by side.
• Include the sjoin and one paragraph on choosing the predicate (contains vs within-distance). That paragraph is what signals you understand spatial joins, not just GeoPandas syntax.

Keep make_data.py and the tests in the repo so anyone can reproduce your numbers from zero.

11. Sources

Rendered from SOURCES.md.

Field	Value
Dataset	Synthetic geospatial layers on a projected 100x100 km grid (Part A 5 zones + 1,500 points; Part B 25 facilities + 2,000 points)
Type	Synthetic
Generated by	data/make_data.py (this repo), seeds 20260740 (A) and 20260741 (B)
Source URL	Not applicable — generated, not fetched
License	Not applicable — synthetic, no third-party rights
Access / generation date	2026-07-18
Redistribution	Permitted — synthetic; output regenerated from the seeded script, never committed
robots.txt checked	Not applicable — no scraping
ToS URL	Not applicable — no external source

Synthetic data teaches the mechanic; it does not prove a finding about the real world. The zone shapes and service areas were designed on purpose — the join technique is the point, not a claim about any real geography.