Geospatial Data Engineering: PostGIS, H3, GeoParquet & Apache Sedona

Geospatial data engineering is what the discipline turns into the moment a latitude and longitude land in your warehouse and you discover none of your usual tools work — because a location is not a number you can GROUP BY, a "customers within five kilometres" query cannot use a B-tree index, and a join of a billion GPS pings against a few million delivery zones is not a join your Spark cluster knows how to run. A pair of coordinates looks like two floats, but the questions people ask of it — does this point fall inside that polygon, how many rides started within this neighbourhood, which store is nearest — are geometry problems, and geometry has its own type system, its own idea of distance, and its own family of indexes.

This guide is the senior-data-engineering walkthrough for building the location stack properly — framed the way interviewers actually probe it: why spatial indexing is a bounding-box problem and not a sorted-key one, how PostGIS turns Postgres into a full spatial engine with geometry/geography types, GiST indexes, and the ST_ function family that make a spatial join cheap, how H3 replaces polygon joins with a hexagonal spatial index so aggregation becomes an ordinary group-by, how GeoParquet stores geometry columnar so GIS tools and query engines can all read the same file, and how Apache Sedona partitions and distributes spatial joins across a Spark cluster when the data no longer fits one machine. Each section pairs a teaching block with a Solution-Tail interview answer — code, a step-by-step trace, an output table, then a concept-by-concept breakdown of why it works.

When you want hands-on reps immediately after reading, drill the spatial indexing practice library →, rehearse point-in-polygon work on the spatial join practice library →, and sharpen the architecture axis with the system design practice library →.

---

On this page

• Why geospatial data engineering is its own discipline
• PostGIS — spatial SQL over Postgres
• H3 — hexagonal hierarchical spatial indexing
• GeoParquet — columnar spatial storage
• Apache Sedona — distributed spatial joins on Spark
• Cheat sheet — geospatial data engineering
• Frequently asked questions
• Practice on PipeCode

---

1. Why geospatial data engineering is its own discipline

Geometry, coordinate systems, and spatial indexing — the three reasons your ordinary tools fail

The one-sentence invariant: geospatial data engineering is a distinct discipline because a location value is a geometry — a point, line, or polygon in a specific coordinate reference system — not a scalar, so the three hard problems are storing that geometry with its CRS attached, asking relationship questions about it (contains, within, nearest) that no equality or range predicate can express, and indexing it with a bounding-box structure (an R-tree / GiST index) that supports a two-phase filter-and-refine query instead of the sorted-key lookup a B-tree gives you. Treat coordinates as two floats and everything silently breaks: distances come out in degrees, "within 5 km" scans the whole table, and a spatial join degrades to a cartesian product.

The three things that make spatial data different.

• Geometry types. A spatial value is one of Point, LineString, Polygon, or their Multi* collections, stored as structured geometry (typically WKB — well-known binary — under the hood, WKT for the human-readable form). A polygon has area and boundaries; a point has neighbours. None of that is expressible as a plain number, which is why spatial engines add a real geometry type rather than two double columns.
• CRS and projections. Every coordinate is meaningless without a coordinate reference system. SRID 4326 (WGS84) is geographic — degrees of latitude/longitude on an ellipsoid — where "distance" is along a curved surface. SRID 3857 (Web Mercator) is projected — metres on a flat plane — where distance is Euclidean but area is distorted away from the equator. Mixing SRIDs, or computing planar distance on 4326 degrees, is the single most common correctness bug in the field.
• Spatial indexing. A B-tree indexes a total order — great for =, <, BETWEEN, useless for "which shapes overlap this box." Spatial engines use an R-tree (in Postgres, a GiST index) that indexes each geometry's bounding box and supports the two-phase model below. The index is the difference between a millisecond point-in-polygon lookup and a full-table scan.

Filter-and-refine — the model every spatial engine uses.

• Phase 1 — filter (cheap, index). The spatial index compares bounding boxes only: it returns every geometry whose minimum bounding rectangle overlaps the query's box. This is fast and index-backed but approximate — a box can overlap while the exact shapes do not.
• Phase 2 — refine (exact, per candidate). The engine runs the exact geometric predicate (ST_Contains, ST_Intersects, ST_DWithin) only on the small candidate set from phase 1. Exact geometry math is expensive, so you want to run it on ten candidates, not ten million rows.
• Why it matters. Every performant spatial query — in PostGIS, in Sedona, in DuckDB — is this two-phase shape. When a spatial query is slow, the cause is almost always that phase 1 was skipped (no spatial index) so the exact predicate ran on the whole table.

The 2026 toolbox — four tools for four jobs.

• PostGIS is the transactional/serving spatial database: exact geometry, GiST indexes, the richest ST_ function library, ideal for point-in-polygon enrichment and serving spatial queries to an app.
• H3 is a hexagonal grid spatial index: it snaps each point to a cell id so aggregation becomes a GROUP BY on a string key — the workhorse for heatmaps, demand grids, and joining datasets by shared cell rather than by expensive geometry.
• GeoParquet is the storage/interchange format: geometry stored columnar inside Parquet with the CRS in file metadata, so lakehouse engines and GIS tools read the same file with predicate pushdown.
• Apache Sedona is the distributed engine: it brings spatial types, indexes, and partitioned spatial join to Spark so a job over billions of geometries runs across a cluster instead of dying on one node.

What interviewers listen for.

• Do you say a location is a geometry, not two numbers, and name the geometry types unprompted? — senior signal.
• Do you name the CRS/SRID and know 4326 is degrees while 3857/UTM is metres — and that distance must be computed in a metric CRS or with geography? — required answer.
• Do you explain a spatial index as a bounding-box R-tree and the filter-and-refine two-phase model? — senior signal.
• Do you pick a grid (H3) over exact geometry when the question is aggregation, not exact containment? — senior signal.
• Do you know when one machine is not enough and reach for a distributed spatial engine with spatial partitioning? — required answer.

Worked example — the CRS mistake that returns wrong distances

Detailed explanation. The most common geospatial bug is computing distance on geographic coordinates as if they were planar. A senior engineer states the CRS first, then either transforms to a metric CRS or uses the geography type so distances come back in metres. Walk the wrong-versus-right calculation for "how far apart are two points."

• The trap. ST_Distance on two geometry(Point, 4326) values returns a number in degrees, which is not a distance anyone can use (a degree of longitude is ~111 km at the equator and ~0 km at the poles).
• Fix A — transform. ST_Transform the geometry to a metric CRS (UTM or 3857) and measure there.
• Fix B — geography. Cast to geography, whose functions compute true distances on the ellipsoid in metres.

Question. Compute the distance between two coordinates correctly, and show why the naive geometry-in-degrees version is wrong.

Input.

Approach	Type / CRS	Unit returned	Correct?
ST_Distance(g1, g2) on 4326 geometry	geometry, degrees	degrees	no
ST_Distance(ST_Transform(...,3857), ...)	geometry, metres	metres (distorted)	approx
ST_Distance(g1::geography, g2::geography)	geography	metres (ellipsoid)	yes

Code.

-- Two points near London and Paris, stored as geographic coordinates (SRID 4326).
WITH pts AS (
  SELECT ST_SetSRID(ST_MakePoint(-0.1278, 51.5074), 4326) AS london,
         ST_SetSRID(ST_MakePoint( 2.3522, 48.8566), 4326) AS paris
)
SELECT
  -- WRONG: distance in DEGREES — a meaningless "3.5" that is not kilometres.
  ST_Distance(london, paris)                                   AS degrees_wrong,
  -- RIGHT (A): cast to geography → true metres on the ellipsoid.
  ST_Distance(london::geography, paris::geography)             AS metres_geography,
  -- RIGHT (B): project to a metric CRS first, then measure in metres.
  ST_Distance(ST_Transform(london, 3857), ST_Transform(paris, 3857)) AS metres_projected
FROM pts;

Step-by-step explanation.

1. Both points are built with ST_SetSRID(..., 4326) so the database knows they are longitude/latitude in degrees — attaching the SRID is step zero; a geometry with SRID 0 (unknown) is a latent bug because functions cannot reason about its units.
2. ST_Distance(london, paris) on plain geometry computes a planar Cartesian distance in the coordinate units, which for 4326 are degrees. The ~3.5 it returns is a straight-line distance in degree-space — not kilometres, and not even proportional to real distance across large spans.
3. Casting to ::geography switches to ellipsoidal math: ST_Distance now returns the true great-circle distance in metres (~343 km London–Paris) — the correct answer with no manual projection.
4. ST_Transform(..., 3857) reprojects the points into Web Mercator metres and measures there; it is close but distorted (Mercator inflates distance away from the equator), so for accuracy you use geography or a local metric CRS like the right UTM zone.
5. The lesson: never compute distance or area on 4326 as if it were planar. State the CRS, then either use geography (easiest, accurate globally) or project to an appropriate metric CRS (fastest for bulk planar work in a known region).

Output.

Column	Value	Meaning
degrees_wrong	~3.5	degrees in coordinate space — unusable
metres_geography	~343,000	true ellipsoidal metres — correct
metres_projected	~455,000	Mercator metres — distorted at this latitude
takeaway	—	measure in geography or a metric CRS, never raw 4326

Rule of thumb. A coordinate is nothing without its CRS. Store data in 4326 for interchange, but compute distance and area either by casting to geography (accurate metres on the ellipsoid) or by ST_Transform-ing into a local metric CRS — never by running planar math on degrees.

Worked example — filter-and-refine: why a spatial index is a bounding-box index

Detailed explanation. Understanding why a spatial index makes a query fast is the difference between an engineer who adds the right index and one who is puzzled that theirs did nothing. A spatial index stores bounding boxes and answers phase 1; the exact ST_ predicate is phase 2. Walk a "points inside a polygon" query through both phases.

• Phase 1 — the && bounding-box overlap. PostGIS's && operator asks "do the bounding boxes overlap?" and is answered directly by the GiST index — cheap, approximate.
• Phase 2 — the exact predicate. ST_Contains(poly, pt) runs the true geometric test, but only on the candidates phase 1 returned.
• The key insight. ST_Contains and friends automatically use the index by adding an internal && filter — but only if a GiST index exists. No index means no phase 1, so the exact test runs on every row.

Question. Explain how an indexed point-in-polygon query narrows millions of points to a handful before running exact geometry, and what happens without the index.

Input.

Phase	Operator	Cost	Backed by
1 filter	&& (bbox overlap)	cheap, approximate	GiST index
2 refine	ST_Contains (exact)	expensive, exact	run on candidates only
no index	ST_Contains on all rows	expensive × N	full scan

Code.

-- The GiST spatial index is what enables phase 1 (bounding-box filter).
CREATE INDEX pings_geom_gix ON pings USING GIST (geom);
CREATE INDEX zones_geom_gix ON zones USING GIST (geom);

-- Indexed point-in-polygon: ST_Contains INTERNALLY adds a && bbox filter,
-- so GiST returns candidate points, then exact geometry runs on just those.
EXPLAIN ANALYZE
SELECT z.zone_id, count(*) AS pings
FROM zones z
JOIN pings p
  ON z.geom && p.geom              -- phase 1: bounding-box overlap (index)
 AND ST_Contains(z.geom, p.geom)   -- phase 2: exact containment (candidates only)
GROUP BY z.zone_id;

Step-by-step explanation.

1. CREATE INDEX ... USING GIST (geom) builds an R-tree over each geometry's bounding box, not its exact shape — the index knows "zone 7 spans this rectangle," which is enough to exclude points nowhere near it.
2. z.geom && p.geom is phase 1: the GiST index answers "which point boxes overlap which zone boxes" without touching exact geometry, cutting millions of points to the few dozen plausibly inside each zone.
3. ST_Contains(z.geom, p.geom) is phase 2: it runs the exact (expensive) polygon test, but only on the candidate pairs phase 1 produced — so exact math runs thousands of times, not billions.
4. In practice you write only ST_Contains(...); PostGIS's planner injects the && filter automatically when a GiST index exists. Writing the && explicitly here just makes the two phases visible — the optimiser does it for you.
5. Drop the GiST index and phase 1 vanishes: ST_Contains must run on the full cartesian product of zones × points, turning a sub-second query into minutes or hours. This is why "my spatial query is slow" is almost always "I have no spatial index."

Output.

Scenario	Phase 1 candidates	Exact tests run	Latency
GiST index present	few per zone	thousands	milliseconds–seconds
No index	none (skipped)	billions (N×M)	minutes–hours
index + selective bbox	tiny	minimal	fastest
takeaway	bbox filter first	refine on candidates	index is mandatory

