Geospatial, but make it fast

Pickups, routes, and taxi zones

We can use the operations we introduced above to analyze taxi trips: match each pickup to a taxi zone, then add up fares by zone.

The examples below use pickups, routes, and zones in Lower Manhattan, New York City. A pickup gives us a point, a route gives us a path, and a zone gives us a polygon. The pickup, route, and zone tables each have a geometry column beside columns such as a trip_id or zone (name of the zone).

New York City's Taxi and Limousine Commission (TLC) uses areas called taxi zones to report pickup and drop-off locations in its public trip data. Our examples use nineteen simplified zone boundaries from the official taxi-zone dataset. The illustrated route follows a section of the eastbound M8 bus, from Hudson Street to Avenue A, using MTA route data. The pickup locations are illustrative.^[3]

You don't need to know your way around Manhattan to follow along!

Spatial filter

Let's start with the pickup location for trip 1042 and find the zone that contains it. We supply that location, compare it with each zone's polygon, and only keep the zone rows that contain the point. Notice how this is just a filter!

In SQL, the condition goes in WHERE: st.within($pickup, boundary). The $pickup parameter is the point we're looking up, and boundary is the geometry column in the zone table storing the boundaries of the different zones. The st. prefix identifies spatial functions.^[4]

The other tabs in the example below show off different spatial conditions. Nearby points finds pickups within 400 meters of a location you can move on the map. Routes crossing zones searches the zone table for polygons that the route intersects. Both of these are filters as well!

-- $pickup: point
SELECT location_id, zone
FROM taxi_zones
WHERE st.within(
  $pickup,
  boundary
);

Spatial join

The point-in-zone filter takes just one location as a query parameter. But what if we want to look up the zone for every pickup in a table?

A spatial join uses the same st.within condition, but both inputs now come from tables. Each matching pair gives us a row with the trip ID from one table and the zone name from the other.

The example starts with nine pickups and three selected zones. All nineteen zones are available: try changing the selection or moving a pickup to see how it changes the result. Note that pickups outside the selected zones won't appear in the results.

As humans, we can simply look at the map and see which zone each pickup is in. A spatial join gives us that result as rows we can use in the rest of a query.

Spatial aggregate

The spatial join tells us which zone each pickup is in. Now we can use that result to answer a question about fares: how much did passengers pay by card in each zone?

We can filter the joined rows to keep card payments, then group them by zone and sum their fares. That combines a spatial join with an ordinary filter and aggregate, all in one SQL query.

For this example, we'll use a larger set of 48 illustrative trips in and around four selected taxi zones. Trips outside those zones don't match the join, so they don't contribute to the totals. Change the payment method or minimum fare to filter the remaining trips, then play the aggregation to follow their fares into the zone totals.

SELECT zones.zone, COUNT(*) AS trip_count,
       SUM(trips.fare_amount) AS total_fares
FROM trips
JOIN taxi_zones AS zones
  ON st.within(trips.pickup, zones.boundary)
WHERE zones.location_id IN (234, 79, 125, 261)
  AND trips.payment_method = 'card'
  AND trips.fare_amount >= 0
GROUP BY zones.zone
ORDER BY zones.zone;

What's inside the geometry column?

A common way to store geometry is Well-Known Binary, or WKB. It stores each geometry as one binary value containing its type and coordinates. For example, a path is represented by a point count followed by the path's (x, y) pairs. This gives systems a common format for exchanging complete geometries, letting a single binary column hold points, paths, and polygons.

However, suppose we only want to read every x coordinate in the geometry column. Later, we're going to use each geometry's minimum and maximum x values to find the left and right edges of its bounding box, so this is not an uncommon query.

WKB stores each x beside its corresponding y, with type identifiers and counts between coordinate sequences. The reader uses those identifiers and counts to locate each sequence. It can skip decoding the y values, but CPUs fetch memory in blocks, so reading the x values also brings in neighboring y bytes.

This is a familiar tradeoff in database design. A row store keeps each row's values together, much as WKB keeps a geometry's type, counts, and coordinates together. That's useful when a query retrieves complete rows. But analytics often reads just a few columns across many rows. When those values share a storage block with other columns, reading that block brings in values the query doesn't need.

A column store, on the other hand, keeps each column's values together, letting the query read only the columns it needs. Here are a few of the benefits for analytical queries:^[5]

• Keeping a column's values together improves locality: each memory block holds more values the query needs, so scans fetch fewer blocks.
• Values in a column often share patterns that compress well, reducing storage space and the number of bytes read.
• SIMD instructions can process several adjacent numeric values at once.^[6]
• Multiple CPU cores can scan separate ranges of a column, each reading adjacent values.^[7]

Essentially, placing WKB values in a column doesn't make the geometry inside them columnar. To make the coordinates columnar as well, we can store x and y in separate numeric arrays. Then the engine can scan every x coordinate without reading the y array.

