# R based Delaunay Triangulation Function for PostGIS using the deldir package

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 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[],
triangle\$y[],
triangle\$z[]
)
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.