Rule of thumb. Every fast spatial query is filter-and-refine: a bounding-box index (GiST/R-tree) narrows candidates, then the exact ST_ predicate refines them. Always create a GiST index on any geometry column you filter or join on — without it, the exact test runs on every row and the query collapses to a scan.

Worked example — choosing the tool: PostGIS vs H3 vs GeoParquet vs Sedona

Detailed explanation. A senior geospatial interview always tests whether you reach for the right tool. The four in this guide are not competitors — they solve different jobs, and naming the job picks the tool. Build the decision table for a mobility company's location stack.

• The jobs. Serve exact spatial queries to an app; aggregate billions of points into a grid; store geometry for many engines to read; run a spatial join too big for one machine.
• The mistake. Using one tool for all four — e.g. forcing every aggregation through PostGIS polygon joins, or trying to serve low-latency app queries from a Spark job.
• The rule. Match the tool to whether the work is exact-geometry serving, grid aggregation, storage/interchange, or distributed compute.

Question. For four common geospatial jobs, name the right tool and the reason, and the wrong tool people often reach for.

Input.

Job	Right tool	Why
Exact point-in-polygon, served to an app	PostGIS	geometry + GiST + ST_ functions, low latency
Aggregate billions of points into a grid	H3	snap to cell → GROUP BY, no polygon join
Store geometry for many engines	GeoParquet	columnar, CRS in metadata, pushdown
Spatial join bigger than one machine	Apache Sedona	spatial partitioning + distributed join

Code.

Geospatial tool selection — pick by the JOB, not by habit.

Q: "Assign each live GPS ping to a delivery zone, serve to the dispatch app"
   -> PostGIS: geometry(zone) + GiST + ST_Contains join. Exact, indexed, low-latency.

Q: "Count rides per ~1 km cell across a year of pings for a heatmap"
   -> H3: latLngToCell(res 8) then GROUP BY cell. No polygon join, trivially parallel.

Q: "Publish the zones + a billion pings so DuckDB, GeoPandas AND Spark can read them"
   -> GeoParquet: columnar geometry (WKB) + CRS metadata + bbox pushdown. Engine-neutral.

Q: "Join 5B pings against 3M polygons nightly — one Postgres box can't"
   -> Apache Sedona: read GeoParquet, spatially partition both sides, distributed range join.

Anti-patterns:
  - PostGIS polygon join for a grid heatmap        -> use H3 (aggregation, not containment)
  - Spark job to serve a single app map query      -> use PostGIS (latency, not scale)
  - Shapefile/GeoJSON as the lake storage format   -> use GeoParquet (columnar, typed, big)
  - H3 where you need EXACT geometry (legal boundary) -> use PostGIS (a grid approximates)

Step-by-step explanation.

1. Exact, low-latency, served queries are PostGIS's home: a geometry column, a GiST index, and ST_Contains/ST_DWithin give millisecond point-in-polygon and radius answers a dispatch app can call per request.
2. Grid aggregation is H3's home: snapping every point to a hexagon id turns "how many here" into a GROUP BY cell that needs no geometry at query time and parallelises perfectly — the wrong shape for PostGIS, whose polygon joins are overkill for a regular grid.
3. Storage-for-many-readers is GeoParquet's home: geometry lives columnar with the CRS in the file, so DuckDB, GeoPandas, and Sedona all read it with predicate pushdown — a job neither a database nor a grid solves.
4. Bigger-than-one-machine spatial joins are Sedona's home: it spatially co-partitions both sides across a cluster and joins locally, which is the only one of the four that scales a raw geometry join horizontally.
5. The senior move is naming the job first — serve, aggregate, store, or distribute — and only then the tool; a real mobility stack uses all four together, GeoParquet feeding Sedona, H3 for grids, PostGIS for serving.

Output.

Job	Reach for	Not
Serve exact spatial queries	PostGIS	Sedona (latency)
Grid / heatmap aggregation	H3	PostGIS polygon join
Multi-engine geometry storage	GeoParquet	shapefile / GeoJSON
Distributed spatial join	Apache Sedona	single-node Postgres

Rule of thumb. The four tools are complementary, not competing: PostGIS serves exact geometry, H3 aggregates on a grid, GeoParquet stores for every engine, and Sedona distributes the heavy joins. Name the job — serve, aggregate, store, or distribute — and the tool picks itself; a mature stack runs all four.

Senior interview question on designing a geospatial pipeline

A senior interviewer often opens with: "A mobility company streams five billion GPS pings a day and owns three million delivery-zone polygons. They need three things: a live dispatch app that maps each ping to its zone in milliseconds, a daily heatmap of demand at roughly one-kilometre resolution, and a nightly enrichment that tags every ping with its zone for the warehouse. Design the geospatial data engineering stack: what geometry and CRS you store, where each spatial index lives, and which tool does each job — serving, aggregation, storage, and the distributed join — and why one database cannot do all of it."

Solution Using PostGIS serving, H3 aggregation, GeoParquet storage, and Sedona for the distributed join

-- Step 1 — PostGIS serves the LIVE dispatch query: zones as indexed geometry.
CREATE TABLE zones (
  zone_id int PRIMARY KEY,
  name    text,
  geom    geometry(MultiPolygon, 4326)     -- geometry + explicit CRS
);
CREATE INDEX zones_gix ON zones USING GIST (geom);   -- R-tree for filter-and-refine

-- Live: which zone contains this ping? Indexed point-in-polygon, milliseconds.
SELECT zone_id FROM zones
WHERE ST_Contains(geom, ST_SetSRID(ST_MakePoint(:lng, :lat), 4326))
LIMIT 1;

# Step 2 — H3 powers the daily heatmap: snap each ping to a hex, then GROUP BY.
import h3
res = 8                                   # ~0.7 km edge cells
df["cell"] = [h3.latlng_to_cell(lat, lng, res) for lat, lng in pings]
heatmap = df.groupby("cell").size()       # demand per cell — no polygon join

# Step 3 — GeoParquet is the lake storage every engine reads (pings + zones).
import geopandas as gpd
zones_gdf.to_parquet("s3://lake/zones.parquet")          # WKB geometry + CRS metadata
pings_gdf.to_parquet("s3://lake/pings/date=2026-08-26/") # partitioned, columnar

# Step 4 — Sedona runs the NIGHTLY distributed enrichment join (5B x 3M).
sedona.read.format("geoparquet").load("s3://lake/pings/date=2026-08-26/") \
      .createOrReplaceTempView("pings")
sedona.read.format("geoparquet").load("s3://lake/zones.parquet") \
      .createOrReplaceTempView("zones")
enriched = sedona.sql("""
  SELECT /*+ BROADCAST(z) */ p.ping_id, z.zone_id       -- broadcast the small side
  FROM pings p JOIN zones z
    ON ST_Contains(z.geom, p.geom)                       -- distributed range join
""")

Step-by-step trace.

Job	Tool	Where the spatial index lives
Live dispatch (ms)	PostGIS	GiST R-tree on zones.geom
Daily heatmap	H3	the hex cell id is the index (a group key)
Lake storage	GeoParquet	per-row-group bbox in the file metadata
Nightly enrich (5B×3M)	Apache Sedona	spatial partitioning + per-partition local index
CRS	4326 stored, metric for distance	attached everywhere, never assumed
Failure avoided	one tool for all	each job on the tool built for it

After the design, the dispatch app hits PostGIS with an indexed ST_Contains and gets a zone in milliseconds; the heatmap is an H3 GROUP BY cell that needs no geometry at query time and parallelises trivially; both pings and zones live in GeoParquet so every engine reads the same CRS-tagged columnar files; and the nightly five-billion-by-three-million enrichment runs in Sedona, which broadcasts the small zone side and spatially partitions the pings so the join is local per partition — a job no single Postgres box could finish. One tool per job, all sharing 4326 geometry.

Output:

Metric	One-database attempt	Right-tool stack
Dispatch query latency	ms (fine in PostGIS)	ms (PostGIS)
Heatmap over a year	slow polygon joins	fast H3 group-by
Storage read by all engines	DB export per engine	one GeoParquet file
5B×3M nightly join	impossible on one box	distributed (Sedona)
CRS correctness	ad-hoc, error-prone	explicit everywhere

Why this works — concept by concept:

• PostGIS for serving — a geometry(_, 4326) column plus a GiST index gives exact, indexed point-in-polygon at app latency, which is precisely the serving job a grid or a batch engine cannot do in milliseconds.
• H3 for aggregation — snapping points to hexagon ids replaces expensive polygon joins with a GROUP BY on a string key, so a year-long heatmap becomes a trivially parallel counting problem instead of billions of containment tests.
• GeoParquet for storage — columnar geometry with the CRS in file metadata makes one file readable by GeoPandas, DuckDB, and Sedona with bbox pushdown, so you store once and every engine reads it natively.
• Sedona for the distributed join — spatial partitioning plus a broadcast of the small side turns a five-billion-by-three-million join into local per-partition work across a cluster, the only way that join finishes at all.
• Cost — one indexed lookup per dispatch request, one GROUP BY per heatmap, one columnar file per dataset, and a partitioned nightly join, versus forcing everything through a single database that melts on the batch join and does slow group-bys. The eliminated cost is the outage and bill of running batch-scale geometry on a serving box — O(log n) indexed serving and O(n) partitioned batch, each on the right engine.

Design
Topic — design
Design problems on geospatial pipelines and system architecture

Practice →

Indexing Topic — indexing Indexing problems on spatial and bounding-box indexes

Practice →

---

2. PostGIS — spatial SQL over Postgres

A geometry column plus a GiST index turns Postgres into a spatial engine

The mental model in one line: PostGIS is the extension that gives Postgres a real geometry/geography type, a GiST **spatial index that indexes each shape's bounding box, and a library of ST_ functions — so a location becomes a first-class column you can ST_Contains, ST_DWithin, and ST_Intersects, and a spatial join (which points fall in which polygons, which stores are within five kilometres) becomes ordinary indexed SQL that runs filter-and-refine automatically.** Model the geometry with the right type and SRID, add a GiST index, and PostGIS serves exact spatial queries at app latency; skip the index or mix SRIDs and the same queries scan the whole table or return nonsense.

geometry vs geography — pick by how you measure.

• geometry. Planar Cartesian math in the column's SRID. Fast, supports the full ST_ library, and correct if the SRID is a metric projection (UTM, 3857) for the region — the default for bulk analytical work.
• geography. Ellipsoidal math on 4326 lon/lat, returning true metres for distance and area globally with no projection step. Slightly slower and a smaller function set, but the right choice when data spans the globe or you want correct distances without choosing a projection.
• The SRID belongs on the column. geometry(Point, 4326) declares both the type and the CRS; the database rejects inserts in the wrong SRID and functions know the units. A bare geometry with SRID 0 is a bug waiting to happen.

The GiST spatial index — the filter half.

• What it indexes. A GiST index over a geometry column stores each shape's bounding box in an R-tree, so overlap queries are index scans, not table scans.
• The && operator. a && b means "bounding boxes overlap" and is the phase-1 filter every spatial predicate uses under the hood.
• Automatic use. ST_Contains, ST_Intersects, and ST_DWithin add the && bbox filter internally, so they use the GiST index automatically — if it exists. The index is mandatory on any geometry you filter or join.

The ST_ function family — the refine half.

• Relationships. ST_Contains, ST_Within, ST_Intersects, ST_Covers — the exact topological tests, run on the candidates the index returned.
• Distance. ST_DWithin(a, b, d) is the indexed "within distance d" predicate (it uses the bbox expanded by d); ST_Distance computes the actual distance. Prefer ST_DWithin for filtering because it uses the index; ST_Distance in a WHERE does not.
• Measures and transforms. ST_Area, ST_Length, ST_Transform (reproject), ST_Buffer, ST_Centroid — the everyday toolkit for computing on and reshaping geometry.

The failure modes senior engineers pre-empt.

• No GiST index. The single biggest PostGIS performance bug: without the index there is no phase-1 filter, so every spatial predicate scans the table. Mitigation: CREATE INDEX ... USING GIST (geom) on every filtered/joined geometry.
• Mixed or missing SRIDs. Operating on two different SRIDs raises an error or silently misbehaves; SRID 0 means "unknown units." Mitigation: declare the SRID on the column and ST_Transform before combining.
• ST_Distance in a WHERE. WHERE ST_Distance(a,b) < 5000 computes distance for every row and cannot use the index. Mitigation: WHERE ST_DWithin(a, b, 5000), which is index-backed.

Common interview probes on PostGIS.

