Instantly share code, notes, and snippets.

# AsgerPetersen/create_voronoi_function.sql

Last active November 17, 2021 15:41
Show Gist options
• Save AsgerPetersen/a439ca7aef7efde094f8f524efeaf150 to your computer and use it in GitHub Desktop.
Fast(er) voronoi polygons in PostGIS
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 -- 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;

# 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