# [postgis-users] Voronoi tessellation

Derek Jones scunacc at yahoo.com
Mon Mar 5 20:55:15 PST 2012

```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[],tile\$y[])

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

if (length(ipoint) > 0)
{
poly[[p]] <- ipoint\$polygon
id[[p]] <- ipoint\$id
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
>
>
>

```

More information about the postgis-users mailing list