Tiling PostGIS Geometries

15 Aug 2013

ZIP Code boundaries in Alaska tend to be rather large. Some contain over a half-million vertices, and span multiple UTM zones.

Alaska ZIPs

It’s expensive to process intersections with these geometries. This problem compounds itself: because the complex geometries are so cover such a large area, their bounding boxes are huge. The largest features have a bounding box that covers the entire state, and will be returned by a bounding box intersection with any other feature:

Alaska ZIP Bounding Boxes

A solution for many processing situations is to split the largest polygons into smaller polygons, increasing the size of the spatial index but allowing it to return fewer polygons that don’t actually intersect. Just breaking apart multipolygons into polygons helps, and in some cases, this can make downstream processing hundreds of times faster. Still, it’s possible to do better. The code below will tile an input geometry according to a specified maximum dimension.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
CREATE OR REPLACE FUNCTION Tiler(size float, geom Geometry) RETURNS
SETOF geometry AS
$$
DECLARE
xlength float;
ylength float;
midpoint float;
splitline Geometry;
splitpoly Geometry;
g Geometry;
s Geometry;
i int;
BEGIN
IF ST_IsCollection(geom) THEN
    FOR i IN SELECT generate_series(1, ST_NumGeometries(geom)) LOOP
        FOR s IN SELECT Tiler(size, ST_GeometryN(geom, i)) LOOP
            RETURN NEXT s;
        END LOOP;
    END LOOP;
ELSE
    xlength := ST_XMax(geom) - ST_XMin(geom);
    ylength := ST_YMax(geom) - ST_YMin(geom);
    IF xlength > size THEN
        midpoint = ST_XMin(geom) + 0.5*xlength;
        splitline = ST_MakeLine(ST_MakePoint(midpoint, -90),
                                ST_MakePoint(midpoint,  90));
    ELSE
        IF ylength > size THEN
            midpoint = ST_YMin(geom) + 0.5*ylength;
            splitline = ST_MakeLine(ST_MakePoint(-180, midpoint),
                                    ST_MakePoint(180,  midpoint));
        END IF;
    END IF;
 
    IF xlength > size OR ylength > size THEN
        splitline := ST_SetSRID(splitline, ST_SRID(geom));
        splitpoly = ST_Split(geom, splitline);
        FOR g IN SELECT (ST_Dump(splitpoly)).geom LOOP
            FOR s IN SELECT Tiler(size, g) LOOP
                RETURN NEXT s;
            END LOOP;
        END LOOP ;
    ELSE
        RETURN NEXT geom;
    END IF;
END IF;
RETURN;
END;
$$
language plpgsql;

Below is the output from the function with a tile dimension of 1 degree latitude/longitude.

Alaska ZIPs (Tiled)

The bounding boxes of the tiled features are below.

Alaska ZIP Bounding Boxes (Tiled)

Published on 15 Aug 2013