Skip to content

Instantly share code, notes, and snippets.

@djq djq/voronoi.sql
Last active Feb 15, 2017

Embed
What would you like to do?
Updated SQL function to create Voronoi diagrams to work with PostGIS 2.0. Original code here: http://punkish.org/Voronoi-Diagrams-In-PostGIS
/*
Feb 5, 2013
Original source from here: http://punkish.org/Voronoi-Diagrams-In-PostGIS
Updated old PostGIS function names to work with PostGIS2.0 (buffer, intersect, intersection and srid)
-djq
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 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
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 TYPE voronoi;
DROP FUNCTION r_voronoi(text, text, text);
CREATE TYPE voronoi AS (id integer, polygon geometry);
CREATE OR REPLACE FUNCTION r_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 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, st_intersection(''SRID=''||st_srid(%2$s)||'';%4$s''::text,''%5$s'') AS polygon FROM %1$s WHERE st_intersects(%2$s::text,''SRID=''||st_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))
' language 'plr';
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.