Monday, April 9, 2012

Faster Voronoi Diagrams in PostGIS

A couple years back, I posted a quick and dirty function for generating voronoi diagrams in PostGIS. As one of the commenters pointed out, this is not the world's most efficient algorithm. So I started looking into implementing a better one myself, but once again, being a virtuous programmer, I borrowed someone else's code, and modified it to work with PostGIS.

This one's in Python, not PL/pgSQL, so plan appropriately.  It is called similarly to my previous version:

select * from voronoi('table_name', 'geom_column') as (id integer, the_geom geometry);

You no longer need to specify an id column, but the downside is that there's no longer a relationship between the id of the input point, and the id of the output polygon. On the other hand, it's much, much, much more efficient. At least a couple orders of magnitude. In testing, my 40,000 points of input only took 20 seconds to calculate, whereas previously it took... well, I don't know, because I didn't want to wait that long.

Code is here:  No guarantees it actually does anything correctly, but my limited testing seems to work. Enjoy.

UPDATE: the code linked above only works with Python 2.x. If you look through the comments, you'll see there's a link to a version of my code for 3.x, but it still has a couple of bugs which I have not had time to investigate.

UPDATE 2: See the comments below about the QGIS code also producing incorrect output with this function. I'll try to investigate in my copious free time.

UPDATE 3: Matthias has published a version that corrects the problems in my initial code. His dropbox link is here:  I've also published it as GIST, here: