GEOGRAPHY and spatial functions in BigQuery

BigQuery's GEOGRAPHY type is a set of points, lines and polygons on a sphere, with edges that follow great circles and distances in meters. ST_GEOGPOINT takes longitude first. This covers WKT and GeoJSON, polygon orientation, spatial joins and clustering.

A GEOGRAPHY value is a set of points, lines and polygons on the Earth’s surface, stored as longitude and latitude on a sphere. Edges between points are the shortest path over the sphere, not straight lines on a map, and every distance and area is in meters and square meters. The most common mistake is the argument order: ST_GEOGPOINT(longitude, latitude), x before y, which is the reverse of how people say coordinates out loud.

Longitude first

ST_GEOGPOINT(lng, lat) accepts any latitude between -90 and 90, so swapping the arguments often succeeds and gives a valid point in the wrong place:

SELECT
  city,
  ST_ASTEXT(ST_GEOGPOINT(lng, lat)) AS right_order,
  ST_ASTEXT(ST_GEOGPOINT(lat, lng)) AS swapped
FROM UNNEST([
  STRUCT('London' AS city, -0.1278 AS lng, 51.5074 AS lat),
  STRUCT('Quito', -78.4678, -0.1807)
])
ORDER BY city

London swapped becomes longitude 51.5, latitude -0.13, a point in the Indian Ocean off Somalia. Nothing fails. The error only appears when the swapped latitude is out of range, as for Tokyo:

SELECT ST_GEOGPOINT(35.6762, 139.6503) AS tokyo

The error message itself suggests checking the order. A longitude outside -180 to 180 isn’t an error: ST_GEOGPOINT(200, 10) is wrapped to POINT(-160 10). Check a few known cities after loading any latitude and longitude columns.

WKT, WKB and GeoJSON

BigQuery reads and writes three text and binary formats. ST_GEOGFROMTEXT parses WKT, ST_GEOGFROMGEOJSON parses GeoJSON, ST_GEOGFROMWKB parses WKB, and ST_GEOGFROM detects which of the text forms it was given. ST_ASTEXT, ST_ASGEOJSON and ST_ASBINARY go the other way. A GEOGRAPHY column exports to CSV or JSON as WKT.

SELECT
  ST_ASTEXT(ST_GEOGPOINT(-121, 41)) AS wkt,
  ST_ASGEOJSON(ST_GEOGPOINT(-121, 41)) AS geojson,
  ST_ASTEXT(ST_GEOGFROMGEOJSON('{"type": "Point", "coordinates": [-121, 41]}')) AS from_geojson

The GeoJSON functions take a geometry object only. A Feature or FeatureCollection is rejected:

SELECT ST_GEOGFROMGEOJSON('{"type": "Feature", "geometry": {"type": "Point", "coordinates": [1, 2]}, "properties": {}}') AS g

To load features, use a newline-delimited GeoJSON file with bq load --json_extension=GEOJSON, which turns geometry into a GEOGRAPHY column and each member of properties into its own column.

Some things the formats can’t carry. Three-dimensional coordinates (a Z in WKT, an altitude in GeoJSON) and M values are rejected, and WKT other than points, lines, polygons and their multi and collection forms is rejected. Parsing is strict about invalid shapes, so wrap untrusted input in SAFE.:

SELECT id, SAFE.ST_GEOGFROMTEXT(wkt) IS NULL AS is_bad
FROM UNNEST([
  STRUCT(1 AS id, 'POINT(1 2)' AS wkt),
  STRUCT(2, 'POINT(1 2 3)'),
  STRUCT(3, 'POLYGON((0 0, 0 10, 10 10, 10 0))')
])
ORDER BY id

Rows 2 and 3 fail for different reasons (a Z coordinate, a ring that isn’t closed). Without SAFE., the whole query fails on the first bad row.

GeoJSON export adds vertices

A GeoJSON geometry has planar edges and a GEOGRAPHY has spherical ones, so ST_ASGEOJSON adds points wherever needed to keep the planar line within 10 meters of the original. A two-point route comes out with hundreds of points:

WITH route AS (SELECT ST_GEOGFROMTEXT('LINESTRING(-118.4085 33.9416, -73.7781 40.6413)') AS g)
SELECT
  ST_NUMPOINTS(g) AS points_in_wkt,
  ARRAY_LENGTH(JSON_QUERY_ARRAY(PARSE_JSON(ST_ASGEOJSON(g), wide_number_mode => 'round'), '$.coordinates')) AS points_in_geojson
FROM route

We had to pass wide_number_mode => 'round' to PARSE_JSON here, because the coordinates BigQuery writes have more digits than a FLOAT64 round-trips.

Edges are great circles, distances are meters

A line from (0 60) to (30 60) is not drawn along the 60th parallel. It follows the great circle through both points, which bulges toward the pole. The point halfway along the parallel is 94.8 km from that line:

