[postgis-users] Voronoi tessellation

[Derek Jones]

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
> 
> 
>