[postgis-users] Voronoi tessellation

[Puneet Kishor]

Thanks Derek,

I was able to make this code work by fiddling with the `sprintf` lines. For some reason, those lines (there are a couple of them below) are breaking. I will take a more detailed look at this tomorrow, as this is going to be a fairly used function for me.


I also saw a PL/PgSQL version of this floating around, but that was really inefficient... I think that ends up doing a cartesian join because it takes forever to run.


Voronoi and Delaunay are useful, and might be worth adding them in core Postgis.


On Mar 5, 2012, at 10:55 PM, Derek Jones wrote:

> Here is what I have that AFAIR works fine. It's been a while since I used the code - I may be digging it out again soon though :-)
> 
> -------------
> 
> 
> 
> 
> create type voronoi as (id integer, polygon geometry);
> 
> create or replace function voronoi(text,text,text) returns setof voronoi as '
>    library(deldir)
> 
>    # select the point x/y coordinates into a data frame...
> 
>    points <- pg.spi.exec( sprintf("select x(%2$s) as x, y(%2$s) as y from %1$s;",arg1,arg2) )
> 
>    # calculate an approprate buffer distance (~10%):
> 
>    buffer = ( ( abs( max( points$x ) - min( points$x ) ) + abs( max( points$y ) - min ( points$y ) ) ) / 2 ) * (0.10)
> 
>    # get EWKB for the overall buffer of the convex hull for all points:
> 
>    buffer <- pg.spi.exec(sprintf("select buffer(convexhull(geomunion(%2$s)),%3$.6f) as ewkb from %1$s;",arg1,arg2,buffer))
> 
>    # the following use of deldir uses high precision and digits to prevent slivers between the output polygons, and uses
>    # a relatively large bounding box with four dummy points included to ensure that points in the peripheral areas of the
>    # dataset are appropriately enveloped by their corresponding polygons:
> 
>    voro = deldir ( points$x, points$y, digits=22, frac=0.00000000000000000000000001, list(ndx=2,ndy=2), rw=c ( min( points$x ) - abs( min( points$x ) - max( points$x ) ), max( points$x ) + abs( min( points$x ) - max( points$x ) ), min( points$y ) - abs( min( points$y ) - max( points$y ) ), max( points$y ) + abs( min( points$y ) - max( points$y ) ) ) )
> 
>    tiles = tile.list(voro)
> 
>    poly = array()
> 
>    id = array()
> 
>    p = 1
> 
>    for (i in 1:length(tiles))
>    {
>      tile = tiles[[i]]
> 
>      curpoly = "POLYGON(("
> 
>      for (j in 1:length(tile$x))
>      {
>        curpoly = sprintf("%s %.6f %.6f,",curpoly,tile$x[[j]],tile$y[[j]])
>      }
> 
>      curpoly = sprintf("%s %.6f %.6f))",curpoly,tile$x[[1]],tile$y[[1]])
> 
>      # this bit will find the original point that corresponds to the current polygon, along with its id
>      # and the SRID used for the point geometry (presumably this is the same for all points)...
>      # this will also filter out the extra polygons created for the
>      # four dummy points, as they will not return a result from this query:
> 
>      ipoint <- pg.spi.exec ( sprintf ( "select %3$s as id, intersection(''SRID=''||srid(%2$s)||'';%4$s'',''%5$s'') as polygon from %1$s where intersects(%2$s,''SRID=''||srid(%2$s)||'';%4$s'');", arg1, arg2, arg3, curpoly, buffer$ewkb[1] ) )
> 
>      if (length(ipoint) > 0)
>      {
>        poly[[p]] <- ipoint$polygon[1]
>        id[[p]] <- ipoint$id[1]
>        p = (p + 1)
>      }
>    }
> 
> Puneet Kishor wrote:
>> On Mar 5, 2012, at 9:44 PM, Derek Jones wrote:
>>> Hi all,
>>> 
>>> I have used an R solution that works well with the plsql to do this. Found here:
>>> 
>>> http://www.bostongis.com/PrinterFriendly.aspx?content_name=postgresql_plr_tut02
>>> 
>>> Needed some mods for my local solution, but helpful.
>>> 
>> Yes, I tried that, but as I noted in my earlier email, I got the following error
>> ERROR:  R interpreter expression evaluation error
>> DETAIL:  Error in pg.spi.exec(sprintf("select x(%2$s) as x, y(%2$s) as y from %1$s;",  :   error in SQL statement : syntax error at or near ";"
>> CONTEXT:  In R support function pg.spi.exec
>> In PL/R function voronoi