All,
Attached is a PL/R-based function that I put together based on Jan
Hartmann's help in the PotGIS list about a year ago for generating
voronoi polygons using the deldir library. Now that I've started with
PL/R, I was able capture this in a single function that will take in the
points from a table or query, produce the voronoi polygons, and return a
list of the polygon geometries along with the id's of their originating
points.
I'm not sure if the approach I used in this example is ideal; it hasn't
been tested on large datasets yet, and I'm sure someone can probably
find a different way to do this (there is also the tripack library in R,
which I haven't tried yet). However, it seems to work fine from what I
can tell. I just thought this might serve as an example of how PL/R can
be handy alongside PostGIS.
Regards,
Mike
P.S.: I did notice the same problem as Regina - once I started trying to
execute this script directly from the text file, I had to make sure it
was in UNIX format for PL/R to parse it properly...even though
everything I'm doing is on Windows at the moment.
--This function uses the deldir library in R to generate voronoi polygons for
an input set of points in a PostGIS table.
--Requirements:
-- R-2.5.0 with deldir-0.0-5 installed
-- PostgreSQL-8.2.x with PostGIS-1.x and PL/R-8.2.0.4 installed
--Usage: select * from voronoi('table','point-column','id-column');
--Where:
-- 'table' is the table or query containing the points to be used for
generating the voronoi polygons,
-- 'point-column' is a single 'POINT' PostGIS geometry type (each point
must be unique)
-- 'id-column' is a unique identifying integer for each of the originating
points (e.g., 'gid')
--Output: returns a recordset of the custom type 'voronoi', which contains the
id of the
-- originating point, and a polygon geometry
drop function voronoi(text,text,text);
drop type voronoi;
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)
}
}
return(data.frame(id,poly))
' language 'plr';
--An example of how this function can be used:
--create a table with some random (but unique) points:
create table geo (id serial primary key);
select addgeometrycolumn('','geo','pt',4326,'POINT',2);
insert into geo (pt) values ('SRID=4326;POINT(1 -4)');
insert into geo (pt) values ('SRID=4326;POINT(3.5 5)');
insert into geo (pt) values ('SRID=4326;POINT(-2 2)');
insert into geo (pt) values ('SRID=4326;POINT(6 -9.25)');
insert into geo (pt) values ('SRID=4326;POINT(0 -3)');
insert into geo (pt) values ('SRID=4326;POINT(7 2)');
--get the voronoi EWKB geometry and originating point IDs
select * from voronoi('geo','pt','id');
--same as above, but instead of a direct table as input, a subquery with
aliases can be used:
select * from voronoi('(select id as id_alias, pt as pt_alias from geo) as
x','x.pt_alias','x.id_alias');
--same as above, but viewing output geometries as EWKT
select id, asewkt(polygon) from voronoi('geo','pt','id');
--join the voronoi output to the originating table of points:
select geo.*, v.polygon as voronoi from geo left join (select * from
voronoi('geo','pt','id')) as v using (id);
--cleanup the sample table:
select dropgeometrycolumn('','geo','pt');
drop table geo;
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
