R based Delaunay Triangulation Function for PostGIS using the deldir package

31 Jul 2020 By Simon Greener

Introduction

Recently I had need of a function in PostgreSQL/PostGIS that generated a set of Delaunay triangles from a set of XYZ points.

I happened upon Regina Obe’s Voronoi implementation that is based on deldir.

Taking this as a basis I installed R and PL/R, and then used Mike Leahy’s original Voronoi function (see PostGIS users newsgroup) as the basis for a new pg/PLSQL function called r_delaunay.

Initial testing showed that it certainly generated Delaunay triangles except that it only produced 2D triangles because the version of deldir I was using did not maintain the Z ordinate in the process.

I contacted the writer of the deldir R package, University of Auckland (NZ) lecturer Dr Rolf Turner about the problem. Rolf changed the inputs to the base deldir function and released the changes as Version0.0-18). These changes are:

deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, sort=TRUE, plotit=FALSE, digits=6, z=NULL, zdum=NULL, …)

Now the input coordinates of the point set being triangulated can be described by:
1) Two arguments x and y (vectors) or
2) A single argument x (a list with components x and y, and possibly z (or any associated auxilliary values or weights).

Using this new version of deldir, with Regina’s voronoi() function as a template, I produced the following:

 DROP FUNCTION IF EXISTS r_delaunay(text);
 CREATE OR REPLACE FUNCTION r_delaunay(p_query text)
 RETURNS SETOF text
 AS
 '   library(deldir)
    # select the point x/y/z coordinates into a data frame
    points <- pg.spi.exec(p_query)
    # 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)
    # the following use of deldir uses high precision and digits to prevent
    # slivers between the output triangles, 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 triangles:
    #    points$x,
    #    points$y,
    #    points$z,
    voro = deldir(
        list( x=points$x, y=points$y, z=points$z),
        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))
        )
    )
    # Do the triangulation
    #
    triangles = triang.list(voro)
    # Now create the output
    #
    poly = array()
    id = array()
    p = 1
    # construct the outgoing WKT now
    #
    for (i in 1:length(triangles)) {
        triangle = triangles[[i]]
        wktpoly = "POLYGON(("
        for (j in 1:length(triangle$x)) {
             wktpoly = sprintf(
                "%s%.8f %.8f %.2f,",
                wktpoly,
                triangle$x[[j]],
                triangle$y[[j]],
                triangle$z[[j]]
             )
        }
        wktpoly = sprintf(
            "%s%.8f %.8f %.2f))",
            wktpoly,
            triangle$x[[1]],
            triangle$y[[1]],
            triangle$z[[1]]
        )
        poly[[p]] <- wktpoly
        p = (p + 1)
    }
    return(data.frame(poly))
 ' LANGUAGE 'plr';