-- A line along the 60th parallel is not what BigQuery draws: edges are geodesics.
WITH line AS (SELECT ST_GEOGFROMTEXT('LINESTRING(0 60, 30 60)') AS g)
SELECT
  ROUND(ST_DISTANCE(g, ST_GEOGPOINT(15, 60))) AS meters_from_the_parallel,
  ROUND(ST_BOUNDINGBOX(g).ymax, 2) AS northmost_latitude,
  ST_DWITHIN(g, ST_GEOGPOINT(15, 60), 1000) AS within_1_km
FROM line

The line’s northmost latitude is 60.85. If your source data was drawn on a flat map (most shapefiles and GeoJSON are), long straight edges mean something different in BigQuery. Add vertices to long edges before loading, or read them with ST_GEOGFROMTEXT(wkt, planar => TRUE), which interprets the edges with planar map semantics. planar and oriented can’t both be TRUE.

Distances use a sphere, not the WGS84 ellipsoid. Google’s reference lists a use_spheroid argument on ST_DISTANCE, ST_AREA, ST_LENGTH and others, but it only accepts FALSE: TRUE fails with “does not yet implement use_spheroid=true”. ST_DISTANCE, ST_LENGTH, ST_PERIMETER and ST_AREA return meters or square meters, and ST_DWITHIN takes its distance in meters:

SELECT
  ROUND(ST_DISTANCE(ST_GEOGPOINT(-0.1278, 51.5074), ST_GEOGPOINT(2.3522, 48.8566)) / 1000, 1) AS london_paris_km,
  ROUND(ST_LENGTH(ST_GEOGFROMTEXT('LINESTRING(-118.4085 33.9416, -73.7781 40.6413)')) / 1000, 1) AS lax_jfk_km,
  ST_DWITHIN(ST_GEOGPOINT(-0.1278, 51.5074), ST_GEOGPOINT(2.3522, 48.8566), 350000) AS within_350_km

ST_DISTANCE between the two airports and ST_LENGTH of the line joining them are the same number. Because the Earth isn’t a sphere, a distance can differ from a survey on the ellipsoid by up to about half a percent.

Polygons: which side is the inside

On a sphere, a ring of points divides the Earth into two regions, and either could be “the polygon”. ST_GEOGFROMTEXT picks the smaller one by default. For polygons bigger than a hemisphere, pass oriented => TRUE, which makes the inside the side on your left as you walk the vertices in order:

WITH box AS (
  SELECT 'clockwise' AS ring, 'POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))' AS wkt
  UNION ALL
  SELECT 'counter-clockwise', 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'
)
SELECT
  ring,
  ROUND(ST_AREA(ST_GEOGFROMTEXT(wkt)) / 1e9) AS default_km2_thousands,
  ROUND(ST_AREA(ST_GEOGFROMTEXT(wkt, oriented => TRUE)) / 1e9) AS oriented_km2_thousands
FROM box
ORDER BY ring

Clockwise as drawn on a map, the ring is read as the whole Earth minus the 10 degree box: about 508.8 million square kilometers. GeoJSON never has this ambiguity, because it is planar. If an area comes out near 510 million km², the polygon was wound the wrong way.

Invalid polygons are errors, with messages like “Edge 0 crosses edge 2” or “Polygon loop should be a closed ring”:

SELECT ST_GEOGFROMTEXT('POLYGON((0 0, 10 10, 10 0, 0 10, 0 0))') AS bowtie

make_valid => TRUE repairs what it can. Here the bow tie becomes a MULTIPOLYGON of two triangles meeting at a new vertex:

SELECT ST_ASTEXT(ST_GEOGFROMTEXT('POLYGON((0 0, 10 10, 10 0, 0 10, 0 0))', make_valid => TRUE)) AS repaired

Predicates and spatial joins

ST_CONTAINS, ST_COVERS, ST_INTERSECTS, ST_WITHIN, ST_COVEREDBY, ST_TOUCHES, ST_EQUALS and ST_DWITHIN are the boolean predicates. The difference between ST_CONTAINS and ST_COVERS is the boundary: a point on the edge of a polygon is covered but not contained.

WITH zone AS (SELECT ST_GEOGFROMTEXT('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))') AS area)
SELECT
  point_name,
  ST_CONTAINS(area, p) AS st_contains,
  ST_COVERS(area, p) AS st_covers,
  ST_INTERSECTS(area, p) AS st_intersects
FROM zone, UNNEST([
  STRUCT('inside' AS point_name, ST_GEOGPOINT(5, 5) AS p),
  STRUCT('on the boundary', ST_GEOGPOINT(10, 5))
])
ORDER BY point_name

A spatial join is a join on one of these predicates:

WITH zones AS (
  SELECT 'north' AS zone, ST_GEOGFROMTEXT('POLYGON((0 5, 10 5, 10 10, 0 10, 0 5))') AS area
  UNION ALL
  SELECT 'south', ST_GEOGFROMTEXT('POLYGON((0 0, 10 0, 10 5, 0 5, 0 0))')
),
stores AS (
  SELECT 'Aldgate' AS store, ST_GEOGPOINT(3, 8) AS loc
  UNION ALL SELECT 'Brixton', ST_GEOGPOINT(4, 1)
  UNION ALL SELECT 'Camden', ST_GEOGPOINT(7, 9)
  UNION ALL SELECT 'Dalston', ST_GEOGPOINT(2, 4)
)
SELECT zone, COUNT(*) AS stores
FROM zones JOIN stores ON ST_CONTAINS(area, loc)
GROUP BY zone
ORDER BY zone

