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
_______________________________________________
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users
