Constructing Voronoi diagram in PostGIS

Question

I'm trying to construct voronoi diagram from grid of points using modified code from here. This is SQL query after my modifications:

DROP TABLE IF EXISTS example.voronoi;
WITH 
    -- Sample set of points to work with
    Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM example."MeshPoints2d"),
    -- Build edges and circumscribe points to generate a centroid
    Edges AS (
    SELECT id,
        UNNEST(ARRAY['e1','e2','e3']) EdgeName,
        UNNEST(ARRAY[
            ST_MakeLine(p1,p2) ,
            ST_MakeLine(p2,p3) ,
            ST_MakeLine(p3,p1)]) Edge,
        ST_Centroid(ST_ConvexHull(ST_Union(-- Done this way due to issues I had with LineToCurve
            ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p1,p2),ST_MakeLine(p2,p3)))),'LINE','CIRCULAR'),15),
            ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p2,p3),ST_MakeLine(p3,p1)))),'LINE','CIRCULAR'),15)
        ))) ct      
    FROM    (
        -- Decompose to points
        SELECT id,
            ST_PointN(g,1) p1,
            ST_PointN(g,2) p2,
            ST_PointN(g,3) p3
        FROM    (
            SELECT (gd).Path id, ST_ExteriorRing((gd).geom) g -- ID andmake triangle a linestring
            FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a -- Get Delaunay Triangles
            )b
        ) c
    )
SELECT ST_SetSRID((ST_Dump(ST_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, ST_ExteriorRing(ST_ConvexHull(v)))))))).geom, 2180)
INTO example.voronoi
FROM (
    SELECT  -- Create voronoi edges and reduce to a multilinestring
        ST_LineMerge(ST_Union(ST_MakeLine(
        x.ct,
        CASE 
        WHEN y.id IS NULL THEN
            CASE WHEN ST_Within(
                x.ct,
                (SELECT ST_ConvexHull(geom) FROM sample)) THEN -- Don't draw lines back towards the original set
                -- Project line out twice the distance from convex hull
                ST_MakePoint(ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * 2),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 2))
            END
        ELSE 
            y.ct
        END
        ))) v
    FROM    Edges x 
        LEFT OUTER JOIN -- Self Join based on edges
        Edges y ON x.id <> y.id AND ST_Equals(x.edge,y.edge)
    ) z

Below - result of my query.

As you can see I get "almost" correct voronoi diagram, except highlighted points which don't have unique polygon. Below is what QGIS algorithm produce and what I would like to obtain from query. Any sugestions where is the problem with my code?

Perhaps you could also compare the result of the SpatiaLite function "VoronojDiagram" gaia-gis.it/gaia-sins/spatialite-sql-latest.html and have a look at the source code in gaia-gis.it/fossil/libspatialite/index. — user30184, Dec 1, 2015 at 18:05
Nice question, I have been looking at the same question you reference, with a view to speeding it up, but keep running out of time. I wasn't aware of this issue with the outer points. — John Powell, Dec 1, 2015 at 19:23
For what it's worth we have ST_Voronoi coming in PostGIS 2.3, Dan Baston will soon be checking in code for it - trac.osgeo.org/postgis/ticket/2259 looks pretty much done just needs to be pulled in. Would be great to have folks testing — LR1234567, Dec 1, 2015 at 21:03
Can you post up the set of points you are using? I would mind doing a bit of testing on this myself — MickyT, Dec 1, 2015 at 21:05

John Powell · Accepted Answer · 2016-12-01 17:13:23Z

Following a suggestion by @LR1234567 to try out the new ST_Voronoi functionality that has been developed by @dbaston, @MickyT 's original amazing answer (as stated in OP's question) and using the original data can now be simplified to:

WITH voronoi (vor) AS 
     (SELECT ST_Dump(ST_Voronoi(ST_Collect(geom))) FROM meshpoints)
SELECT (vor).path, (vor).geom FROM voronoi;

which results in this output, identical to the OP's question.

However, this suffers from the same issue, now fixed in MickyT's answer, that points on the concave hull do not get an enclosing polygon (as would be expected from the algorithm). I fixed this issue with a query with the following series of steps.