• "How does a spatial index work?" — GiST/R-tree over bounding boxes; phase-1 && filter, then exact ST_ refine.
• "geometry or geography?" — geometry (planar, fast) in a metric SRID; geography (ellipsoidal metres) for global data or easy correct distances.
• "How do you find everything within 5 km?" — ST_DWithin(a, b, 5000) (index-backed), not ST_Distance < 5000.
• "Why is my spatial query slow?" — almost always no GiST index, so the exact predicate scans every row.

Worked example — a geometry column, a GiST index, and a point-in-polygon query

Detailed explanation. The canonical PostGIS setup: a table of polygons and a table of points, each geometry indexed with GiST, then an indexed ST_Contains join that assigns every point to its containing polygon. Assign delivery pings to zones.

• The tables. zones(geom MultiPolygon), pings(geom Point), both SRID 4326.
• The indexes. A GiST index on each geometry column.
• The query. ST_Contains(zone, ping) — indexed point-in-polygon.

Question. Model zones and pings, index both, and count how many pings fall inside each zone with an indexed spatial join.

Input.

Piece	Value
Zones	geometry(MultiPolygon, 4326) + GiST
Pings	geometry(Point, 4326) + GiST
Predicate	ST_Contains(zone.geom, ping.geom)
Output	pings per zone

Code.

-- 1. Model geometry with an explicit type and SRID (units are known).
CREATE TABLE zones (
  zone_id int PRIMARY KEY,
  name    text,
  geom    geometry(MultiPolygon, 4326)
);
CREATE TABLE pings (
  ping_id bigint PRIMARY KEY,
  geom    geometry(Point, 4326)
);

-- 2. A GiST spatial index on each geometry column → enables the phase-1 bbox filter.
CREATE INDEX zones_gix ON zones USING GIST (geom);
CREATE INDEX pings_gix ON pings USING GIST (geom);

-- 3. Indexed point-in-polygon spatial join: ST_Contains auto-adds the && bbox filter.
SELECT z.zone_id, z.name, count(p.ping_id) AS pings
FROM zones z
LEFT JOIN pings p
  ON ST_Contains(z.geom, p.geom)   -- GiST filter, then exact containment per candidate
GROUP BY z.zone_id, z.name
ORDER BY pings DESC;

Step-by-step explanation.

1. Each geom column declares its type and SRID 4326, so the planner and every ST_ function know the coordinates are lon/lat degrees — the SRID is metadata the query engine actually uses, not decoration.
2. CREATE INDEX ... USING GIST (geom) builds the R-tree of bounding boxes on both tables; this is the precondition for the join to be anything other than a full cartesian scan.
3. ST_Contains(z.geom, p.geom) reads as the exact test, but PostGIS injects an internal z.geom && p.geom bbox filter first, so GiST prunes each zone to the handful of pings whose boxes overlap it before the expensive polygon test runs.
4. The LEFT JOIN keeps zones with zero pings (they show 0), and GROUP BY counts matches per zone — an ordinary aggregate over what is, mechanically, an indexed spatial join.
5. Remove either GiST index and the planner falls back to testing every zone against every ping — the same query, correct results, but orders of magnitude slower. The index is what makes it a serving-grade query.

Output.

zone_id	name	pings
12	Downtown	48,201
7	Harbour	21,940
3	Airport	15,002
21	Suburb North	0

Rule of thumb. Model geometry with an explicit type and SRID, put a GiST index on every geometry you filter or join, and write the exact predicate (ST_Contains) — PostGIS adds the bounding-box filter for you. The index turns a point-in-polygon join from a scan into a millisecond lookup.

Worked example — a within-radius spatial join with geography for correct metres

Detailed explanation. "Find everything within N metres" is a spatial join people get subtly wrong: they use ST_Distance < N (no index, and in degrees if the type is geometry-4326). The correct, indexed, metre-accurate form is ST_DWithin on geography. Find every customer within 5 km of each store.

• The trap. WHERE ST_Distance(store, customer) < 5000 — scans every pair and, on 4326 geometry, compares degrees to 5000.
• The fix. ST_DWithin(store::geography, customer::geography, 5000) — index-backed and in true metres.
• The index. A GiST index (usable on the geography cast) provides the phase-1 filter.

Question. For each store, list customers within five kilometres, using an indexed, metre-accurate radius join.

Input.

Aspect	Wrong	Right
Predicate	ST_Distance(a,b) < 5000	ST_DWithin(a,b,5000)
Units on 4326 geometry	degrees	metres (via geography)
Uses index	no (scan)	yes (bbox+d filter)
Cost	O(N×M)	O(candidates)

Code.

-- Store the point as geometry(4326) but a GiST index also serves the geography cast.
CREATE INDEX stores_gix    ON stores    USING GIST (geom);
CREATE INDEX customers_gix ON customers USING GIST (geom);

-- Within-radius spatial join: ST_DWithin is INDEX-BACKED and, on geography, in METRES.
SELECT s.store_id, c.customer_id,
       ST_Distance(s.geom::geography, c.geom::geography) AS metres
FROM stores s
JOIN customers c
  ON ST_DWithin(s.geom::geography, c.geom::geography, 5000)   -- <= 5 km, indexed
ORDER BY s.store_id, metres;

Step-by-step explanation.

1. ST_DWithin(a, b, 5000) is the indexed radius predicate: internally it expands the bounding box by the distance and uses GiST for phase 1, so only nearby candidates reach the exact distance test — unlike ST_Distance(...) < 5000, which computes distance for every pair.
2. Casting both operands to ::geography makes the 5000 mean 5000 metres on the ellipsoid, the correct interpretation; the same call on raw geometry(4326) would treat 5000 as degrees, matching essentially the whole planet.
3. GiST on the geometry column still accelerates the geography predicate here (PostGIS can use the bbox), so the join stays index-backed rather than degrading to a scan.
4. The selected ST_Distance(...::geography) reports the actual metres for each surviving pair, so you get both the filter (within 5 km) and the measurement (how far) in one pass.
5. The ORDER BY s.store_id, metres yields, per store, its customers nearest-first — the shape a "stores and their nearby customers" feature needs, produced by one indexed radius join instead of an application loop computing distances.

Output.

store_id	customer_id	metres
1	5567	412
1	9021	3,880
2	3345	1,204
2	7789	4,950

Rule of thumb. For radius queries use ST_DWithin(a, b, metres) on geography — it is index-backed and metre-accurate — never ST_Distance(...) < metres, which scans every pair and, on 4326 geometry, silently compares degrees. Filter with ST_DWithin; measure the survivors with ST_Distance.

Worked example — nearest-neighbour with the KNN <-> operator

Detailed explanation. "Which store is nearest to each customer" is a k-nearest-neighbour query, and PostGIS has an index-assisted operator for it: <-> (distance-ordered) drives a GiST index scan that returns rows in distance order without computing every distance. Use it with ORDER BY ... LIMIT for true index-backed nearest-neighbour.

• The naive way. Cross join, ST_Distance, ORDER BY, LIMIT 1 — computes every pair.
• The KNN way. ORDER BY store.geom <-> customer.geom LIMIT k — GiST walks the index nearest-first.
• The pattern. A LATERAL join runs the per-customer KNN lookup against the indexed store table.

Question. For each customer, find the single nearest store using the index-assisted KNN operator rather than an all-pairs distance computation.

Input.

Approach	Mechanism	Cost
Cross join + ST_Distance	compute all pairs	O(N×M)
<-> + ORDER BY + LIMIT	GiST nearest-first walk	O(N log M)
LATERAL per customer	index scan per row	O(N log M)

Code.

CREATE INDEX stores_gix ON stores USING GIST (geom);   -- KNN needs the GiST index

-- For each customer, the nearest store: <-> drives a distance-ordered index scan.
SELECT c.customer_id, s.store_id,
       ST_Distance(c.geom::geography, s.geom::geography) AS metres
FROM customers c
CROSS JOIN LATERAL (
  SELECT store_id, geom
  FROM stores
  ORDER BY stores.geom <-> c.geom   -- KNN: index returns rows nearest-first
  LIMIT 1                            -- just the closest (use k for k-nearest)
) s
ORDER BY c.customer_id;

Step-by-step explanation.

1. stores.geom <-> c.geom is the KNN distance operator; combined with ORDER BY ... LIMIT, GiST walks the R-tree returning stores in increasing distance from the customer, so it can stop after k without scoring the rest.
2. The LATERAL subquery runs that nearest-first index scan once per customer, correlating on c.geom — this is the mechanism that keeps a "nearest store to every customer" query O(N log M) instead of O(N×M).
3. LIMIT 1 returns the single closest store; changing it to LIMIT k yields the k nearest, still index-driven — the same operator serves nearest-one and nearest-k.
4. The outer ST_Distance(...::geography) reports the actual metres, computed only for the surviving nearest rows, not for every pair — measurement follows the index-narrowed selection.
5. Note the ordering operator uses planar <-> on the geometry for the index walk while distance is reported in geography metres; for high-accuracy nearest-neighbour across large areas you project to a local metric CRS so the index ordering matches true distance.

Output.

customer_id	store_id	metres
101	4	233
102	1	1,904
103	4	508
104	7	2,220

Rule of thumb. For nearest-neighbour, use the <-> KNN operator with ORDER BY ... LIMIT k and a GiST index — it walks the index nearest-first and stops early, turning an all-pairs distance computation into a per-row index scan. Wrap it in a LATERAL join to run one KNN lookup per driving row.

Senior interview question on PostGIS spatial joins at scale

A senior interviewer might ask: "You have three million delivery-zone polygons and a table of tens of millions of GPS pings in Postgres, and you need to enrich every ping with its containing zone and also answer 'stores within 5 km' for a live feature. Design it in PostGIS: the geometry types and CRS you choose, the indexes you build, how the point-in-polygon enrichment stays fast, how the radius query stays index-backed and metre-accurate, and the mistakes that would turn either query into a full-table scan."

Solution Using geometry columns, GiST indexes, indexed ST_Contains, and ST_DWithin on geography

-- 1. Geometry with explicit type + CRS; GiST index on every column we filter/join.
CREATE TABLE zones  (zone_id int PRIMARY KEY, geom geometry(MultiPolygon, 4326));
CREATE TABLE pings  (ping_id bigint PRIMARY KEY, geom geometry(Point, 4326), zone_id int);
CREATE TABLE stores (store_id int PRIMARY KEY, geom geometry(Point, 4326));

CREATE INDEX zones_gix  ON zones  USING GIST (geom);
CREATE INDEX pings_gix  ON pings  USING GIST (geom);
CREATE INDEX stores_gix ON stores USING GIST (geom);

-- 2. Point-in-polygon enrichment: indexed ST_Contains assigns each ping a zone.
UPDATE pings p
SET zone_id = z.zone_id
FROM zones z
WHERE ST_Contains(z.geom, p.geom)          -- GiST filter + exact refine, per ping
  AND p.zone_id IS DISTINCT FROM z.zone_id; -- only touch rows that change

-- 3. Live "stores within 5 km": index-backed + metre-accurate via geography.
SELECT s.store_id,
       ST_Distance(s.geom::geography, ST_SetSRID(ST_MakePoint(:lng,:lat),4326)::geography) AS m
FROM stores s
WHERE ST_DWithin(s.geom::geography,
                 ST_SetSRID(ST_MakePoint(:lng,:lat),4326)::geography, 5000)
ORDER BY m;

-- 4. Prove the plan uses the index (the difference between ms and a full scan).
EXPLAIN (ANALYZE, BUFFERS)
SELECT count(*) FROM zones z JOIN pings p ON ST_Contains(z.geom, p.geom);
-- Expect: "Index Scan using pings_gix" / a nested loop over GiST candidates,
--         NOT a "Seq Scan" on pings.

Step-by-step trace.

Concern	Component	Effect
Units known	geometry(_, 4326)	functions know it is lon/lat degrees
Phase-1 filter	GiST on every geom	bbox candidates, not full scan
Enrichment	indexed ST_Contains	each ping tagged with its zone
Radius, correct	ST_DWithin on geography	index-backed, true metres
Verify	EXPLAIN ANALYZE	confirms index scan, not seq scan
Idempotence	IS DISTINCT FROM	re-runs touch only changed rows

After deployment, every geometry column has a GiST index, so the enrichment UPDATE runs ST_Contains only on bbox-overlapping ping/zone pairs and tags tens of millions of pings in one pass; the live radius query uses ST_DWithin on geography, so it is both index-backed and measured in real metres; and EXPLAIN ANALYZE confirms an index scan rather than a sequential scan — the check that separates a millisecond query from a table scan. The two classic bugs (no GiST index, ST_Distance in the WHERE) are structurally avoided.

