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;


    SELECT v.id, v.polygon
    FROM r_voronoi(
        '(SELECT id, the_geom FROM table LIMIT 10) AS p',
    ) 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

      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

          R-2.5.0 with deldir-0.0-5 installed
          PostgreSQL-8.2.x with PostGIS-1.x and PL/R- installed

      Usage: SELECT * FROM r_voronoi('table','point-column','id-column');

      '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',
    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


        # select the point x/y coordinates into a data frame
        points <- pg.spi.exec(
                "SELECT ST_X(%2$s) AS x, ST_Y(%2$s) AS y FROM %1$s;",

        # 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(
                "SELECT ST_Buffer(ST_Convexhull(ST_Union(%2$s)),%3$.6f) AS ewkb 
                FROM %1$s;",

        # 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(
                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 = sprintf(
                "%s %.6f %.6f))",

            # 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(
                    "SELECT %3$s AS id, 
                       AS polygon 
                    FROM %1$s 
                    WHERE intersects(%2$s,''SRID=''||srid(%2$s)||'';%4$s'');",
            if (length(ipoint) > 0) {
                poly[[p]] <- ipoint$polygon[1]
                id[[p]] <- ipoint$id[1]
                p = (p + 1)

Close off with SQL

    ' language 'plr';