Calculate the concave hull of the input points -- points on the concave hull are those that have unbounded polygons in the output Voronoi diagram.
Find the original points on the concave hull (yellow points in diagram 2 below).
Buffer the concave hull (the buffer distance is arbitrary and could perhaps be found more optimally from input data?).
Find the closest points on the buffer of concave hull nearest to points in step 2. These are show as green in the diagram below.
Add these points to original data set
Calculate the Voronoi diagram of this combined data set. As can be seen in the 3rd diagram, the points on the hull now have enclosing polygons.

Diagram 2 showing points on concave hull (yellow) and closest points to buffer on hull (green). Diagram 2.

The query could obviously be simplified/compressed, but I left it is this form as a series of CTEs, as it is easier to follow the steps sequentially that way. This query runs on the original data set in milliseconds (11ms average on a dev server) whereas MickyT's answer using ST_Delauney runs in 4800ms on the same server. DBaston claims that another order of magnitude speed enhancement can be gained from building against the current trunk of GEOS, 3.6dev, due to enhancements in triangulation routines.

WITH 
  conv_hull(geom) AS 
        (SELECT ST_Concavehull(ST_Union(geom), 1) FROM meshpoints), 
  edge_points(points) AS 
        (SELECT mp.geom FROM meshpoints mp, conv_hull ch 
        WHERE ST_Touches(ch.geom, mp.geom)), 
  buffered_points(geom) AS
        (SELECT ST_Buffer(geom, 100) as geom FROM conv_hull),
  closest_points(points) AS
        (SELECT 
              ST_Closestpoint(
                   ST_Exteriorring(bp.geom), ep.points) as points,
             ep.points as epoints 
         FROM buffered_points bp, edge_points ep),
  combined_points(points) AS
        (SELECT points FROM closest_points 
        UNION SELECT geom FROM meshpoints),
  voronoi (vor) AS 
       (SELECT 
            ST_Dump(
                  ST_Voronoi(
                    ST_Collect(points))) as geom 
        FROM combined_points)
 SELECT 
     (vor).path[1] as id, 
     (vor).geom 
 FROM voronoi;

Diagram 3 showing all points now enclosed in a polygon

Note: Currently ST_Voronoi involves building Postgis from source (version 2.3, or trunk) and linking against GEOS 3.5 or above.

Edit: I've just looked at Postgis 2.3 as it is installed on Amazon Web Services, and it seems the function name is now ST_VoronoiPolygons.

No doubt this query/algorithm could be improved. Suggestions welcome.

@dbaston. Wondering if you have any comments on this approach? — John Powell, Dec 19, 2015 at 10:03
well, all of the points do get an enclosing polygon, it's just that it's disproportionately large. If and how these should be made smaller is pretty subjective, and without knowing exactly what is desired for the outer polygons it's hard to know what the "best" way is. Yours seems a fine method to me. I've used a less sophisticated method that's similar in spirit to yours, dropping extra points along a buffered convex hull boundary at a fixed spacing determined by the avg point density. — dbaston, Dec 23, 2015 at 18:55
@dbaston. Thanks, just making sure I wasn't missing something obvious. An algorithm for shrinking the outer polygons to something more in line with the size of the inner ones (in my case postcode areas) is something I will have to think some more about. — John Powell, Dec 24, 2015 at 11:11
@John Barça Thanks for another great solution. Calculations speed is more than satisfying with this approach. Unfortunately I would like to use this algorithm inside my QGIS plugin, and it have to work with PostGIS 2.1+ out of the box. But surely I will use this solution after PostGIS 2.3 official release. Anyway thank you for such comprehensive answers. :) — DamnBack, Dec 28, 2015 at 18:20
@DamnBack. You are very welcome. I needed this for work and your question actually helped me a lot, as I had no idea about ST_Voronoi coming out, and the older solutions are much slower (as you have noticed). It was fun figuring it out too :-) — John Powell, Dec 28, 2015 at 20:03

MickyT · Accepted Answer · 2015-12-02 02:57:24Z

While this fixes the immediate issue with the query for the data in question, I am not happy with it as a solution for general usage and I will revisit this and the previous answer when I can.

The issue was that the original query was using a convex hull on the voronoi edges to determine the exterior edge for the voronoi polygon. This meant that some of the voronoi edges weren't reaching the exterior when they should have been. I looked at using the concave hull functionality, but it didn't really work as I wanted.
As a quick fix I have changed convex hull to be built around the closed set of voronoi edges plus a buffer around the original edges. The voronoi edges that don't close are extended out a large distance to try and make sure that they intersect the exterior and are use to build polygons. You may want to play around with the buffer parameters.

