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