Voronoi Diagrams In PostGIS (with PL/R)
A Voronoi diagram represents proximity information about a set of points. The points on the Voronoi diagram are equidistant to two or more sites. Per Voronoi Diagrams and Delaunay Triangulation
Voronoi diagrams were first discussed by Peter Lejeune-Dirichlet in 1850. But it was more than a half of a century later in 1908 that these diagrams were written about in a paper by Voronoi, hence the name Voronoi Diagrams. The Voronoi cells/polygons are sometimes also called Thiessen Polytopes or Dirichlet Regions.
Ever wanted to create a Voronoi tessellation for your points? Well, so have I, for our Macrostrat points. Fortunately for us, others already figured out how to do this. I have dusted off the original code by Mike Leahy ever so slightly, and added a little bit more explanation, but it is good to go. It does require Postgres/PostGIS with PL/R. I am using Pg 9.0.x, PostGIS 1.5.3, and the latest PL/R as installed with help of MacPorts. Here is how to use the function once it is installed.
SELECT v.id, v.polygon
FROM r_voronoi('table', 'the_geom', 'id') As v;
or
SELECT v.id, v.polygon
FROM r_voronoi(
'(SELECT id, the_geom FROM table LIMIT 10) AS p',
'p.the_geom',
'p.id'
) As v;
The code for the function is below. All my contribution is under a CC0 License Waiver. I claim no responsibility for your successes or failures.
/*
This function uses the deldir library in R to generate voronoi
polygons for an input set of points in a PostGIS table.
Original function by Mike Leahy mgleahy at alumni.uwaterloo.ca
http://postgis.refractions.net/pipermail/postgis-users/2007-June/016102.html
Mar 6, 2012
Cleaned up formatting, and updated a teensy bit to modern times
Puneet Kishor punkish@eidesis.org
All my contribution is released under a CC0 License Waiver
Effectively in the Public Domain
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 r_voronoi('table','point-column','id-column');
Where:
'table' is the table or query (see below) containing the points to
be usedfor 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
Alternative usage: Instead of the name of a table, pass a query
that returns a result set. Make sure to us an ALIAS for the
result set.
SELECT * FROM r_voronoi(
'(SELECT id, the_geom FROM table LIMIT 10) AS p',
'p.the_geom',
'p.id'
)
*/
DROP FUNCTION r_voronoi(text, text, text);
DROP TYPE voronoi;
CREATE TYPE voronoi AS (id integer, polygon geometry);
CREATE OR REPLACE FUNCTION r_voronoi(text, text, text)
RETURNS SETOF voronoi AS '
R code below
library(deldir)
# select the point x/y coordinates into a data frame
points <- pg.spi.exec(
sprintf(
"SELECT ST_X(%2$s) AS x, ST_Y(%2$s) AS y FROM %1$s;",
arg1,
arg2
)
)
# calculate an approprate buffer distance (~10%):
buffer_distance = (
(
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_set <- pg.spi.exec(
sprintf(
"SELECT ST_Buffer(ST_Convexhull(ST_Union(%2$s)),%3$.6f) AS ewkb
FROM %1$s;",
arg1,
arg2,
buffer_distance
)
)
# 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
# construct the outgoing WKT now
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_set$ewkb[1]
)
)
if (length(ipoint) > 0) {
poly[[p]] <- ipoint$polygon[1]
id[[p]] <- ipoint$id[1]
p = (p + 1)
}
}
return(data.frame(id,poly))
Close off with SQL
' language 'plr';