Fast(er) voronoi polygons in PostGIS

create_voronoi_function.sql

-- plpython needs to be enabled. Furthermore scipy needs to be installed in the python used. 
-- CREATE LANGUAGE plpythonu;

CREATE OR REPLACE FUNCTION septima.voronoi(geom geometry, boundingpointdistance double precision)
  RETURNS SETOF geometry_dump AS
$BODY$
from scipy.spatial import Voronoi
import numpy

class GeoVoronoi:
    def __init__(self, points, boundingpointdistance = None):
        """Points are list tuples of (x,y)
        boundingpointdistance: Expand points bbox by this distance and add four points forming an enclosing box"""
        self.points = numpy.array(points)
        self.vor = None

        if boundingpointdistance:
            xmin, ymin = self.points.min(0)
            xmax, ymax = self.points.max(0)
            dx = dy = float(boundingpointdistance)
            xmin, ymin = (xmin - dx, ymin - dy)
            xmax, ymax = (xmax + dx, ymax + dy)
            # points
            ll = [xmin, ymin]
            ul = [xmin, ymax]
            ur = [xmax, ymax]
            lr = [xmax, ymin]
            plpy.info( "Adding enclosing points: ", ll, ul, ur, lr )
            self.points = numpy.append(self.points, [ll, ul, ur, lr], axis=0 )

    def process(self):
        self.vor = Voronoi( self.points )

    def polygons(self):
        # See http://docs.scipy.org/doc/scipy/reference/tutorial/spatial.html
        if not self.vor:
            self.process()
        # Get region belonging to each input point
        npregions = numpy.asarray( self.vor.regions )
        for pointix, region in enumerate( npregions[ self.vor.point_region ] ):
            region = numpy.asarray( region )
            if numpy.all(region >= 0):
                # Closed region
                yield self.vor.vertices[ region ]
            else:
                # region has ridges extending to infinity
                plpy.info( "Warning: infinite regions are not implemented" )

# Main function

result = plpy.execute("SELECT ST_srid('{0}'::geometry) as srid".format( geom ))
srid = result[0]['srid']

sql =   "SELECT st_x(geom) as x, st_y(geom) as y "\
        "FROM (SELECT (st_dumppoints('{0}'::geometry)).geom) as foo".format( geom )
input = plpy.execute(sql)
# input[i]['x'] holds x coord of point i

points = [ (row['x'], row['y']) for row in input  ]

gvor = GeoVoronoi( points, boundingpointdistance )
gvor.process()
part = 0
for p in gvor.polygons():
    part = part + 1
    wkt = "MULTIPOINT({0})".format( ','.join([ str(c[0]) + " " + str(c[1]) for c in p]) )
    #plpy.debug(wkt)
    sql = "SELECT st_convexhull(st_geometryfromtext('%s', %d)) as geom" % (wkt,srid)
    #arr = "ARRAY[ {0} ]".format( ','.join([ "st_makepoint(" + str(c[0]) + "," + str(c[1] +")") for c in p]) )
    #sql = "SELECT st_convexhull(st_setsrid(st_makeline({0}), {1})) as voronoi_region".format([wkt,srid])
    #plpy.debug(sql)
    result = plpy.execute(sql)
    yield ( [part] , result[0]['geom'])



$BODY$
LANGUAGE plpythonu 
IMMUTABLE
PARALLEL SAFE
COST 10000000;

test_voronoi_function.md

Voronoi in PostGIS

This code was implemented some years ago when I found the native PostGIS voronoi function to be too slow on my data. At that time this custom function was way faster than the PostGIS one on larger datasets. I have done any comparisons on newer PostGIS versions, so things may have changed since then.

Also - this code solved my problems - it may not work for you. Use it on your own responsibility :-)

Examples

Function returns set of geometry_dump as known from PostGIS function st_dump.

If the second param is null the voronoi cells along the edge of the diagram will be unbounded. These are not returned by the function.

If the second param is not null, the bbox of the points is expanded by this value and four points added at the corners of the expanded bbox. this ensures that the original points are all enclosed in finite cells.

SET client_min_messages TO DEBUG;
SELECT * FROM septima.voronoi('MULTIPOINT ( 2 0, 2 3, 0 2, 2 0 )' :: geometry, 100);

SELECT path[1], st_astext(geom ::geometry(polygon)) FROM
(SELECT * FROM septima.voronoi('MULTIPOINT ( 2 0, 2 3, 0 2, 2 0 )' :: geometry, 100)) as foo;

Input must be one single geometry. Calculating voronoi for a set of point records can be done by using st_collect.

--CREATE TABLE test_voronoi AS
SELECT path[1], geom ::geometry(polygon,25832) 
FROM septima.voronoi(
	(
		SELECT st_collect(geom) as mpoint 
		FROM adresses 
		WHERE municipalitycode = '0101'
	)
)