Output:

Metric	Naive (no index / ST_Distance)	PostGIS (GiST + ST_DWithin)
Enrichment of N pings	O(N×M) scan	O(N log M) indexed
Radius query	full scan, degrees	index scan, metres
Correctness of distance	wrong (degrees)	correct (geography metres)
Plan	Seq Scan	Index Scan (verified)
Re-run cost	rewrites all rows	only changed rows

Why this works — concept by concept:

• Typed geometry with an SRID — declaring geometry(_, 4326) attaches the CRS to the column, so functions know the units and the database rejects mismatched inserts, eliminating the silent degrees-versus-metres class of bug.
• GiST index for filter-and-refine — the R-tree over bounding boxes provides the phase-1 filter that ST_Contains/ST_DWithin use automatically, turning both the enrichment join and the radius query from scans into indexed candidate lookups.
• ST_DWithin on geography — the radius predicate is both index-backed (bbox expanded by the distance) and metre-accurate (ellipsoidal), the two properties ST_Distance < d on raw geometry fails to provide at once.
• EXPLAIN as the proof — reading the plan for Index Scan versus Seq Scan is how a senior engineer verifies the spatial index is used rather than assuming it, catching the missing-index regression before production does.
• Cost — indexed candidate lookups and a metre-correct radius filter, versus an all-pairs scan computing wrong-unit distances. The eliminated cost is the full-table scan a missing GiST index forces — O(log n) index navigation per probe instead of O(n) per query, on data sized for serving.

Indexing
Topic — indexing
Indexing problems on GiST and spatial index design

Practice →

Joins Topic — joins Join problems on spatial and point-in-polygon joins

Practice →

---

3. H3 — hexagonal hierarchical spatial indexing

Snap every point to a hexagon id, and spatial aggregation becomes an ordinary group-by

The mental model in one line: H3 is a hierarchical hexagonal grid that tiles the earth into cells at sixteen resolutions, where every location maps to a single 64-bit cell id — so instead of an expensive polygon **spatial join, you latLngToCell each point to its hexagon and GROUP BY the id, neighbours are a cheap gridDisk (k-ring) call, and changing resolution is a cellToParent/cellToChildren walk up and down the hierarchy — turning "how many events happened near here" from a geometry problem into a string-key aggregation that any engine can parallelise.** Hexagons beat the square cells of a geohash because every neighbour is equidistant and there is no distortion of adjacency, which is exactly what demand grids, heatmaps, and spatial rollups need.

What H3 is and why hexagons.

• A hierarchical hex grid. H3 covers the globe with hexagonal cells at resolutions 0 (coarse, ~1,000 km) through 15 (fine, ~0.5 m edge). Each cell has a stable id (a 64-bit integer, usually shown as a 15-character hex string like 8a2a1072b59ffff).
• Hexagons vs geohash squares. In a hex grid every cell has six neighbours all at the same centre-to-centre distance; a geohash/square grid has neighbours at two different distances (edge vs corner) and its rectangles distort with latitude. Uniform adjacency makes hexagons better for gradients, flow, and k-ring smoothing.
• The id is the index. Because the cell id encodes location, grouping by it is spatial indexing — no R-tree, no geometry at query time. That is what makes H3 aggregation trivially distributable.

The core operations.

• latLngToCell(lat, lng, res). Map a coordinate to its cell id at a resolution — the "index this point" call you run once per row at ingest.
• cellToLatLng / cellToBoundary. Recover a cell's centre or its hexagon polygon — for rendering a heatmap or exporting cells as geometry.
• gridDisk(cell, k) (k-ring). All cells within k steps of a cell — the neighbourhood operator for "demand within ~k cells," smoothing, and coverage.
• gridDistance / gridPathCells. Grid distance and path between cells — cheap integer operations on the grid, not geometry.

The hierarchy — change resolution without touching geometry.

• cellToParent(cell, res). The coarser cell that contains this one — roll a fine grid up to a coarse one by re-grouping on the parent id, no re-aggregation from raw points.
• cellToChildren(cell, res). The finer cells inside a coarse one — drill down.
• Compaction. compactCells replaces a full set of fine cells with the smallest mixed-resolution set covering the same area, shrinking storage for coverage/region sets.
• The caveat. A parent hex is approximately, not exactly, the union of its children (hexagons cannot perfectly nest), so H3 hierarchy is for aggregation and indexing, not for exact area subdivision.

The failure modes senior engineers pre-empt.

• Wrong resolution. Too coarse hides the signal; too fine explodes cardinality and storage. Mitigation: pick resolution from the analysis grain (city ≈ res 6–7, neighbourhood ≈ 8–9, building ≈ 11–12) and store the id, re-rolling up with cellToParent.
• Pentagons and edge effects. H3 has twelve pentagon cells per resolution (icosahedron vertices); most workloads never touch them, but distance/area math near them is special-cased. Mitigation: be aware for global metrics; ignore for city-scale work.
• Using a grid where exact geometry is required. A hex only approximates a real boundary (a legal district, a coastline). Mitigation: use H3 for aggregation/joins-by-cell and PostGIS geometry for exact containment.

Common interview probes on H3.

• "Why hexagons over a geohash?" — uniform equidistant neighbours, no latitude distortion, better for gradients and k-rings.
• "How does H3 make a spatial join cheap?" — both datasets get a cell id; you join/group on the id, a string equality, not geometry.
• "How do you change resolution?" — cellToParent/cellToChildren; roll up by re-grouping on the parent id.
• "When is H3 the wrong tool?" — when you need exact containment against a real polygon; a hex only approximates it.

Worked example — latLngToCell: index points to hexagons and aggregate

Detailed explanation. The core H3 pattern: assign each point a cell id at ingest, then aggregate by the id. Build a demand grid — rides per hexagon — from raw pickup coordinates, with no geometry at query time.

• The step. latlng_to_cell(lat, lng, res) per row → a cell id column.
• The aggregation. GROUP BY cell → count per hexagon.
• The render. cell_to_boundary(cell) only when you need the polygon for a map.

Question. Turn a table of ride pickups into a per-hexagon demand count at ~1 km resolution, then show how you would render the top cells.

Input.

Piece	Value
Input	pickups (lat, lng)
Resolution	8 (~0.46 km² cell)
Index op	latlng_to_cell
Aggregate	GROUP BY cell → count

Code.

import h3
import pandas as pd

pickups = pd.DataFrame({"lat": [...], "lng": [...]})   # raw ride pickups

# 1. INDEX: snap every pickup to its H3 cell id at resolution 8 (~0.46 km^2).
RES = 8
pickups["cell"] = [
    h3.latlng_to_cell(lat, lng, RES) for lat, lng in zip(pickups.lat, pickups.lng)
]

# 2. AGGREGATE: demand per hexagon is now an ordinary group-by on a STRING key.
demand = (pickups.groupby("cell").size()
          .reset_index(name="rides")
          .sort_values("rides", ascending=False))

# 3. RENDER (only when needed): recover each hot cell's hexagon polygon for a map.
demand["boundary"] = [h3.cell_to_boundary(c) for c in demand.cell]  # list of latlng

Step-by-step explanation.

1. h3.latlng_to_cell(lat, lng, 8) maps each pickup to the id of the resolution-8 hexagon that contains it — a pure function of the coordinate, computed once at ingest and stored as a plain string column.
2. Resolution 8 gives ~0.46 km² hexagons (~460 m edge), a good neighbourhood grain for city demand; choosing the resolution is choosing the analysis grain, and it is the one real modelling decision here.
3. groupby("cell").size() is the whole aggregation: because the cell id encodes location, counting rides per hexagon is a string-key group-by with no spatial predicate, no index, and no join — which is why it parallelises across any engine trivially.
4. The result is a demand grid keyed by hexagon; sorting by rides surfaces hotspots directly, and nothing in the query path ever touched a polygon.
5. cell_to_boundary(cell) is called only for the cells you actually draw, converting each hot cell id back into its hexagon vertices for the map — geometry re-enters only at render time, not during aggregation.

Output.

cell	rides
8a2a1072b59ffff	1,240
8a2a1072b5affff	980
8a2a1072b587fff	754
8a2a1072b5b7fff	611

Rule of thumb. Index each point with latlng_to_cell once at ingest, store the cell id, and let every "how many near here" question become a GROUP BY cell. Geometry only re-enters at render time via cell_to_boundary; choose the resolution to match the analysis grain, because that is the one decision that matters.

Worked example — gridDisk (k-ring) neighbourhood aggregation

Detailed explanation. Many questions are not "how many in this cell" but "how many in this cell and its neighbourhood" — coverage, smoothing, "demand around a driver." gridDisk(cell, k) returns every cell within k grid steps, so neighbourhood aggregation is a lookup plus a group-by. Compute demand within one ring of each candidate location.

• The op. grid_disk(cell, k) → the set of cells ≤ k steps away.
• The use. Sum a per-cell metric over each cell's k-ring for a smoothed surface.
• The cost. A k-ring has 3k(k+1)+1 cells — cheap integer grid math, no geometry.

Question. Given per-cell demand, compute a smoothed "demand within 1 ring" value for each cell using k-ring neighbourhoods.

Input.

Aspect	Value
Per-cell metric	rides by cell
Neighbourhood	grid_disk(cell, k=1)
Cells in a 1-ring	7 (centre + 6)
Output	smoothed demand per cell

Code.

import h3

# demand: dict cell -> rides (from the previous group-by)
def demand_in_ring(demand: dict, k: int = 1) -> dict:
    smoothed = {}
    for cell in demand:
        # k-ring: this cell plus everything within k grid steps (7 cells for k=1).
        ring = h3.grid_disk(cell, k)
        # Sum the metric over the neighbourhood — missing neighbours contribute 0.
        smoothed[cell] = sum(demand.get(c, 0) for c in ring)
    return smoothed

smoothed = demand_in_ring(demand_dict, k=1)

# "Is this cell a hotspot including its surroundings?" is now a dict lookup.
hottest = sorted(smoothed.items(), key=lambda kv: kv[1], reverse=True)[:5]

Step-by-step explanation.

1. h3.grid_disk(cell, 1) returns the seven cells of the 1-ring (the centre plus its six equidistant neighbours) — because hexagon neighbours are uniform, this ring is a clean, isotropic neighbourhood with no edge-versus-corner ambiguity a square grid would have.
2. Summing demand.get(c, 0) over the ring produces a smoothed value: a cell surrounded by busy cells scores high even if its own count is modest, which is exactly the coverage/heat signal you want for siting or dispatch.
3. demand.get(c, 0) treats absent neighbours as zero, so cells on the edge of the data are handled without special-casing — the grid math is defined everywhere.
4. A k-ring has 3k(k+1)+1 cells (7 for k=1, 19 for k=2), so widening the neighbourhood is a controlled, integer-bounded cost — no geometry, no distance computation, just id enumeration.
5. Because every operation is on cell ids, the whole smoothing pass is a hash-join/lookup problem that distributes across a cluster identically to the plain group-by — the neighbourhood analysis inherits H3's parallelism.

Output.

cell	own rides	ring demand (k=1)
8a2a1072b59ffff	1,240	3,910
8a2a1072b5affff	980	3,455
8a2a1072b587fff	754	2,102
8a2a1072b5b7fff	611	1,880

Rule of thumb. Use grid_disk(cell, k) to turn neighbourhood questions — coverage, smoothing, "demand around here" — into a k-ring enumeration plus a group-by on cell ids. Hexagons make the ring isotropic (all neighbours equidistant), and the whole computation stays integer grid math with no geometry, so it parallelises like any key aggregation.

Worked example — cellToParent: hierarchical rollup and compaction

Detailed explanation. A powerful H3 property: once points are indexed at a fine resolution, you can roll up to any coarser resolution without re-reading the raw data, by mapping each fine cell to its parent and re-grouping. And a coverage set of fine cells can be compacted into a smaller mixed-resolution set. Roll a res-9 demand grid up to res-7, and compact a coverage set.

• Rollup. cell_to_parent(cell, res) → coarse id; re-GROUP BY the parent.
• Drilldown. cell_to_children(cell, res) → the fine cells inside a coarse one.
• Compaction. compact_cells(cells) → smallest mixed-res set covering the same area.

Question. Aggregate a fine-resolution demand grid up to a coarser resolution by re-grouping on parent cells, and compact a coverage set for storage.

Input.

Operation	Call	Result
Rollup res 9 → 7	cell_to_parent(c, 7)	coarse cell id
Re-aggregate	GROUP BY parent	demand per coarse cell
Compact coverage	compact_cells(set)	fewer, mixed-res cells