Testing this we get:

 SELECT 'SRID=4326;' ||
 r_delaunay('
 WITH xyz AS (
 SELECT 403 as Z, ST_GeomFromText(''POINT (-90.2618913358086 33.2184179805337)'',4326) as geomXY UNION ALL
 SELECT 332 as Z, ST_GeomFromText(''POINT (-90.2534717690058 33.2235960587115)'',4326) as geomXY UNION ALL
 SELECT 393 as Z, ST_GeomFromText(''POINT (-90.2526790106307 33.2219983038141)'',4326) as geomXY UNION ALL
 SELECT 349 as Z, ST_GeomFromText(''POINT (-90.2548360643114 33.2220152472292)'',4326) as geomXY UNION ALL
 SELECT 375 as Z, ST_GeomFromText(''POINT (-90.2569931215599 33.2220321532301)'',4326) as geomXY UNION ALL
 SELECT 329 as Z, ST_GeomFromText(''POINT (-90.2591702280608 33.2202360124111)'',4326) as geomXY UNION ALL
 SELECT 242 as Z, ST_GeomFromText(''POINT (-90.2570132117573 33.2202191464654)'',4326) as geomXY UNION ALL
 SELECT 311 as Z, ST_GeomFromText(''POINT (-90.2548561990173 33.2202022431079)'',4326) as geomXY UNION ALL
 SELECT 296 as Z, ST_GeomFromText(''POINT (-90.2526991898448 33.2201853023387)'',4326) as geomXY UNION ALL
 SELECT 375 as Z, ST_GeomFromText(''POINT (-90.2527193654874 33.2183723003489)'',4326) as geomXY UNION ALL
 SELECT 331 as Z, ST_GeomFromText(''POINT (-90.2548763301557 33.2183892384722)'',4326) as geomXY UNION ALL
 SELECT 246 as Z, ST_GeomFromText(''POINT (-90.2570332983911 33.2184061391864)'',4326) as geomXY UNION ALL
 SELECT 292 as Z, ST_GeomFromText(''POINT (-90.2591902701899 33.2184230024914)'',4326) as geomXY UNION ALL
 SELECT 228 as Z, ST_GeomFromText(''POINT (-90.2602038624534 33.2165408049195)'',4326) as geomXY UNION ALL
 SELECT 295 as Z, ST_GeomFromText(''POINT (-90.2580468665807 33.2165321889039)'',4326) as geomXY UNION ALL
 SELECT 305 as Z, ST_GeomFromText(''POINT (-90.2558898738654 33.2165235354584)'',4326) as geomXY UNION ALL
 SELECT 396 as Z, ST_GeomFromText(''POINT (-90.2537328843112 33.2165148445831)'',4326) as geomXY UNION ALL
 SELECT 329 as Z, ST_GeomFromText(''POINT (-90.2537432444696 33.2147017821533)'',4326) as geomXY UNION ALL
 SELECT 213 as Z, ST_GeomFromText(''POINT (-90.255900189497 33.2147104705514)'',4326) as geomXY UNION ALL
 SELECT 218 as Z, ST_GeomFromText(''POINT (-90.2580571376854 33.2147191215224)'',4326) as geomXY UNION ALL
 SELECT 226 as Z, ST_GeomFromText(''POINT (-90.2602140890308 33.2147277350661)'',4326) as geomXY UNION ALL
 SELECT 310 as Z, ST_GeomFromText(''POINT (-90.2591415442112 33.2133037851107)'',4326) as geomXY UNION ALL
 SELECT 228 as Z, ST_GeomFromText(''POINT (-90.254879330931 33.2133111535788)'',4326) as geomXY
 )
 SELECT ST_X(A.geomXY) as x,ST_Y(A.geomXY) as y,A.z
   FROM xyz a;') as triangles;
 -- Results
 "SRID=4326;POLYGON((-90.26189134 33.21841798 403.00,-90.25917023 33.22023601 329.00,-90.25699312 33.22203215 375.00,-90.26189134 33.21841798 403.00))"
 "SRID=4326;POLYGON((-90.26189134 33.21841798 403.00,-90.25919027 33.21842300 292.00,-90.25917023 33.22023601 329.00,-90.26189134 33.21841798 403.00))"
 "SRID=4326;POLYGON((-90.26189134 33.21841798 403.00,-90.26020386 33.21654080 228.00,-90.25919027 33.21842300 292.00,-90.26189134 33.21841798 403.00))"
 "SRID=4326;POLYGON((-90.26189134 33.21841798 403.00,-90.26021409 33.21472774 226.00,-90.26020386 33.21654080 228.00,-90.26189134 33.21841798 403.00))"
 "SRID=4326;POLYGON((-90.25347177 33.22359606 332.00,-90.25483606 33.22201525 349.00,-90.25267901 33.22199830 393.00,-90.25347177 33.22359606 332.00))"
 "SRID=4326;POLYGON((-90.25347177 33.22359606 332.00,-90.25699312 33.22203215 375.00,-90.25483606 33.22201525 349.00,-90.25347177 33.22359606 332.00))"
 "SRID=4326;POLYGON((-90.25267901 33.22199830 393.00,-90.25483606 33.22201525 349.00,-90.25269919 33.22018530 296.00,-90.25267901 33.22199830 393.00))"
 "SRID=4326;POLYGON((-90.25267901 33.22199830 393.00,-90.25269919 33.22018530 296.00,-90.25271937 33.21837230 375.00,-90.25267901 33.22199830 393.00))"
 "SRID=4326;POLYGON((-90.25483606 33.22201525 349.00,-90.25699312 33.22203215 375.00,-90.25485620 33.22020224 311.00,-90.25483606 33.22201525 349.00))"
 "SRID=4326;POLYGON((-90.25483606 33.22201525 349.00,-90.25485620 33.22020224 311.00,-90.25269919 33.22018530 296.00,-90.25483606 33.22201525 349.00))"
 "SRID=4326;POLYGON((-90.25699312 33.22203215 375.00,-90.25917023 33.22023601 329.00,-90.25701321 33.22021915 242.00,-90.25699312 33.22203215 375.00))"
 "SRID=4326;POLYGON((-90.25699312 33.22203215 375.00,-90.25701321 33.22021915 242.00,-90.25485620 33.22020224 311.00,-90.25699312 33.22203215 375.00))"
 "SRID=4326;POLYGON((-90.25917023 33.22023601 329.00,-90.25703330 33.21840614 246.00,-90.25701321 33.22021915 242.00,-90.25917023 33.22023601 329.00))"
 "SRID=4326;POLYGON((-90.25917023 33.22023601 329.00,-90.25919027 33.21842300 292.00,-90.25703330 33.21840614 246.00,-90.25917023 33.22023601 329.00))"
 "SRID=4326;POLYGON((-90.25701321 33.22021915 242.00,-90.25487633 33.21838924 331.00,-90.25485620 33.22020224 311.00,-90.25701321 33.22021915 242.00))"
 "SRID=4326;POLYGON((-90.25701321 33.22021915 242.00,-90.25703330 33.21840614 246.00,-90.25487633 33.21838924 331.00,-90.25701321 33.22021915 242.00))"
 "SRID=4326;POLYGON((-90.25485620 33.22020224 311.00,-90.25271937 33.21837230 375.00,-90.25269919 33.22018530 296.00,-90.25485620 33.22020224 311.00))"
 "SRID=4326;POLYGON((-90.25485620 33.22020224 311.00,-90.25487633 33.21838924 331.00,-90.25271937 33.21837230 375.00,-90.25485620 33.22020224 311.00))"
 "SRID=4326;POLYGON((-90.25271937 33.21837230 375.00,-90.25487633 33.21838924 331.00,-90.25373288 33.21651484 396.00,-90.25271937 33.21837230 375.00))"
 "SRID=4326;POLYGON((-90.25271937 33.21837230 375.00,-90.25373288 33.21651484 396.00,-90.25374324 33.21470178 329.00,-90.25271937 33.21837230 375.00))"
 "SRID=4326;POLYGON((-90.25487633 33.21838924 331.00,-90.25703330 33.21840614 246.00,-90.25588987 33.21652354 305.00,-90.25487633 33.21838924 331.00))"
 "SRID=4326;POLYGON((-90.25487633 33.21838924 331.00,-90.25588987 33.21652354 305.00,-90.25373288 33.21651484 396.00,-90.25487633 33.21838924 331.00))"
 "SRID=4326;POLYGON((-90.25703330 33.21840614 246.00,-90.25919027 33.21842300 292.00,-90.25804687 33.21653219 295.00,-90.25703330 33.21840614 246.00))"
 "SRID=4326;POLYGON((-90.25703330 33.21840614 246.00,-90.25804687 33.21653219 295.00,-90.25588987 33.21652354 305.00,-90.25703330 33.21840614 246.00))"
 "SRID=4326;POLYGON((-90.25919027 33.21842300 292.00,-90.26020386 33.21654080 228.00,-90.25804687 33.21653219 295.00,-90.25919027 33.21842300 292.00))"
 "SRID=4326;POLYGON((-90.26020386 33.21654080 228.00,-90.25805714 33.21471912 218.00,-90.25804687 33.21653219 295.00,-90.26020386 33.21654080 228.00))"
 "SRID=4326;POLYGON((-90.26020386 33.21654080 228.00,-90.26021409 33.21472774 226.00,-90.25805714 33.21471912 218.00,-90.26020386 33.21654080 228.00))"
 "SRID=4326;POLYGON((-90.25804687 33.21653219 295.00,-90.25590019 33.21471047 213.00,-90.25588987 33.21652354 305.00,-90.25804687 33.21653219 295.00))"
 "SRID=4326;POLYGON((-90.25804687 33.21653219 295.00,-90.25805714 33.21471912 218.00,-90.25590019 33.21471047 213.00,-90.25804687 33.21653219 295.00))"
 "SRID=4326;POLYGON((-90.25588987 33.21652354 305.00,-90.25374324 33.21470178 329.00,-90.25373288 33.21651484 396.00,-90.25588987 33.21652354 305.00))"
 "SRID=4326;POLYGON((-90.25588987 33.21652354 305.00,-90.25590019 33.21471047 213.00,-90.25374324 33.21470178 329.00,-90.25588987 33.21652354 305.00))"
 "SRID=4326;POLYGON((-90.25374324 33.21470178 329.00,-90.25590019 33.21471047 213.00,-90.25487933 33.21331115 228.00,-90.25374324 33.21470178 329.00))"
 "SRID=4326;POLYGON((-90.25590019 33.21471047 213.00,-90.25805714 33.21471912 218.00,-90.25487933 33.21331115 228.00,-90.25590019 33.21471047 213.00))"
 "SRID=4326;POLYGON((-90.25805714 33.21471912 218.00,-90.26021409 33.21472774 226.00,-90.25914154 33.21330379 310.00,-90.25805714 33.21471912 218.00))"
 "SRID=4326;POLYGON((-90.25805714 33.21471912 218.00,-90.25914154 33.21330379 310.00,-90.25487933 33.21331115 228.00,-90.25805714 33.21471912 218.00))"

Of course, there are many ways to pass the XYZ values into an r_delaunay function: the use of a string containing a valid SQL query that returns a set of rows of X, Y and Z values is but one.

I hope this is of help to someone.

Leave a Reply

Your email address will not be published. Required fields are marked *