[postgis-users] A function for using PL/R with deldir to produce voronoi polygons in PostGIS

[Mike Leahy]

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.
-------------- next part --------------
--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;