WITH 
    -- Sample set of points to work with
    Sample AS (SELECT ST_SetSRID(ST_Union(geom), 0) geom FROM MeshPoints2d),
    -- Build edges and circumscribe points to generate a centroid
    Edges AS (
    SELECT id,
        UNNEST(ARRAY['e1','e2','e3']) EdgeName,
        UNNEST(ARRAY[
            ST_MakeLine(p1,p2) ,
            ST_MakeLine(p2,p3) ,
            ST_MakeLine(p3,p1)]) Edge,
        ST_Centroid(ST_ConvexHull(ST_Union(-- Done this way due to issues I had with LineToCurve
            ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p1,p2),ST_MakeLine(p2,p3)))),'LINE','CIRCULAR'),15),
            ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p2,p3),ST_MakeLine(p3,p1)))),'LINE','CIRCULAR'),15)
        ))) ct      
    FROM    (
        -- Decompose to points
        SELECT id,
            ST_PointN(g,1) p1,
            ST_PointN(g,2) p2,
            ST_PointN(g,3) p3
        FROM    (
            SELECT (gd).Path id, ST_ExteriorRing((gd).geom) g -- ID andmake triangle a linestring
            FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a -- Get Delaunay Triangles
            )b
        ) c
    )
SELECT ST_SetSRID((ST_Dump(ST_Polygonize(ST_Node(ST_LineMerge(ST_Union(v, (SELECT ST_ExteriorRing(ST_ConvexHull(ST_Union(ST_Union(ST_Buffer(edge,20),ct)))) FROM Edges))))))).geom, 2180) geom
INTO voronoi
FROM (
    SELECT  -- Create voronoi edges and reduce to a multilinestring
        ST_LineMerge(ST_Union(ST_MakeLine(
        x.ct,
        CASE 
        WHEN y.id IS NULL THEN
            CASE WHEN ST_Within(
                x.ct,
                (SELECT ST_ConvexHull(geom) FROM sample)) THEN -- Don't draw lines back towards the original set
                -- Project line out twice the distance from convex hull
                ST_MakePoint(ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * 200),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 200))
            END
        ELSE 
            y.ct
        END
        ))) v
    FROM    Edges x 
        LEFT OUTER JOIN -- Self Join based on edges
        Edges y ON x.id <> y.id AND ST_Equals(x.edge,y.edge)
    ) z;

Thanks for explanation and quick fix for the problem! It works with my data (a little bit slower due to ST_Union(ST_Buffer(geom))), but i will continue testing with other sets of points. Meanwhile I will waiting for as you said - more general solution. :) — DamnBack, Dec 2, 2015 at 9:59

LR1234567 · Accepted Answer · 2015-12-02 05:09:54Z

3

If you have access to PostGIS 2.3, give the new ST_Voronoi function a try, recently committed:

http://postgis.net/docs/manual-dev/ST_Voronoi.html

There are precompiled builds for windows - http://postgis.net/windows_downloads/

answered Dec 2, 2015 at 5:09

LR1234567

10.3k1 gold badge22 silver badges28 bronze badges

Thanks for info that there is upcomming built-in ST_Voronoi function - I will check it out. Unfortunately I need solution that works on PostGIS 2.1+ versions, so @MickyT query is the closest to my needs at the moment.
– DamnBack
Dec 2, 2015 at 10:12
@LR1234567. Does this require any particular version of GEOS. I have time tomorrow to test 2.3 and ST_Voronoi.
– John Powell
Dec 3, 2015 at 14:37
Requires GEOS 3.5
– LR1234567
Dec 3, 2015 at 16:19

Add a comment |

Stack Exchange Network

Constructing Voronoi diagram in PostGIS

3 Answers 3

Your Answer

Not the answer you're looking for? Browse other questions tagged
postgis
sql
voronoi-thiessen
or ask your own question.

Linked

Hot Network Questions

Constructing Voronoi diagram in PostGIS

3 Answers 3

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged postgissqlvoronoi-thiessen or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
postgis
sql
voronoi-thiessen
or ask your own question.