Because these are ordinary numeric arrays, we can use the database's existing compression techniques on them. But to get the benefit during computation, the geometry operations need to work with that layout too.

Coordinates are columns too

For points, we only need two arrays: one for the x coordinates and another for the y coordinates. Row i gets its point from x[i] and y[i].

A path, or more formally a LineString, has a variable number of points, so we also need to know which coordinates belong to each path. We keep the x and y arrays and add integer offsets that mark where each path starts and ends. A four-point path uses offsets [0, 4], which means "read coordinates starting at position 0, stopping before position 4". If the next path has three points, the column's offsets become [0, 4, 7].

Polygons add another level: coordinates form rings, and rings form polygons. The first ring encloses the area; any further rings cut holes out of it. Each ring repeats its first coordinate at the end to close the boundary.

Because a polygon can have several rings, we need offsets at both levels: ring offsets locate coordinates in the x and y arrays, and polygon offsets locate rings. In the example below, we use a small four-vertex polygon on a 0-to-10 grid so we can follow every coordinate. Select or move a vertex to see where it lives in the arrays and in WKB:

The starting polygon has four vertices. Repeating the first one closes its ring (A, B, C, D, A), giving us five stored coordinate pairs. The ring offsets [0, 5] mark that coordinate range, and the polygon offsets [0, 1] say that this polygon contains one ring.

Turn on Include hole to add a second ring (E, F, G, E). Now we have nine coordinate pairs: the ring offsets become [0, 5, 9], and the polygon offsets become [0, 2]. Moving a vertex changes its coordinates but leaves these offsets alone, because the number of coordinates in each ring hasn't changed.

The offsets do more than keep the rings together. They also let us go straight to part of a geometry: the hole starts at coordinate index 5, so its third vertex is at index 7. We don't have to scan through the exterior ring first! This is random access within a geometry. With WKB, a count precedes each ring, so we have to follow the counts of the earlier rings to locate a later one.

SQL still treats the whole polygon as one value. Underneath, it's a list of rings, and each ring is a list of coordinate pairs stored in the x and y arrays. The GeoArrow specification describes this layout.

Keeping offsets in their own arrays lets us compress them without changing the coordinates. The demo above starts with 32-bit offsets, which take four bytes each. With the hole included, the largest offset is only nine, so each offset fits in one byte. Switch to 8-bit offsets to shrink the offset arrays from 20 bytes to 5 (though do note that the x and y arrays stay the same size).^[8]

Because these are numeric arrays, we can additionally use the database's existing integer and floating-point compression techniques.

Reject the obvious misses

Suppose a taxi pickup point is further north than an entire zone. Comparing the point's y coordinate with ymax, the top edge of the zone's bounding box, is enough to rule out a match, without inspecting any of the rings.

The same reasoning works in all four directions. A point outside the bounding box must also be outside the polygon. At most four coordinate comparisons can spare us a point-in-polygon check.

The edge case (as shown above) is that a point inside the box may still be inside one of the polygon's holes, or in the area between the polygon and the bounding box. So checking a point against a zone's bounding box may result in a false positive (the box test passes even though the point is outside the zone), so we still need to do a point-in-polygon check.

Don't check every pair

So far, we've used bounding boxes to skip reading zones that cannot contain a pickup. But in the spatial join from part one, we need to find the matching zones for every pickup in a table.

An R-tree can make these lookups faster by ruling out groups of zones at once. It applies the same bounding-box idea at several levels: we put the zone boxes into groups, enclose each group in another box, and repeat until one box encloses the whole tree.

If a pickup is outside a group's box, none of the zones in that group can contain it, so we can skip their individual boxes too! Otherwise, we continue searching inside that group.

We can build the tree once for the zone table and reuse it for every pickup in the join.^[12]

The example below uses a small tree over all nineteen taxi zones. Its groups are separate from the storage chunks above. Try moving pickup 1042, then play the search to see which groups we skip and which zones still need a polygon check.

At the starting position, the group boxes let us skip seven of the nineteen zone boxes. Checking the remaining twelve leaves two candidates, Flatiron and Union Square, whose boxes contain the pickup.^[13]

Notice how Flatiron's box contains the pickup even though its polygon doesn't. This is the same kind of false positive we saw earlier. Checking both polygons leaves just (1042, "Union Sq") in the result.

For the remaining polygon checks, we can work directly with our coordinate arrays. The polygon offsets tell us which rings to read, and the ring offsets tell us where their coordinates are in the x and y arrays. That gives the containment test the edges it needs to check.^[14]

Columnar geo in Spiral

In most systems, a polygon is still just one opaque value in a column. Spiral is able to take advantage of the internal structure of the values, allowing it to store, compress, query, and compute over geospatial data much more efficiently and effectively.

If you're working with geospatial data, get in touch! We'd love to hear what you're trying to query, and we can likely help you do it faster.