ST_Subdivide(): a quick intro

About ST_Subdivide()

• this function will receive an input geometry (may well be a very huge one).
• each Linestring or Polygon found within the input geometry will be then processed:
  • all Linestrings or Polygons using a number of Vertices lesser or equal than the given threshold will be returned as they are.
  • but all Linestrings or Polygons using an exceeding number of Vertices will be recursively split, until they'll be represented by a collection of elementary parts using no more than the required number of vertices.
• at the end of the process the function will return a collection containing all elementary parts (a MultiLinestring or a MultiPolygon depending on the nature of the input geometry).

Basic Examples

SQL query	Visual Sample
SELECT ST_Subdivide(geometry) FROM tuscany; in this first call no optional argument is specified, and thus ST_Subdivide() will implicitly assume the standard threshold of max 128 Vertices. As the side figure shows, the returned result is a collection of about 4,300+ elementary parts.
SELECT ST_Subdivide(geometry, 512) FROM tuscany; in this second call the optional argument is explicitly set, and ST_Subdivide() is now assuming a threshold of max 512 Vertices. As the side figure shows, the returned result is a collection of about 1,800+ elementary parts.
SELECT ST_Subdivide(geometry, 2048) FROM tuscany; in this third and final call the threshold is set to max 2,048 Vertices. As shown by the figure, the returned result is now a collection containing just 1,200+ elementary parts. Short conclusion ST_Subdivide() is effectively capable to subdivide some nasty huge Geometry into a collection of many smaller parts much more reasonable to be handled. And it's reasonably flexible and easy to be customized as required; you just have to set the most appropriate max. Vertices threshold.

SELECT ST_IsValid(ST_Subdivide(geometry)) FROM tuscany;
---------
0

Advanced Examples

• once again, the same administrative boundary of Tuscany we've already used before.
• in the role of the POINT dataset we'll use the civici Shapefile (House Numbers) freely available for download from here (search for Grafo Stradale - Road Network).
  This too is a huge dataset containing about 1.5 million points.

A first naive attempt

SELECT Count(*)
FROM civici AS c
LEFT JOIN tuscany AS t ON (ST_Intersects(c.geometry, t.geometry) = 1);


SELECT Count(*)
FROM civici AS c
LEFT JOIN tuscany AS t ON (ST_Intersects(c.geometry, t.geometry) = 1 
     AND c.rowid IN (
         SELECT rowid FROM SpatialIndex
         WHERE f_table_name = 'civici' AND search_frame = t.geometry));

• I was forced to prematurely abort the first test when I discovered that during the first half an hour it had made only a very little progress.
  Accordingly to my estimated extrapolation, waiting for completion will had presumably required between 4 and 6 hours.
  Definitely this isn't a practical solution for our problem.
• the second test was, if possible, even slower than the first one.
  And there is very good reason explaining for this; all Points to be checked fall inside the BBOX of Tuscany.
  So the Spatial Index isn't of any possible help.
  And even worst, querying the Spatial Index imposes a further overhead, thus requiring a longer time for absolutely nothing.
• Short conclusion: computing many million times the intersections between a Point and a Polygon defined by some hundredth thousands Vertices is a damn hard computational problem.
  And consequently it requires a very long time, so long to be absolutely unpractical.

A second attempt based on ST_Subdivide() and VirtualElementary

CREATE TABLE xx (id INTEGER PRIMARY KEY);
SELECT AddGeometryColumn('xx', 'geom', 3003, 'MULTIPOLYGON', 'XY');
INSERT INTO xx VALUES (1, NULL);

• it has just two columns, a Primary Key and a MultiPolygon Geometry.
• and it's intended to store just a single row identified by id=1

UPDATE xx SET geom = (SELECT ST_Subdivide(geometry, 2048) FROM tuscany) WHERE id = 1;
SELECT Count(*)
FROM ElementaryGeometries AS t
JOIN civici AS c ON (ST_Intersects(c.geometry, t.geometry) = 1
     AND c.rowid IN (
         SELECT rowid FROM SpatialIndex
         WHERE f_table_name = 'civici' AND search_frame = t.geometry))
WHERE f_table_name = 'xx' AND origin_rowid = 1;
----------------
1480991 (in 38.094 seconds)


UPDATE xx SET geom = (SELECT ST_Subdivide(geometry, 512) FROM tuscany) WHERE id = 1;
SELECT Count(*)
FROM ElementaryGeometries AS t
JOIN civici AS c ON (ST_Intersects(c.geometry, t.geometry) = 1
     AND c.rowid IN (
         SELECT rowid FROM SpatialIndex
         WHERE f_table_name = 'civici' AND search_frame = t.geometry))
WHERE f_table_name = 'xx' AND origin_rowid = 1;
----------------
1480991 (in 12.617 seconds)


UPDATE xx SET geom = (SELECT ST_Subdivide(geometry) FROM tuscany) WHERE id = 1;
SELECT Count(*)
FROM ElementaryGeometries AS t
JOIN civici AS c ON (ST_Intersects(c.geometry, t.geometry) = 1
     AND c.rowid IN (
         SELECT rowid FROM SpatialIndex
         WHERE f_table_name = 'civici' AND search_frame = t.geometry))
WHERE f_table_name = 'xx' AND origin_rowid = 1;
----------------
1480991 (in 9.534 seconds)

• ST_Subdivide definitely is the tool of the trade we were searching for.
  It effectively allows to resolve a computationally hard problem in a surprisingly quick time.
  Passing from many long hours to just a handful of seconds is a simply astonishing improvement.
  
• subdividing the initial huge Polygon into many simpler/smaller parts has several beneficial effects:
  • the time required to compute a point-to-polygon intersection quickly increases with the number of Vertices.
    Smaller polygons are obviously much more faster to be evaluated than bigger ones.
  • smaller polygons doesn't shadow the Spatial Index, that can be efficiently queried.
• using VirtualElementary in order to process all sub-parts one by one surely imposes some further complexity to our SQL queries.
  But at the end of the day it confirms to be a very efficient mechanism.
• Note: SpatiaLite is not PostGIS, and SQLite3 is not PostgreSQL; their respective internal architectures significantly differ under many aspects.
  More specifically, SQLite is completely unable to support arrays or multi-row results, so using VirtualElementary is an absolutely necessary workaround.

A final experimental attempt based on a pure SQL approach

WITH RECURSIVE magic(n, geom) AS (
    VALUES(1, (SELECT ST_Subdivide(geometry, 2048) FROM tuscany))
    UNION ALL
    SELECT n + 1, geom
    FROM geom_ind
    WHERE ST_GeometryN(geom, n) IS NOT NULL
)
SELECT n, ST_GeometryN(geom, n) FROM magic;

• we simply have to write a Recursive SQL Query in order to get a multi-row resultset, with a distinct row for each elementary part.
• it effectively works, but seems to be slower than ElementaryGeometries.
• and even more relevant, I was unable to write a complete Query determining all intersections between Points and Polygons in a reasonable time.
  I can simply suppose that the Query Optimizer of SQLite will end up in ignoring the Spatial Index, thus completely vanishing the usefulness of this approach.
  May well be that someone else could succeed where I failed; any further third-party investigation on this specific area is surely welcome.