ZIP Code boundaries in Alaska tend to be rather large. Some contain over a half-million vertices, and span multiple UTM zones.
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:
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.
The bounding boxes of the tiled features are below.