Subdivide and Conquer Any Geometry

Let’s write a geospatial function that BigQuery does not currently provide: analogue of PostGIS’s ST_Subdivide. We’ll use SQL UDF, and learn a few tricks along the way.

ST_Subdivide breaks a shape into pieces. There are various reasons to do it, sometimes it is needed to ensure each shape has low number of vertices, sometimes to make all the shapes below some threshold of size or area. It takes a (usually) spatially large and complex shape as argument and returns an array of smaller and simpler pieces.

We’ll just divide input into a fixed number of pieces, to avoid extra complexity of bringing recursive calls into the picture. Thanks StackExchange GIS user Max for asking this question.

First, let’s define overall strategy. A simple approach seems to be picking some point at “the middle” of the shape, and splitting the shape using “vertical” and “horizontal” lines passing through this point into four pieces. We can repeat this to further split each of the resulting geometries.

To split, we’ll construct the shapes defining each of the quadrant, then create four shapes as intersections of the original shape with each quadrant:

[ST_Intesection(geometry, top_left), ST_Intesection(geometry, top_right), ...]

Well, we live on a spherical Earth, so it is a bit more complicated. Let’s look into details.

Division point

The two obvious choices are either centroid

st_x(st_centroid(g)) x, st_y(st_centroid(g)) y 

or middle of the bounding box:

select (bb.xmin + bb.xmax) / 2 as x, (bb.ymin + bb.ymax) / 2 as y
from (select st_boundingbox(g)) bb 

We’ll use centroid below, but the middle of bounding box is a good choice too, and might be preferable in some cases.

Quadrants

We could probably define the four quadrants directly, but it seems easier to split globe into hemispheres instead: west, east, north and south. Pairwise intersections of those give us four quarters of the globe — ”north-west”, etc.

To construct west and east hemispheres, we’ll split the globe using the plane that passes through the selected point (x, y) and poles, shown below. We cannot connect north and south poles directly, as an edge between two antipodal points is undefined. Thus we add another point “behind” the globe. A point antipodal to our origin point works, its coordinates are (x+180, -y). North pole is (*, 90) where * is any value — its longitude does not matter, and south pole is (*, -90).

Splitting the globe into west and east hemispheres

Note that the same four vertices describe both the west and the east hemispheres! To reliably distinguish these two we have to use oriented versions of the polygon constructor. The polygon should be to the left as we walk along the line describing the boundary, so starting with our point if we move to the North, then antipodal point, South and back — we get West hemisphere, otherwise — East one. Unlike WKT, ST_MakePolygon[Oriented] does not require repeating the implicit last vertex, so we can use

st_makepolygonoriented([st_makeline(
    [st_geogpoint(x, y), st_geogpoint(x, 90), 
     st_geogpoint(x+180, -y), st_geogpoint(x, -90)])]) as west

The east hemisphere will be the same, but we visit the vertices in the opposite order: our center point, South pole, antipodal point, North pole.

Now let’s move to defining north and south hemispheres. Here our great circle crosses equator at longitudes x±90, and then meets our antipodal point again:

Splitting the globe into north and south hemispheres

North hemisphere will be

st_makepolygonoriented([st_makeline(
    [st_geogpoint(x, y), st_geogpoint(x+90, 0), 
     st_geogpoint(x+180, -y), st_geogpoint(x-90, 0)])]) as north

Let’s put it all together:

create or replace function tmp.st_subdivide4(g GEOGRAPHY) 
returns ARRAY<GEOGRAPHY>
as ((
  with centroid as (
    select st_x(st_centroid(g)) x, st_y(st_centroid(g)) y 
  ),
  hemis as (
    select 
      st_makepolygonoriented([st_makeline(
          [st_geogpoint(x, y), st_geogpoint(x, 90), 
           st_geogpoint(x+180, -y), st_geogpoint(x, -90)])]) west,
      st_makepolygonoriented([st_makeline(
          [st_geogpoint(x, y), st_geogpoint(x, -90), 
           st_geogpoint(x+180, -y), st_geogpoint(x, 90)])]) east,
      st_makepolygonoriented([st_makeline(
          [st_geogpoint(x, y), st_geogpoint(x+90, 0),
           st_geogpoint(x+180, -y), st_geogpoint(x-90, 0)])]) north,
      st_makepolygonoriented([st_makeline(
          [st_geogpoint(x, y), st_geogpoint(x-90, 0),
           st_geogpoint(x+180, -y), st_geogpoint(x+90, 0)])]) south
    from centroid
  )
  select [st_intersection(g, st_intersection(west, north)), 
          st_intersection(g, st_intersection(east, north)),
          st_intersection(g, st_intersection(west, south)),
          st_intersection(g, st_intersection(east, south))]
  from hemis
));

If we need to divide further:

create or replace function tmp.st_subdivide16(g GEOGRAPHY) 
returns ARRAY<GEOGRAPHY>
as ((
  select array_concat_agg(tmp.st_subdivide4(g1))
  from unnest(tmp.st_subdivide4(g)) g1
));

Test

Let’s visualize the result using GeoViz, adding rand() to each piece so we can color the resulting pieces differently:

select g, rand() r 
from unnest(tmp.st_subdivide16(st_geogfromgeojson(...))) g

Large shape divided into 16 pieces, rendered by BigQuery GeoViz

Caveat

This code does not work if the centroid of the input shape happens to be at the North or South pole. A more reliable version should check this condition, and use another set of hemispheres, e.g. one pair with boundaries along 0-and 180-meridian, and one pair with boundaries along 90 and -90-meridian.