 Home Current PostGIS Blog Articles

Search

Browse

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

Monday April 23 2012 at 00:18

KeywordsR postgis postgresql triangles delaunay XYZ triangulation
Summary

This article shows how to use the deldir R package within Postgresql to generate delaunay triangles that maintain the Z attribute on the vertices of the constructed triangles.

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:

1. DROP FUNCTION IF EXISTS r_delaunay(text);
2. CREATE OR REPLACE FUNCTION r_delaunay(p_query text)
3. RETURNS SETOF text
4. AS
5. '   library(deldir)
6.    # select the point x/y/z coordinates into a data frame
7.    points <- pg.spi.exec(p_query)
8.    # calculate an approprate buffer distance (~10%):
9.    buffer_distance = (
10.        (
11.            abs(max(points\$x) - min(points\$x)) +
12.            abs(max(points\$y) - min(points\$y))
13.        ) / 2
14.    ) * (0.10)
15.    # the following use of deldir uses high precision and digits to prevent
16.    # slivers between the output triangles, and uses a relatively large bounding
17.    # box with four dummy points included to ensure that points in the
18.    # peripheral areas of the dataset are appropriately enveloped by their
19.    # corresponding triangles:
20.    #    points\$x,
21.    #    points\$y,
22.    #    points\$z,
23.    voro = deldir(
24.        list( x=points\$x, y=points\$y, z=points\$z),
25.        digits=22,
26.        frac=0.00000000000000000000000001,
27.        list(ndx=2,ndy=2),
28.        rw=c(
29.            min(points\$x) - abs(min(points\$x) - max(points\$x)),
30.            max(points\$x) + abs(min(points\$x) - max(points\$x)),
31.            min(points\$y) - abs(min(points\$y) - max(points\$y)),
32.            max(points\$y) + abs(min(points\$y) - max(points\$y))
33.        )
34.    )
35.    # Do the triangulation
36.    #
37.    triangles = triang.list(voro)
38.    # Now create the output
39.    #
40.    poly = array()
41.    id = array()
42.    p = 1
43.    # construct the outgoing WKT now
44.    #
45.    for (i in 1:length(triangles)) {
46.        triangle = triangles[[i]]
47.        wktpoly = "POLYGON(("
48.        for (j in 1:length(triangle\$x)) {
49.             wktpoly = sprintf(
50.                "%s%.8f %.8f %.2f,",
51.                wktpoly,
52.                triangle\$x[[j]],
53.                triangle\$y[[j]],
54.                triangle\$z[[j]]
55.             )
56.        }
57.        wktpoly = sprintf(
58.            "%s%.8f %.8f %.2f))",
59.            wktpoly,
60.            triangle\$x[],
61.            triangle\$y[],
62.            triangle\$z[]
63.        )
64.        poly[[p]] <- wktpoly
65.        p = (p + 1)
66.    }
67.    return(data.frame(poly))
68. ' LANGUAGE 'plr';

Testing this we get:

1. SELECT 'SRID=4326;' ||
2. r_delaunay('
3. WITH xyz AS (
4. SELECT 403 as Z, ST_GeomFromText(''POINT (-90.2618913358086 33.2184179805337)'',4326) as geomXY UNION ALL
5. SELECT 332 as Z, ST_GeomFromText(''POINT (-90.2534717690058 33.2235960587115)'',4326) as geomXY UNION ALL
6. SELECT 393 as Z, ST_GeomFromText(''POINT (-90.2526790106307 33.2219983038141)'',4326) as geomXY UNION ALL
7. SELECT 349 as Z, ST_GeomFromText(''POINT (-90.2548360643114 33.2220152472292)'',4326) as geomXY UNION ALL
8. SELECT 375 as Z, ST_GeomFromText(''POINT (-90.2569931215599 33.2220321532301)'',4326) as geomXY UNION ALL
9. SELECT 329 as Z, ST_GeomFromText(''POINT (-90.2591702280608 33.2202360124111)'',4326) as geomXY UNION ALL
10. SELECT 242 as Z, ST_GeomFromText(''POINT (-90.2570132117573 33.2202191464654)'',4326) as geomXY UNION ALL
11. SELECT 311 as Z, ST_GeomFromText(''POINT (-90.2548561990173 33.2202022431079)'',4326) as geomXY UNION ALL
12. SELECT 296 as Z, ST_GeomFromText(''POINT (-90.2526991898448 33.2201853023387)'',4326) as geomXY UNION ALL
13. SELECT 375 as Z, ST_GeomFromText(''POINT (-90.2527193654874 33.2183723003489)'',4326) as geomXY UNION ALL
14. SELECT 331 as Z, ST_GeomFromText(''POINT (-90.2548763301557 33.2183892384722)'',4326) as geomXY UNION ALL
15. SELECT 246 as Z, ST_GeomFromText(''POINT (-90.2570332983911 33.2184061391864)'',4326) as geomXY UNION ALL
16. SELECT 292 as Z, ST_GeomFromText(''POINT (-90.2591902701899 33.2184230024914)'',4326) as geomXY UNION ALL
17. SELECT 228 as Z, ST_GeomFromText(''POINT (-90.2602038624534 33.2165408049195)'',4326) as geomXY UNION ALL
18. SELECT 295 as Z, ST_GeomFromText(''POINT (-90.2580468665807 33.2165321889039)'',4326) as geomXY UNION ALL
19. SELECT 305 as Z, ST_GeomFromText(''POINT (-90.2558898738654 33.2165235354584)'',4326) as geomXY UNION ALL
20. SELECT 396 as Z, ST_GeomFromText(''POINT (-90.2537328843112 33.2165148445831)'',4326) as geomXY UNION ALL
21. SELECT 329 as Z, ST_GeomFromText(''POINT (-90.2537432444696 33.2147017821533)'',4326) as geomXY UNION ALL
22. SELECT 213 as Z, ST_GeomFromText(''POINT (-90.255900189497 33.2147104705514)'',4326) as geomXY UNION ALL
23. SELECT 218 as Z, ST_GeomFromText(''POINT (-90.2580571376854 33.2147191215224)'',4326) as geomXY UNION ALL
24. SELECT 226 as Z, ST_GeomFromText(''POINT (-90.2602140890308 33.2147277350661)'',4326) as geomXY UNION ALL
25. SELECT 310 as Z, ST_GeomFromText(''POINT (-90.2591415442112 33.2133037851107)'',4326) as geomXY UNION ALL
26. SELECT 228 as Z, ST_GeomFromText(''POINT (-90.254879330931 33.2133111535788)'',4326) as geomXY
27. )
28. SELECT ST_X(A.geomXY) as x,ST_Y(A.geomXY) as y,A.z
29.   FROM xyz a;') as triangles;
30. -- Results
31. "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))"
32. "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))"
33. "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))"
34. "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))"
35. "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))"
36. "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))"
37. "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))"
38. "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))"
39. "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))"
40. "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))"
41. "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))"
42. "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))"
43. "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))"
44. "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))"
45. "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))"
46. "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))"
47. "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))"
48. "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))"
49. "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))"
50. "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))"
51. "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))"
52. "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))"
53. "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))"
54. "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))"
55. "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))"
56. "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))"
57. "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))"
58. "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))"
59. "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))"
60. "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))"
61. "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))"
62. "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))"
63. "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))"
64. "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))"
65. "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.          ### Comment

 Name Remember E-mail Message Textile Help 