Break down one extremely large Polygon into many smaller polygons using BigQuery only

Question

I have a giant polygon which is part of the UK's open source Flood Zone 2 dataset, found here.

At the moment, it is breaking our data pipeline due to its sheer size. As far as I am aware, it is the only one like this, and it wasn't historically always like this.

All of our data conversion happens in BigQuery, therefore we must come up with an SQL statement that is capable of breaking this large polygon into hundreds or thousands of smaller polygons.

I think normal PostGIS has ST_Subdivide function for this. However, this function does not exist in BigQuery. I am not sure of an alternative. We run this pipeline monthly which downloads the data from an API before running the conversion, so it has to be reusable (i.e. locally converting it once with QGIS won't be a long term solution).

How can I do this?

This is a link to download the troublesome table row as JSON.

Answer 1

It looked like fun exercise. I created a SQL function that divides a given geography into 4 sections, with centroid as the common point. You can call it a fixed number of times or in a loop until you reach the desired number of vertices.

The idea is to split the globe into 4 pieces using two orthogonal planes. Both planes pass through the centroid of our shape and its antipodal point. One plane passes through poles, dividing Earth into "west" and "east" part around the centroid of our shape, another orthogonal to it divides into "north" and "south" parts.

We then intersect the given polygon with these hemispheres to produce four smaller pieces.

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
  ),
  halfs as (
    select 
      north, south, 
      st_intersection(g, west) as west, 
      st_intersection(g, east) as east
    from hemis
  )
  select [st_intersection(west, north), st_intersection(west, south), 
          st_intersection(east, north), st_intersection(east, south)]
  from halfs
));

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
));

I visualized the result using (rand added so I can color the resulting pieces differently):

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

P.S. I turned this function into a post, it has more explanations and slight query tweaks:

https://mentin.medium.com/subdivide-and-conquer-any-geometry-ca4f0a4b8491