# Voronoi Diagrams In PostGIS (with PL/R)

Invalid Date

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[],
tile\$y[]
)

# 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
)
)
if (length(ipoint) > 0) {
poly[[p]] <- ipoint\$polygon
id[[p]] <- ipoint\$id
p = (p + 1)
}
}
return(data.frame(id,poly))
``````

Close off with SQL

``````
' language 'plr';
``````