(Update: A faster version of this code is here: Faster Voronoi Diagrams in PostGIS)
I found myself in a situation where I needed to create a Voronoi diagram from a set of points using PostGIS. This isn't built in, but a little googling turned up a function written by Mike Leahy. I adapted it and made it a bit more generic, so it shouldn't require any editing to make it work for many (most?) cases.
Call it as: select * from voronoi(table_name text,table_primary_key text,geom_col text) as (id id_type, the_geom geometry)
For example:
select * from voronoi('oregon_cities', 'id', 'the_geom_centroids') as (id integer, the_geom geometry);
The 'as' clause is required because we don't know and don't want to assume what type your primary key might be. But, because we're returning a record type, the query planner wants to know what type to expect. If you're always using the same type, say an integer for your primary key, then it would be easy to define a type thusly:
create type my_type as (id integer, the_geom geometry);
Then change the code at line 24 from 'rec record;' to 'rec record%my_type;' and you will no longer need the 'as' clause.
Code at this URL: http://db.tt/aInuiJ2. Enjoy.
hello!!
ReplyDeletethanks for your script! I have been using it successfully on postgres 9.0. I just had to modify the script using the st_union() function instead of the geomunion()one. Hope it will be usefull for you!
Sorry for my english, I'm a french customer.
Hi!
ReplyDeleteThis is awesome! I'm very glad to see a free option to the ESRI Voronoi utility. Works very well, except I've noticed that sometimes, polygons are missing from the results. Any idea why that happens?
Thanks!
Hey I am trying this function for a dataset with 10.000 points, it fills up more than 50 Gigabytes of disk space and then is cancelled because lackof disk space. Is there any idea how to decrease the disk psace use for this function a bit?
ReplyDeleteThanks!
It's true it's not the most space efficient algorithm. The only suggestion I can think of off the top of my head is to do it "batches" of a thousand or so, grouped by the_geom, and then merge the output of each batch together somehow.
DeleteThe algorithm is pretty naïve. I suppose it wouldn't be too challenging to implement Fortune's algorithm, but I'd probably want to modify an existing implementation (such as this one from QGIS: https://svn.osgeo.org/qgis/trunk/qgis/python/plugins/fTools/tools/voronoi.py) to work inside postgres.
That actually doesn't look like too much of a challenge, and if I get a chance, I might try to do that.
OK, I did exactly that. I'll have a post up on it shortly. I just did a test on 7900 points and it took 7.5 seconds, and 40,000 points took 20 seconds.
DeleteWas wondering if you could make your modified code available? I've come across the other Voronoi example on BostonGIS but it dependso n PL/R and the R stats packages. I'd like to work with your example as it sounds like it doesn't have that dependency.
ReplyDeleteThank
Hi Paul,
DeleteThe code is available. See http://geogeek.garnix.org/2012/04/faster-voronoi-diagrams-in-postgis.html
If you look through the comments you'll see that it only works with Python 2.x – the 3.x version still has some bugs I haven't had time to fix.
After trying the plpythonu Voronoi function and producing artefacts and overlapping and self-overlapping voronoi polygons all the time (same happens if you use the ftools voronoi tool in QGIS, as the plpython code is derived from this), I tried out the PL/R voronoi function from http://punkish.org/Voronoi-Diagrams-In-PostGIS and found out that after some small adaptions to postgis 2 (make srid()-->st_srid() etc), the function produced a flawless set of voronoi polygons for 1000 points in about 13 seconds on a Core2duo CPU.
ReplyDelete