BigQuery optimizes spatial joins for INNER JOIN and CROSS JOIN with ST_DWITHIN, ST_INTERSECTS, ST_CONTAINS, ST_WITHIN, ST_COVERS, ST_COVEREDBY, ST_EQUALS and ST_TOUCHES. Outer joins, anti-joins and a negated predicate are not optimized, and ST_DWITHIN is optimized only when the distance is a constant. Google also advises storing GEOGRAPHY values in the table instead of building them with ST_GEOGPOINT in the join.

What you can’t do with a GEOGRAPHY

A GEOGRAPHY can’t be grouped, sorted, compared with = or used with DISTINCT. Use ST_EQUALS for equality and group on ST_ASTEXT(loc) or a geohash:

SELECT loc, COUNT(*) AS visits
FROM UNNEST([ST_GEOGPOINT(1, 1), ST_GEOGPOINT(1, 1), ST_GEOGPOINT(2, 2)]) AS loc
GROUP BY loc

SELECT ST_ASTEXT(loc) AS loc, COUNT(*) AS visits
FROM UNNEST([ST_GEOGPOINT(1, 1), ST_GEOGPOINT(1, 1), ST_GEOGPOINT(2, 2)]) AS loc
GROUP BY 1
ORDER BY 1

ORDER BY fails with “ORDER BY does not support expressions of type GEOGRAPHY” and = with “Equality is not defined for arguments of type GEOGRAPHY”. To combine geographies, use the aggregate ST_UNION_AGG. CAST(geo AS STRING) is also an error; use ST_ASTEXT.

Grid cells: geohash and S2

Two ways to turn a place into a groupable key. ST_GEOHASH(point, length) gives a base-32 string for a rectangle, and ST_GEOGPOINTFROMGEOHASH goes back to the cell’s center:

SELECT
  ST_GEOHASH(ST_GEOGPOINT(-0.1278, 51.5074), 6) AS geohash_6,
  ST_ASTEXT(ST_GEOGPOINTFROMGEOHASH('gcpvj0')) AS cell_center

BigQuery also has the S2 functions S2_CELLIDFROMPOINT(point, level => n) and S2_COVERINGCELLIDS(geo, ...), which return INT64 cell IDs. The covering function returns the cells that cover a shape, up to max_cells. You can use the cell IDs for your own bucketing and joins:

SELECT
  S2_CELLIDFROMPOINT(ST_GEOGPOINT(-0.1278, 51.5074), level => 10) AS cell_level_10,
  ARRAY_LENGTH(S2_COVERINGCELLIDS(ST_GEOGFROMTEXT('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'), max_cells => 4)) AS cells_covering_box

Google’s guide on grid systems for spatial analysis covers when a cell join beats a geometry join.

Clustering on a GEOGRAPHY column

A GEOGRAPHY column can be a clustering column but not a partitioning column. Clustering on it sorts blocks by their location, so a query with a spatial predicate can skip blocks far from the search area.

We tried this on a public table of 81 million places, clustered by its geometry column. The dry run is an upper bound and doesn’t show the saving:

SELECT COUNT(id) AS places
FROM `bigquery-public-data.overture_maps.place`

SELECT COUNT(id) AS places
FROM `bigquery-public-data.overture_maps.place`
WHERE ST_DWITHIN(geometry, ST_GEOGPOINT(-0.1278, 51.5074), 1000)

Adding the predicate raises the estimate to 6.35 GB because the query now also reads geometry. We then ran the filtered queries for real, with the cache off:

The query read 0.2% of the estimate. Pruning follows the size of the search area: a 1,000 km radius reads a fifth of the table. ST_DISTANCE(...) < n was pruned the same way as ST_DWITHIN in this test, though the predicate documentation only names ST_DWITHIN.

Two consequences. With on-demand pricing you pay for the 14.1 MB, not the estimate. But a cost limit is checked against the estimate: a maximum_bytes_billed of 4 GB made the query fail with “6354370560 or higher required”, although it would have processed 14 MB. See controlling costs and clustering.

Which columns to cluster on is covered in partitioning vs clustering. For geography the rule is simple: if you filter or join by location, cluster by the geography column.

Other functions worth knowing

• ST_BUFFER(geo, meters) and ST_BUFFERWITHTOLERANCE grow or shrink a shape.
• ST_SIMPLIFY(geo, tolerance_meters) removes vertices, and ST_SNAPTOGRID rounds them.
• ST_BOUNDINGBOX returns a STRUCT with xmin, ymin, xmax, ymax, not a geography. ST_ASTEXT on it is a type error. On a great-circle edge, ymax can be north of every vertex, as in the 60.85 above.
• ST_CENTROID, ST_CLOSESTPOINT, ST_INTERSECTION, ST_UNION, ST_DIFFERENCE and ST_UNION_AGG build new geographies. ST_ISEMPTY and ST_DIMENSION test them.
• ST_REGIONSTATS reads raster data through Earth Engine.