3

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.

1 Answer 1

2

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

image divided into 16 pieces

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

4
  • Wow this is perfect! Thank you so much
    – Max
    Feb 20 at 11:37
  • 1
    Please up vote the answer as accepted answer. It helps other users find solutions for similar problems. Feb 20 at 15:54
  • Sigh, this is one of those genius answers that deserves 10 upvotes and instead we get "user unregistered". May 8 at 10:38
  • Thanks @MichaelEntin! Jun 7 at 10:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.