Code.

import h3
import pandas as pd

# fine: DataFrame with res-9 `cell` and `rides` (already indexed + aggregated).

# 1. ROLL UP res 9 -> res 7 WITHOUT re-reading raw points: map to parent, re-group.
fine["parent7"] = [h3.cell_to_parent(c, 7) for c in fine.cell]
coarse = fine.groupby("parent7")["rides"].sum().reset_index(name="rides")

# 2. COMPACT a coverage set (e.g. all res-9 cells of a service area) for storage.
service_area_cells = set(fine.cell)                 # many uniform res-9 cells
compacted = h3.compact_cells(list(service_area_cells))  # fewer, mixed-resolution cells
print(len(service_area_cells), "->", len(compacted))    # e.g. 40,000 -> 6,200

Step-by-step explanation.

1. cell_to_parent(c, 7) returns the resolution-7 cell that (approximately) contains the resolution-9 cell — a pure id-to-id function, so the rollup needs only the existing aggregated grid, never the billions of raw points behind it.
2. Re-grouping on parent7 and summing rides produces the coarse demand grid; because aggregation is associative over the hierarchy, summing children into parents is exact for counts even though the hexagons only approximately nest.
3. This is the operational win: you index once at a fine resolution and serve every coarser zoom level by cheap re-grouping, instead of re-scanning source data per zoom — the pattern behind interactive multi-resolution maps.
4. compact_cells replaces a large set of same-resolution cells with the smallest mixed-resolution set covering the same area (a full parent's worth of children collapses to the single parent), shrinking a coverage/region definition dramatically for storage and membership tests.
5. The caveat stays in view: hierarchy and compaction are exact for set coverage and additive counts, but a parent hexagon is not geometrically the exact union of its children, so you do not use this for precise area math — that is PostGIS's job.

Output.

Operation	Before	After
Rollup res 9 → 7	fine grid (many cells)	coarse grid (fewer)
Re-read raw points	not needed	not needed
Coverage set size	40,000 res-9 cells	6,200 mixed-res cells
Count correctness	—	exact (additive)

Rule of thumb. Index at a fine resolution once, then serve every coarser view by mapping to cell_to_parent and re-grouping — never re-scan raw data per zoom level. Use compact_cells to shrink coverage sets to a mixed-resolution minimum, remembering that H3 hierarchy is exact for counts and coverage but only approximate for geometry.

Senior interview question on H3 for large-scale spatial aggregation

A senior interviewer might ask: "You need a demand heatmap over a year of pings — tens of billions of points — at neighbourhood resolution, with the ability to zoom out to city level, and you want it to join cleanly against a weather dataset that is also point data. A PostGIS polygon-join approach is too slow. Design an H3-based pipeline: how you index the points, why aggregation is now a group-by, how you serve multiple zoom levels, and how joining two datasets by cell replaces a geometry join."

Solution Using latLngToCell indexing, group-by aggregation, cellToParent rollup, and join-by-cell

# 1. INDEX at ingest: every ping gets a res-9 cell id (a plain string column).
import h3
RES_FINE = 9
pings = pings.withColumn(
    "cell", h3_latlng_to_cell(pings.lat, pings.lng, RES_FINE)   # UDF / native H3
)

-- 2. AGGREGATION is a GROUP BY on the id — no geometry, trivially distributable.
CREATE TABLE demand_res9 AS
SELECT cell, count(*) AS rides
FROM pings
GROUP BY cell;

-- 3. MULTI-ZOOM: roll res 9 up to city-level res 6 by re-grouping on the parent.
CREATE TABLE demand_res6 AS
SELECT h3_cell_to_parent(cell, 6) AS cell, sum(rides) AS rides
FROM demand_res9
GROUP BY h3_cell_to_parent(cell, 6);

-- 4. JOIN-BY-CELL replaces a spatial join: weather is indexed to the SAME grid.
--    weather_res9(cell, temp_c) built with the same latLngToCell(res 9).
SELECT d.cell, d.rides, w.temp_c
FROM demand_res9 d
JOIN weather_res9 w USING (cell);        -- a STRING-key join, not a geometry join

Step-by-step trace.

Stage	Operation	Why it scales
Index	latlng_to_cell(res 9) per ping	pure function, done once at ingest
Aggregate	GROUP BY cell	string-key group-by, no geometry
Multi-zoom	cell_to_parent(6) re-group	no re-scan of raw points
Cross-dataset join	join on cell	equality join replaces spatial join
Render	cell_to_boundary on hot cells	geometry only at draw time
Grain choice	res 9 fine, roll up on demand	one index, many zoom levels

After the design, every ping is stamped with a res-9 cell id at ingest, so the year-long heatmap is a GROUP BY cell that any distributed engine runs in parallel with no spatial index or polygon test; city-level views come from re-grouping on cell_to_parent(6) rather than re-reading tens of billions of points; and the weather join collapses from a geometry join to a USING (cell) string-equality join because both datasets share the same grid. Geometry appears only when the hot cells are drawn.

Output:

Metric	PostGIS polygon join	H3 grid
Aggregation	billions of ST_Contains tests	one GROUP BY cell
Zoom levels	re-run per grain	re-group on parent
Cross-dataset join	geometry join	string-key equality join
Parallelism	index-bound, node-bound	trivially distributable
Geometry at query time	always	never (only at render)

Why this works — concept by concept:

• latLngToCell indexing — mapping each point to a cell id once at ingest moves all the spatial work to write time, so every read is a plain key aggregation with no geometry or index on the hot path.
• Aggregation as group-by — because the id encodes location, "how many near here" is a GROUP BY cell, which distributes across a cluster exactly like any string-key aggregation, unlike a polygon join that is bound by index and node.
• cellToParent rollup — coarser zoom levels are produced by re-grouping fine cells on their parent id, so one fine index serves every resolution without ever re-scanning the raw points.
• Join-by-cell — indexing two datasets to the same grid turns a cross-dataset spatial join into a string-equality join, the cheapest join a query engine can run, replacing per-pair geometry math entirely.
• Cost — write-time indexing plus read-time key aggregation and equality joins, versus billions of exact containment tests per query. The eliminated cost is the polygon join itself — O(n) group-by on a distributable key instead of O(n·m) geometry tests, at the price of a grid approximating exact boundaries.

Data transformation
Topic — data-transformation
Data transformation problems on gridding and spatial aggregation

Practice →

Data processing Topic — data-processing Data processing problems on group-by and large-scale aggregation

Practice →

---

4. GeoParquet — columnar spatial storage

Store geometry columnar with the CRS in the file, and every engine can read it

The mental model in one line: GeoParquet is ordinary Parquet with a standardised geo metadata block — the geometry column is stored as WKB, the coordinate reference system travels in the file's key-value metadata, and each row group carries a bounding box — so a single columnar file is self-describing, readable by GeoPandas, DuckDB, **Apache Sedona, and GDAL alike, supports predicate pushdown that skips row groups whose bbox misses the query, and finally gives the data lake a spatial format that scales past the row limits and type poverty of shapefile and GeoJSON.** Cluster the rows by space (a Hilbert curve) so the per-row-group bboxes are tight, and a spatial filter reads only the relevant fragments instead of the whole dataset.

What GeoParquet adds to Parquet.

• A geometry column, encoded. Geometry is stored as WKB (well-known binary) in a normal Parquet column, so existing Parquet readers see the bytes and geo-aware readers decode shapes.
• CRS in the metadata. The file's geo metadata records the geometry column(s), the encoding, and the CRS (e.g. PROJJSON for EPSG:4326) — so a coordinate's meaning travels with the data, the thing shapefile sidecars and GeoJSON leave ambiguous.
• Per-row-group bounding boxes. Like any Parquet column, geometry can carry per-row-group statistics; GeoParquet writers record a bbox so readers can skip groups — the basis of spatial predicate pushdown.
• Columnar and cloud-native. Column pruning, compression, and range reads over object storage come for free from Parquet — a spatial dataset now behaves like the rest of the lakehouse.

Why not shapefile or GeoJSON.

• Shapefile. Multi-file, ~2 GB limit, 10-character column names, no CRS certainty without a sidecar .prj — a 1990s interchange format, not lake storage.
• GeoJSON. Human-readable but row-oriented text: huge, slow to parse, no columnar pushdown, and it inflates numeric precision. Fine for a handful of features, unusable at billions of rows.
• GeoParquet. Columnar, typed, compressed, self-describing with CRS, and splittable — the format an engine can push predicates into and read in parallel.

Spatial clustering for pushdown.

• The problem. If rows are in random/arrival order, every row group's bbox spans the whole dataset, so a spatial filter cannot skip anything — pushdown does nothing.
• The fix — space-filling curve order. Sort rows by a Hilbert (or Z-order/geohash) curve so spatially close geometries land in the same row group; each group's bbox becomes tight and local.
• The payoff. A "features in this area" query reads only the row groups whose bbox intersects the area — often a few percent of the file — the columnar analogue of a spatial index.

The failure modes senior engineers pre-empt.

• No spatial clustering. Writing in arrival order makes bboxes overlap everything, so a spatial read scans the whole file. Mitigation: Hilbert/Z-order the rows before writing; size row groups so bboxes stay tight.
• Missing/incorrect CRS metadata. A file without CRS metadata forces every reader to guess the units — the same degrees-versus-metres bug, now baked into storage. Mitigation: always write the CRS; standardise on 4326 for interchange.
• Mixed geometry encodings / versions. Readers differ across GeoParquet versions and encodings (WKB vs native/GeoArrow). Mitigation: pin a version, prefer WKB for portability, validate with a second engine.

Common interview probes on GeoParquet.

• "Why GeoParquet over shapefile/GeoJSON?" — columnar, typed, compressed, CRS-in-file, splittable, pushdown-capable, no size limits.
• "How does spatial pushdown work?" — per-row-group bbox stats; the reader skips groups whose bbox misses the query, so you cluster rows by space to make bboxes tight.
• "Where does the CRS live?" — in the file's geo metadata, so meaning travels with the data.
• "Who can read it?" — GeoPandas, DuckDB spatial, Sedona, GDAL/OGR — one file, many engines.

Worked example — write GeoParquet from GeoPandas, read it in DuckDB

Detailed explanation. The everyday interoperability loop: write a GeoDataFrame to GeoParquet (CRS and geometry preserved) and read it straight back in a completely different engine, DuckDB's spatial extension, with a spatial filter. Publish zones and query them from SQL.

• Write. GeoDataFrame.to_parquet(...) emits GeoParquet with CRS metadata.
• Read. DuckDB spatial reads the same file and runs ST_ predicates.
• The point. No conversion, no export — one file, two engines.

Question. Write a GeoDataFrame of zones to GeoParquet, then read it in DuckDB and run a spatial filter — showing the file is engine-neutral.

Input.

Step	Tool	Operation
Write	GeoPandas	to_parquet (CRS preserved)
Store	object storage	one .parquet file
Read	DuckDB spatial	read_parquet + ST_ filter
Result	—	same geometry, different engine

Code.

# WRITE — GeoPandas emits GeoParquet with the CRS baked into the file metadata.
import geopandas as gpd
zones = gpd.read_file("zones.geojson").set_crs("EPSG:4326")
zones.to_parquet("s3://lake/zones.parquet")     # WKB geometry + geo metadata (CRS)

-- READ — DuckDB's spatial extension reads the SAME file, no conversion step.
INSTALL spatial; LOAD spatial;

-- The geometry column and its CRS come from the file; run a spatial predicate.
SELECT zone_id, name
FROM read_parquet('s3://lake/zones.parquet')
WHERE ST_Contains(
        geom,
        ST_Point(-0.1278, 51.5074)   -- a point; DuckDB decodes the WKB geometry
      );

Step-by-step explanation.

1. set_crs("EPSG:4326") records the CRS on the GeoDataFrame, and to_parquet writes it into the file's geo metadata alongside the WKB-encoded geometry column — so the file is self-describing, not dependent on a sidecar.
2. The output is a single columnar Parquet file in object storage: compressed, splittable, and indistinguishable from any other lake table except for the extra geo metadata.
3. DuckDB's spatial extension reads that file directly with read_parquet; it decodes the WKB geometry and honours the CRS from the metadata — no shapefile import, no GeoJSON parse, no manual reprojection.
4. ST_Contains(geom, ST_Point(...)) runs a real spatial predicate in DuckDB against the geometry that GeoPandas wrote — the same bytes, interpreted identically, because both engines follow the GeoParquet spec.
5. This is the interoperability payoff: the producer (Python/GeoPandas) and the consumer (SQL/DuckDB) share one file with no conversion, which is exactly what a shapefile or GeoJSON hand-off cannot do cleanly at scale.

Output.

zone_id	name
12	Downtown

Rule of thumb. Write geometry once as GeoParquet with its CRS, and every engine — GeoPandas, DuckDB, Sedona, GDAL — reads the same file natively with no conversion. The CRS lives in the file, so meaning travels with the data instead of in a fragile sidecar or a reader's assumption.

Worked example — Hilbert-order clustering and bbox pushdown

Detailed explanation. Pushdown only helps if the row-group bboxes are tight, and that requires clustering rows by space before writing. Sort by a Hilbert curve so spatially near rows share a row group, then a spatial filter skips most of the file. Compare a randomly ordered write to a Hilbert-ordered one.

• The lever. Order rows by a space-filling curve (Hilbert) on the geometry's centroid/cell.
• The effect. Each row group's bbox becomes small and local instead of spanning everything.
• The payoff. A bounded-area query reads only intersecting row groups.

Question. Show why Hilbert-ordering geometries before writing GeoParquet enables bbox pushdown, and estimate the read reduction versus random order.

Input.

Layout	Row-group bbox	A local query reads
Random / arrival order	spans whole dataset	all row groups
Hilbert / Z-order	tight, local	only intersecting groups
Effect	pushdown useless	pushdown skips most

Code.

import geopandas as gpd
from shapely import Point

gdf = gpd.read_parquet("s3://lake/pings.parquet")   # billions of point rows

# 1. Compute a Hilbert-curve key so spatially close rows sort together.
#    (GeoPandas exposes hilbert_distance on the geometry's bounds/centroid.)
gdf["hilbert"] = gdf.geometry.hilbert_distance()     # 1-D key preserving 2-D locality
gdf = gdf.sort_values("hilbert")

# 2. Write with bounded row-group size so each group covers a small AREA.
gdf.to_parquet(
    "s3://lake/pings_clustered.parquet",
    row_group_size=100_000,           # smaller groups -> tighter per-group bbox
    write_covering_bbox=True,         # store per-row-group bbox for pushdown
)

-- 3. A bounded-area query now SKIPS row groups whose bbox misses the window.
--    DuckDB / Sedona use the per-row-group bbox to prune before reading data.
SELECT count(*)
FROM read_parquet('s3://lake/pings_clustered.parquet')
WHERE ST_Within(geom, ST_MakeEnvelope(-0.20, 51.45, -0.05, 51.55));  -- a map viewport
-- Reads only the row groups intersecting the envelope, not the whole file.

Step-by-step explanation.

1. hilbert_distance() maps each 2-D geometry to a 1-D key along a Hilbert space-filling curve, which has the property that points close on the curve are close in space — so sorting by it clusters spatially near rows together.
2. Sorting by hilbert before writing means consecutive rows (and therefore each row group) cover a small, contiguous patch of space rather than a random scatter across the whole extent.
3. row_group_size=100_000 with write_covering_bbox=True produces many row groups, each carrying a tight bounding box in its metadata — tight because the rows within it are spatially local.
4. The ST_Within(geom, envelope) query lets the reader compare the query envelope to each row group's bbox and skip every group that does not intersect — often 90%+ of the file for a small viewport — reading only the relevant fragments.
5. Without the Hilbert sort, every row group's bbox would span nearly the whole dataset, every bbox would intersect the query, and the reader would fall back to scanning everything — pushdown present but useless. Clustering is what makes the bbox stats selective.

Output.

Layout	Row groups read (small viewport)	I/O
Random order	~all	full scan
Hilbert order	few intersecting	a few % of file
Effect	no skipping	heavy skipping
Analogue	table scan	spatial index

Rule of thumb. Cluster rows by a space-filling curve (Hilbert/Z-order) before writing GeoParquet and keep row groups modest, so each group's bbox is tight and a bounded-area query skips the groups it does not need. Pushdown is only as good as your spatial clustering — random order makes the bbox statistics worthless.

Worked example — one file, three engines: the interoperability win

Detailed explanation. The strategic value of GeoParquet is that the same file feeds an ad-hoc Python analyst, a SQL engine, and a distributed Spark job with no format conversion — each reads geometry and CRS natively. Show one file consumed by GeoPandas, DuckDB, and Sedona.

• Producer. One pings.parquet written once.
• Consumers. GeoPandas (exploration), DuckDB (SQL/BI), Sedona (distributed join).
• The invariant. All three decode the same WKB geometry and read the same CRS.

Question. Read a single GeoParquet file from three engines and note that geometry and CRS are consistent across all of them, with no conversion.

Input.

Engine	Role	Read call
GeoPandas	interactive analysis	read_parquet
DuckDB spatial	SQL / BI	read_parquet + ST_
Apache Sedona	distributed compute	format("geoparquet")

Code.

# Consumer A — GeoPandas: interactive exploration, geometry + CRS restored.
import geopandas as gpd
g = gpd.read_parquet("s3://lake/pings.parquet")
print(g.crs, g.geometry.geom_type.value_counts())     # EPSG:4326, Point

-- Consumer B — DuckDB: SQL/BI over the SAME file, no import.
LOAD spatial;
SELECT count(*) FROM read_parquet('s3://lake/pings.parquet')
WHERE ST_Within(geom, ST_MakeEnvelope(-0.2, 51.4, 0.0, 51.6));

# Consumer C — Apache Sedona: distributed read of the SAME file for a big join.
df = sedona.read.format("geoparquet").load("s3://lake/pings.parquet")
df.createOrReplaceTempView("pings")   # now joinable at cluster scale

Step-by-step explanation.

1. GeoPandas reads the file and reconstructs the CRS (EPSG:4326) and geometry types from the geo metadata, so an analyst starts exploring immediately with no reprojection or type inference.
2. DuckDB reads the identical file with its spatial extension and runs SQL ST_ predicates against the same geometry — the BI/SQL persona gets the data with zero hand-off from the Python producer.
3. Sedona loads the same path with format("geoparquet") into a distributed DataFrame, ready for a cluster-scale spatial join — the heavy-compute persona, again from the one file.
4. Because all three follow the GeoParquet spec, the WKB geometry decodes identically and the CRS is read (not guessed) everywhere, so results are consistent across engines — the property that makes GeoParquet a genuine interchange format, not just "Parquet with a geometry blob."
5. Contrast the alternative: a shapefile would need conversion and column-name mangling, a GeoJSON would be reparsed and reprecisioned per engine, and CRS would depend on sidecars or assumptions — GeoParquet removes all of that friction with one self-describing file.

Output.

Engine	Reads geometry	Reads CRS	Conversion needed
GeoPandas	yes (WKB → shapely)	yes (from metadata)	none
DuckDB	yes (ST_ on WKB)	yes	none
Apache Sedona	yes (distributed)	yes	none
net	consistent everywhere	consistent	zero

Rule of thumb. Treat GeoParquet as the interchange contract of the geospatial lakehouse: write once, and GeoPandas, DuckDB, and Sedona all read the same geometry and CRS with no conversion. One self-describing file replaces per-engine exports, format conversions, and CRS guesswork.

Senior interview question on designing geospatial lake storage

A senior interviewer might ask: "Design the storage layer for a geospatial lakehouse holding billions of point observations and millions of polygons, consumed by a Python analytics team, a DuckDB-backed BI layer, and nightly Spark jobs. Choose the format and justify it over shapefile/GeoJSON, explain how you lay the data out so spatial filters do not scan everything, how the CRS is guaranteed across engines, and how the same files feed all three consumers."

Solution Using GeoParquet, Hilbert clustering with bbox pushdown, CRS-in-file, and multi-engine reads

# 1. FORMAT: GeoParquet — columnar, typed, compressed, CRS-in-file, splittable.
#    Points partitioned by date; polygons a single clustered file.
import geopandas as gpd
zones.set_crs("EPSG:4326").to_parquet("s3://lake/zones.parquet")

# 2. LAYOUT: Hilbert-cluster points so per-row-group bboxes are tight (pushdown).
pings["hilbert"] = pings.geometry.hilbert_distance()
(pings.sort_values("hilbert")
      .to_parquet("s3://lake/pings/",
                  partition_cols=["date"],       # partition prune by day
                  row_group_size=100_000,         # + bbox prune within a day
                  write_covering_bbox=True))

-- 3. PUSHDOWN: a viewport query prunes by partition (date) AND row-group bbox.
SELECT count(*) FROM read_parquet('s3://lake/pings/date=2026-08-26/*.parquet')
WHERE ST_Within(geom, ST_MakeEnvelope(-0.2, 51.45, -0.05, 51.55));

# 4. MULTI-ENGINE: the same files feed Python, SQL/BI, and distributed Spark.
g   = gpd.read_parquet("s3://lake/zones.parquet")                       # analysts
# DuckDB: SELECT ... FROM read_parquet('s3://lake/pings/...')           # BI
sdf = sedona.read.format("geoparquet").load("s3://lake/pings/")        # nightly jobs

Step-by-step trace.

Decision	Choice	Why
Format	GeoParquet	columnar, typed, CRS-in-file, no size limit
Partitioning	by date	prune whole days before any I/O
Clustering	Hilbert order	tight row-group bboxes → spatial pushdown
Pushdown	partition + bbox	viewport reads a few % of a day
CRS	in file metadata	every engine reads, none guesses
Consumers	GeoPandas / DuckDB / Sedona	one file set, three personas

After the design, geometry is stored as GeoParquet with the CRS in every file, so shapefile size limits and GeoJSON parsing costs are gone and no reader has to guess units; points are partitioned by date and Hilbert-clustered within each day, so a viewport query prunes to one day and then to the handful of row groups whose bbox intersects the window; and the identical files feed the Python, DuckDB, and Sedona consumers with no conversion. Storage is written once and read by everyone.

Output:

Metric	Shapefile / GeoJSON	GeoParquet lakehouse
Size / row limits	~2 GB / bloated text	none (columnar, compressed)
CRS guarantee	sidecar / assumed	in file, read by all engines
Spatial filter	full scan	partition + bbox pushdown
Multi-engine read	convert per engine	one file, native everywhere
Analyst / BI / Spark	separate exports	shared files

Why this works — concept by concept:

• GeoParquet as the format — columnar, compressed, typed storage with the CRS in the file removes shapefile's limits and GeoJSON's parse cost while making geometry a first-class lake column engines can push predicates into.
• Partition then cluster — partitioning by date prunes whole days before any read, and Hilbert-clustering within a day makes row-group bboxes tight, so a viewport query is pruned twice and reads a tiny fraction of the data.
• bbox pushdown — per-row-group bounding-box statistics let the reader skip non-intersecting groups, giving the file the selectivity of a spatial index without a separate index structure.
• CRS in the file — writing the coordinate reference system into the metadata means every engine reads the true units instead of assuming, killing the degrees-versus-metres bug at the storage layer.
• Cost — write-once columnar files with partition and bbox pruning read by every engine, versus per-engine exports and full scans of row-oriented formats. The eliminated cost is repeated conversion and whole-dataset scans — O(selected fragments) reads instead of O(dataset), shared across every consumer.

Data processing
Topic — data-processing
Data processing problems on columnar formats and predicate pushdown

Practice →

Optimization Topic — optimization Optimization problems on file layout, clustering, and I/O pruning

Practice →

---

5. Apache Sedona — distributed spatial joins on Spark

Co-partition by space, index each partition, then join — or watch it become a cartesian blow-up

The mental model in one line: Apache Sedona extends Spark (and Flink) with spatial types, ST_ functions, and — the part that actually matters — spatial partitioning, so a distributed **spatial join works by co-partitioning both datasets into the same spatial grid (a quad-tree or KDB-tree), building a local spatial index inside each partition, and joining only within co-located partitions, which is the only way a join of billions of points against millions of polygons avoids the O(n·m) cartesian explosion a naive crossJoin would produce.** Read GeoParquet in, use ST_Contains/ST_Intersects for range joins and distance functions for kNN, broadcast the small side when one exists, and the same geometry work that overwhelms a single Postgres box finishes across a cluster.

What Sedona is.

• Spatial types and functions on Spark. Sedona adds geometry columns and the ST_ function family to Spark SQL and the DataFrame API, so spatial logic runs inside Spark's distributed execution engine.
• Spatial partitioning and indexing. The differentiator: Sedona can repartition data into a spatial grid (KDB-tree/quad-tree) and build an in-partition R-tree, so joins are local, not global.
• Reads the lake, speaks the grid. Sedona reads GeoParquet natively and interoperates with H3, so it slots directly onto the storage and gridding from the previous sections.

The distributed spatial-join problem.

• Why a naive join explodes. A spatial predicate like ST_Contains(poly, pt) is not an equality, so Spark cannot hash-join on it; without spatial partitioning it degrades to a crossJoin — every point compared to every polygon, O(n·m).
• The three-step fix. (1) Partition both sides into the same spatial grid so a point and the polygons that might contain it land on the same node; (2) index each partition with a local R-tree; (3) join only within co-located partitions, running filter-and-refine locally.
• Broadcast when one side is small. If the polygons fit in memory, broadcast them to every partition and skip repartitioning the huge side — the spatial analogue of a broadcast hash join.

Range joins and kNN joins.

• Range join. ST_Contains, ST_Intersects, ST_DWithin — "which points fall in which polygons," the point-in-polygon enrichment at cluster scale.
• kNN / distance join. Nearest-neighbour at scale (e.g. ST_DistanceSphere / ST_Distance with a k-nearest join), for "nearest facility to each event."
• Predicate drives the strategy. Sedona picks a spatial-partitioned join or a broadcast join based on the predicate and sizes; you help it by broadcasting a small side and reading pre-clustered GeoParquet.

The failure modes senior engineers pre-empt.

• No spatial partitioning → cartesian blow-up. The default failure: a spatial join with no co-partitioning is O(n·m) and never finishes. Mitigation: let Sedona spatially partition, or broadcast the small side; confirm the plan is not a crossJoin.
• Partition skew. Dense areas (a city centre) put most geometries in one partition, so one task runs forever. Mitigation: a spatial partitioner that balances by density (KDB-tree adapts to data), and adequate partition count.
• Forgetting to broadcast a small side. Repartitioning a billion-row side to join a tiny polygon set wastes a shuffle. Mitigation: broadcast the small side and leave the big side in place.

Common interview probes on Sedona.

• "Why can't Spark hash-join a spatial predicate?" — it is not equality; without spatial partitioning it becomes a cartesian cross join.
• "How does Sedona make a spatial join scale?" — co-partition both sides into a spatial grid, build a local index per partition, join within co-located partitions.
• "When do you broadcast?" — when one side (usually polygons) is small enough to fit in memory; broadcast it and skip repartitioning the big side.
• "What causes a slow Sedona job?" — no spatial partition (cross join) or skew concentrating geometries in one partition.

Worked example — a distributed range join (points in polygons)

Detailed explanation. The workhorse Sedona job: a point-in-polygon range join across a cluster. Read both sides from GeoParquet, let Sedona spatially partition, and join with ST_Contains. Enrich billions of pings with their zone.

• Read. Both pings and zones from GeoParquet (geometry + CRS preserved).
• Join. ST_Contains(zone, ping) — a range join Sedona partitions spatially.
• Result. Each ping tagged with its zone, computed distributed.

Question. Write a Sedona SQL range join that assigns each ping its containing zone across a cluster, and explain how it avoids a cartesian product.

Input.

Piece	Value
Left	pings (billions, points)
Right	zones (millions, polygons)
Predicate	ST_Contains(z.geom, p.geom)
Strategy	spatial partition + local index

Code.

from sedona.spark import SedonaContext
sedona = SedonaContext.create(spark)

# 1. Read both sides from GeoParquet — geometry + CRS come from the files.
sedona.read.format("geoparquet").load("s3://lake/pings/date=2026-08-26/") \
      .createOrReplaceTempView("pings")
sedona.read.format("geoparquet").load("s3://lake/zones.parquet") \
      .createOrReplaceTempView("zones")

# 2. Range join: Sedona spatially co-partitions both sides and indexes each
#    partition, so ST_Contains runs LOCALLY (filter-and-refine), not as a cross join.
enriched = sedona.sql("""
  SELECT p.ping_id, z.zone_id
  FROM pings p
  JOIN zones z
    ON ST_Contains(z.geom, p.geom)     -- range predicate: point-in-polygon
""")
enriched.write.format("geoparquet").save("s3://lake/pings_enriched/date=2026-08-26/")

Step-by-step explanation.

1. Reading GeoParquet gives Sedona typed geometry columns with the CRS already set, so no parsing or reprojection precedes the join — the storage layer feeds the compute layer directly.
2. When Sedona sees a spatial predicate (ST_Contains), it plans a spatial join: it repartitions both pings and zones into the same spatial grid so a ping and the zones that could contain it are co-located on one executor.
3. Inside each partition Sedona builds a local R-tree over the polygons, so the point-in-polygon test is filter-and-refine locally — bounding-box filter then exact ST_Contains — exactly the two-phase model from section 1, now per partition.
4. Because the join only compares co-located geometries, total work is roughly O(n log m) across the cluster instead of the O(n·m) a crossJoin would do — the spatial partitioning is what turns an impossible join into a linear-ish distributed one.
5. Writing the result back to partitioned GeoParquet keeps the enriched output in the same lake format, ready for the next engine — the pipeline stays columnar and engine-neutral end to end.

Output.

ping_id	zone_id
900001	12
900002	7
900003	12
900004	3

Rule of thumb. For a distributed point-in-polygon join, read GeoParquet and let Sedona spatially co-partition both sides and index each partition — the ST_Contains then runs filter-and-refine locally instead of as a cross join. The spatial partition is the difference between an O(n log m) distributed join and an O(n·m) explosion.

Worked example — a broadcast spatial join for a small polygon side

Detailed explanation. When one side is small — a few thousand zones against billions of pings — repartitioning the huge side is wasteful. Broadcast the small polygon side to every executor and join in place, the spatial analogue of a broadcast hash join. Broadcast the zones.

• The condition. The polygon side fits in executor memory.
• The move. BROADCAST(z) hint (or broadcast()), so zones are copied to every task.
• The win. No shuffle of the billion-row side; each partition joins locally against the in-memory zones.

Question. Rewrite the point-in-polygon join to broadcast the small zone side, and explain when this beats spatial repartitioning.

Input.

Aspect	Spatial-partition join	Broadcast join
Small side fits in memory	not required	required
Shuffle of big side	yes (repartition)	no
Best when	both sides large	one side small
Mechanism	co-partition grid	copy small side to all tasks

Code.

# Small side: a few thousand zone polygons → broadcast to every executor,
# so the billion-row pings side is NEVER shuffled/repartitioned.
enriched = sedona.sql("""
  SELECT /*+ BROADCAST(z) */ p.ping_id, z.zone_id
  FROM pings p
  JOIN zones z
    ON ST_Contains(z.geom, p.geom)
""")

# Equivalent DataFrame form with an explicit broadcast + local spatial index.
from pyspark.sql.functions import broadcast
enriched_df = pings.join(
    broadcast(zones),                       # zones copied to each task
    ST_Contains(zones.geom, pings.geom)     # local point-in-polygon per partition
)

Step-by-step explanation.

1. The /*+ BROADCAST(z) */ hint tells Sedona/Spark to ship the entire zones dataset to every executor, so each task holds a full local copy of the polygons in memory.
2. Because every task already has all the polygons, the huge pings side stays exactly where it is — no repartition, no shuffle of billions of rows, which is usually the most expensive part of a distributed join.
3. Each task builds a local R-tree over the broadcast polygons and runs filter-and-refine ST_Contains for its slice of pings — the same two-phase test, but with zero data movement on the big side.
4. This beats spatial repartitioning precisely when the small side fits in memory: you trade a small broadcast (thousands of polygons) for eliminating a massive shuffle (billions of points), which is almost always the right trade when it is available.
5. The guardrail is memory: if the "small" side is not actually small, the broadcast blows up executor memory and you fall back to spatial-partitioned join — so the senior move is to check the small side's size before broadcasting.

Output.

Strategy	Shuffle of pings	When to use
Broadcast zones	none	zones fit in memory
Spatial repartition	yes	both sides large
Naive cross join	—	never (O(n·m))
chosen here	broadcast	small polygon side

Rule of thumb. When the polygon side is small enough to fit in memory, broadcast it and join in place — you trade a tiny broadcast for eliminating a shuffle of the billion-row side. Reserve spatial repartitioning for when both sides are large; never let either fall back to a cross join.

Worked example — a distributed kNN / distance join

Detailed explanation. "Nearest facility to every event" is a kNN join, and at cluster scale it needs the same spatial-partition-then-index approach, plus a distance function. Find, for each ping, the nearest of a set of facilities using a Sedona distance join.

• The question. For each point on the left, the nearest (or k nearest) on the right.
• The mechanism. Spatial partitioning co-locates candidates; a local index ranks by distance.
• The function. ST_DistanceSphere/ST_Distance measures; ST_DWithin can bound the search.

Question. Write a Sedona distance join that finds facilities within 2 km of each ping (a bounded nearest-neighbourhood), and note how it scales.

Input.

Piece	Value
Left	pings (points)
Right	facilities (points)
Predicate	ST_DWithin/distance within 2 km
Strategy	spatial partition + distance rank

Code.

# Bounded distance join: every ping paired with facilities within 2 km,
# ordered by distance — the scalable form of "nearest facility to each event".
nearby = sedona.sql("""
  SELECT p.ping_id, f.facility_id,
         ST_DistanceSphere(p.geom, f.geom) AS metres      -- true metres on the sphere
  FROM pings p
  JOIN facilities f
    ON ST_DWithin(p.geom, f.geom, 2000)                   -- <= 2 km, spatially joined
""")

# Reduce to the single nearest facility per ping (k = 1) with a window.
from pyspark.sql import Window
from pyspark.sql.functions import row_number
w = Window.partitionBy("ping_id").orderBy("metres")
nearest = (nearby.withColumn("rk", row_number().over(w))
                 .filter("rk = 1").drop("rk"))

Step-by-step explanation.

1. ST_DWithin(p.geom, f.geom, 2000) is a bounded distance predicate — it restricts the join to pairs within 2 km — which lets Sedona spatially partition on that radius and compare only co-located candidates instead of all pairs.
2. Sedona co-partitions pings and facilities into a spatial grid sized to the search radius, so each ping only meets the facilities that could plausibly be within 2 km, keeping the join local and bounded.
3. ST_DistanceSphere computes the true great-circle distance in metres for the surviving candidate pairs, so the reported distances are correct without a projection step.
4. The row_number() window partitioned by ping_id ordered by distance reduces the bounded neighbourhood to the single nearest facility (k=1); widening to rk <= k yields k-nearest — the distance join plus a window is the scalable kNN pattern.
5. This scales because the expensive part — comparing every ping to every facility — is pruned by the spatial partition and the ST_DWithin bound; only the local, bounded candidate sets reach the distance computation and ranking.

Output.

ping_id	facility_id	metres
900001	51	412
900002	88	1,203
900003	51	980
900004	12	1,740

Rule of thumb. Express nearest-neighbour at scale as a bounded distance join (ST_DWithin on a radius) plus a per-key window ordered by ST_DistanceSphere, so Sedona spatially partitions on the radius and only ranks local candidates. An unbounded all-pairs distance computation does not scale; the radius bound and spatial partition are what make kNN distributable.

Senior interview question on distributed spatial joins at scale

A senior interviewer might ask: "You must join five billion GPS pings against three million zone polygons every night to enrich the pings with their zone, and a single-node PostGIS job cannot finish in the window. Design it in Apache Sedona: how you read the data, how you make the spatial join scale instead of exploding into a cartesian product, when you broadcast versus spatially partition, how you handle skew from dense cities, and how the output stays in the lakehouse."

Solution Using GeoParquet reads, spatial partitioning, a broadcast small side, skew handling, and a GeoParquet write

# 1. READ both sides from GeoParquet (geometry + CRS from the files, columnar, pruned).
from sedona.spark import SedonaContext
sedona = SedonaContext.create(spark)
pings = sedona.read.format("geoparquet").load("s3://lake/pings/date=2026-08-26/")
zones = sedona.read.format("geoparquet").load("s3://lake/zones.parquet")
pings.createOrReplaceTempView("pings"); zones.createOrReplaceTempView("zones")

# 2. JOIN that SCALES: broadcast the small (3M) polygon side; the 5B pings never shuffle.
#    Sedona builds a local R-tree over the broadcast polygons → filter-and-refine per task.
enriched = sedona.sql("""
  SELECT /*+ BROADCAST(z) */ p.ping_id, z.zone_id
  FROM pings p JOIN zones z
    ON ST_Contains(z.geom, p.geom)      -- range join, NOT a crossJoin
""")

# 3. If polygons are too big to broadcast, spatially partition BOTH sides instead.
#    (KDB-tree adapts to density → mitigates city-centre skew by balancing partitions.)
sedona.conf.set("sedona.join.gridtype", "kdbtree")
sedona.conf.set("sedona.join.numpartition", "2048")   # enough partitions to spread skew

# 4. WRITE back to partitioned GeoParquet — output stays in the lakehouse for all engines.
enriched.write.format("geoparquet").mode("overwrite") \
        .save("s3://lake/pings_enriched/date=2026-08-26/")

# 5. VERIFY the plan is a spatial/broadcast join, not a cartesian product.
enriched.explain()   # expect BroadcastIndexJoin / spatial join — NOT CartesianProduct

Step-by-step trace.

Concern	Mechanism	Effect
Read	GeoParquet	typed geometry + CRS, columnar, pruned
Scale the join	broadcast small side	5B pings never shuffled
Local test	per-task R-tree	filter-and-refine per partition
Fallback	KDB-tree spatial partition	when polygons too big to broadcast
Skew	more partitions + adaptive grid	dense cities spread across tasks
Verify	explain()	confirms spatial join, not cartesian

After the design, both sides are read from GeoParquet with geometry and CRS intact; because the three-million polygon side fits in memory it is broadcast, so the five-billion-row ping side is never shuffled and each task runs a local filter-and-refine ST_Contains against an in-memory R-tree; if the polygons were too large, Sedona instead co-partitions both sides with a density-adaptive KDB-tree and enough partitions to spread city-centre skew; the enriched output is written back to partitioned GeoParquet; and explain() confirms a spatial/broadcast join rather than a CartesianProduct — the check that the join actually scales.

Output:

Metric	Single-node PostGIS	Sedona distributed
5B×3M nightly join	misses the window	finishes across the cluster
Join complexity	O(n·m) if unpartitioned	O(n log m), co-partitioned
Big-side shuffle	n/a	none (broadcast small side)
Skew	one hot node	spread (adaptive partitions)
Output	DB table	partitioned GeoParquet (lakehouse)

Why this works — concept by concept:

• GeoParquet reads — reading typed geometry with the CRS already set feeds Sedona directly from the lake with columnar pruning, so the compute layer starts from clean spatial data with no parse or reprojection step.
• Broadcast the small side — copying the three-million polygon side to every executor eliminates the shuffle of the five-billion-row side, turning the join into local filter-and-refine per task — the cheapest distributed spatial join when one side fits in memory.
• Spatial partitioning fallback — when neither side is broadcastable, co-partitioning both into a density-adaptive KDB-tree grid co-locates candidate geometries so the join is local and O(n log m) instead of the O(n·m) cartesian a naive join produces.
• Skew handling — an adaptive grid plus enough partitions spreads dense areas across tasks, so one busy city centre does not become a single straggler that defines the job's runtime.
• Cost — a broadcast plus local indexed joins (or a balanced spatial repartition), versus an impossible single-node scan or a cluster-wide cartesian product. The eliminated cost is the O(n·m) explosion and the single-machine ceiling — O(n log m) distributed filter-and-refine, verified by the query plan.

Joins
Topic — joins
Join problems on distributed and broadcast spatial joins

Practice →

Optimization
Topic — optimization
Optimization problems on partitioning, skew, and shuffle avoidance

Practice →

---

Cheat sheet — geospatial data engineering

• A location is a geometry, not two numbers. Model it as Point/LineString/Polygon with an explicit CRS/SRID. Store in 4326 (degrees) for interchange; compute distance/area in geography (ellipsoidal metres) or a projected metric CRS — never planar math on 4326 degrees. The degrees-versus-metres bug is the field's most common.
• A spatial index is a bounding-box index. A B-tree indexes a total order and is useless for "overlaps this box." Spatial engines use an R-tree/GiST that runs filter-and-refine: phase 1 filters by bounding box (&&), phase 2 refines with the exact ST_ predicate on candidates. A slow spatial query almost always means no spatial index.
• PostGIS = serving. geometry(_, 4326) + CREATE INDEX ... USING GIST (geom); ST_Contains for point-in-polygon; ST_DWithin(a::geography, b::geography, metres) for indexed metre-accurate radius (never ST_Distance < d); geom <-> point ORDER BY ... LIMIT k for indexed kNN. Verify with EXPLAIN — Index Scan, not Seq Scan.
• H3 = aggregation. latlng_to_cell(lat, lng, res) at ingest turns "how many near here" into GROUP BY cell; grid_disk(cell, k) (k-ring) for neighbourhoods; cell_to_parent(cell, res) to roll up zoom levels without re-scanning; compact_cells to shrink coverage sets. Hexagons beat a geohash: equidistant neighbours, no latitude distortion. A grid approximates — use PostGIS for exact boundaries.
• Resolution picks the grain. City ≈ res 6–7, neighbourhood ≈ 8–9, building ≈ 11–12. Index fine once; roll up with cell_to_parent. Too coarse hides signal; too fine explodes cardinality.
• GeoParquet = storage/interchange. Columnar Parquet + geometry as WKB + CRS in file metadata + per-row-group bbox. Reads in GeoPandas, DuckDB, Sedona, GDAL with no conversion. Beats shapefile (2 GB limit, sidecar CRS) and GeoJSON (row-oriented text, no pushdown).
• Cluster for pushdown. Sort rows by a Hilbert/Z-order curve before writing so row-group bboxes are tight; then a viewport query skips non-intersecting groups (a few % of the file). Random order makes bbox stats useless. Partition by date and cluster by space for double pruning.
• Apache Sedona = distributed joins. A spatial predicate is not equality, so without spatial partitioning it becomes a cartesian cross join. Fix: co-partition both sides into a spatial grid (KDB/quad-tree), build a per-partition local index, join locally. Broadcast the small side (polygons) when it fits in memory so the huge side never shuffles.
• Range vs kNN join. Range: ST_Contains/ST_Intersects/ST_DWithin (point-in-polygon, within-radius). kNN: a bounded distance join (ST_DWithin radius) + a per-key window ordered by ST_DistanceSphere. Never an unbounded all-pairs distance.
• Skew is the distributed enemy. Dense areas concentrate geometries in one partition. Mitigate with a density-adaptive grid (KDB-tree) and enough partitions; verify the plan is a spatial/broadcast join, not a CartesianProduct.
• Tool by job. Serve exact geometry → PostGIS. Aggregate on a grid → H3. Store for many engines → GeoParquet. Distribute a heavy join → Apache Sedona. They are complementary; a real stack runs all four, with GeoParquet feeding Sedona and H3 grids serving heatmaps.

Frequently asked questions

What is geospatial data engineering?

Geospatial data engineering is the discipline of storing, indexing, transforming, and serving location data — coordinates, lines, and polygons — at scale, and it is distinct from ordinary data engineering because a location is a geometry in a specific coordinate reference system, not a scalar you can GROUP BY or index with a B-tree. The three hard problems it exists to solve are representing geometry with its CRS attached, answering relationship questions (contains, intersects, nearest, within a distance) that no equality or range predicate expresses, and indexing shapes with a bounding-box structure (an R-tree) that supports the filter-and-refine query model. In practice it means choosing the right tool for each job — PostGIS to serve exact spatial queries, H3 to aggregate points on a hexagonal grid, GeoParquet to store geometry every engine can read, and Apache Sedona to run spatial joins across a cluster — and getting the CRS and the spatial index right so queries are both correct and fast.

PostGIS vs H3 — when do I use each?

Use PostGIS when you need exact geometry and low-latency serving: point-in-polygon against a real boundary, within-radius search, nearest-neighbour, or anything an application calls per request. It gives you true geometry types, GiST spatial indexes, and the full ST_ function library, so the answer is exact and index-fast. Use H3 when the question is aggregation on a grid: heatmaps, demand surfaces, or joining two point datasets by shared cell. Snapping each point to a hexagon id turns spatial aggregation into a GROUP BY on a string key with no geometry at query time, which parallelises trivially and is far cheaper than repeated polygon joins. The trade-off is precision: an H3 cell only approximates a real boundary, so you would not use it to decide whether a point is inside a legal district. Many stacks use both — H3 for the analytical grid, PostGIS for exact serving — and they interoperate cleanly.

Why can't I just use a B-tree index for location queries?

A B-tree indexes a total order — it is built for =, <, >, and BETWEEN on a single sortable key — and location is inherently two-dimensional with no meaningful total order, so a B-tree cannot answer "which shapes overlap this box" or "what is within five kilometres." Sorting by latitude tells you nothing about longitude, so a range scan on one axis still has to check every row on the other. Spatial engines instead use an R-tree (in PostGIS, a GiST index) that indexes each geometry's bounding box and supports the two-phase filter-and-refine model: phase one uses the index to return every shape whose bounding box overlaps the query box (cheap, approximate), and phase two runs the exact geometric predicate only on those candidates. Without that spatial index, an ST_Contains or ST_DWithin has no phase-one filter and must run the expensive exact test on every row — which is why a missing GiST index is the single most common cause of a slow spatial query.

What is GeoParquet and why not shapefile or GeoJSON?

GeoParquet is ordinary Parquet plus a standardised geo metadata block: the geometry is stored as well-known binary (WKB) in a normal column, the coordinate reference system travels in the file's metadata, and each row group carries a bounding box for predicate pushdown. That makes it columnar, compressed, splittable, self-describing, and readable by GeoPandas, DuckDB, Apache Sedona, and GDAL with no conversion — the geometry a data lake actually needs. Shapefile is a multi-file 1990s format with a ~2 GB limit, 10-character column names, and a CRS that lives in a fragile sidecar; GeoJSON is human-readable but row-oriented text that is huge, slow to parse, and offers no columnar pushdown. Neither scales to billions of rows or supports the predicate pushdown a query engine wants. With GeoParquet you also get real spatial selectivity by clustering rows along a Hilbert curve so per-row-group bounding boxes are tight and a bounded-area query skips most of the file.

When do I need Apache Sedona instead of PostGIS?

Reach for Apache Sedona when the spatial join no longer fits on one machine — the classic case being billions of points joined against millions of polygons in a batch window a single PostGIS box cannot meet. PostGIS is excellent for serving and for joins that fit one node's memory and disk, but a spatial predicate is not an equality, so a huge join cannot be hash-joined and, without spatial partitioning, degrades to a cartesian product. Sedona solves this by distributing the work: it co-partitions both datasets into the same spatial grid (a quad-tree or KDB-tree), builds a local spatial index inside each partition, and joins only co-located geometries, and it can broadcast a small polygon side so the huge side never shuffles. It also reads GeoParquet natively and speaks H3, so it drops onto your lake storage and grid directly. The rule of thumb: PostGIS to serve, Sedona to crunch cluster-scale batch joins.

What CRS or SRID should I store my data in?

Store your data in EPSG:4326 (WGS84 longitude/latitude) as the interchange default — it is what almost every source emits, what GeoParquet and web maps expect, and what keeps data portable across engines. But do not compute distance or area on 4326 as if it were planar, because its units are degrees, not metres: a degree of longitude is ~111 km at the equator and shrinks to zero at the poles. When you need real measurements, either cast to the geography type (which does ellipsoidal math and returns true metres globally with no projection choice) or ST_Transform into a projected metric CRS appropriate to your region — a UTM zone for local accuracy, or Web Mercator (3857) for web-map alignment while accepting its distortion away from the equator. The discipline is simple: store in 4326, always attach the SRID to the column, and switch to a metric CRS or geography the moment a query measures distance or area.

Practice on PipeCode

• Drill the spatial indexing practice library → for the GiST/R-tree, bounding-box, and filter-and-refine problems that make PostGIS and Sedona joins fast instead of full scans.
• Rehearse point-in-polygon and within-radius work on the spatial join practice library → for the range-join, broadcast-join, and kNN patterns that turn a cartesian blow-up into a distributed join.
• Sharpen the architecture axis with the system design practice library → for the serve-vs-aggregate-vs-store-vs-distribute decisions a geospatial stack must get right.
• Stack the prerequisites against PipeCode's broader 450+ data-engineering catalogue to anchor the CRS, spatial-index, H3-grid, and GeoParquet-layout patterns against real graded inputs — geometry, projections, gridding, and distributed joins.

Lock in geospatial data engineering muscle memory

Docs explain PostGIS, H3, GeoParquet, and Sedona. PipeCode drills explain the decision — when a `spatial index` is a bounding box and not a sorted key, when an H3 `GROUP BY cell` beats a polygon join, when to cluster GeoParquet by a Hilbert curve, and when to broadcast the small side of a distributed `spatial join`. Pipecode.ai is Leetcode for Data Engineering — geospatial practice tuned for the production trade-offs senior data engineers actually face.

Practice spatial indexing problems →
Practice spatial join problems →