Go to content Go to navigation and search

Home

Current PostGIS Blog Articles


Search

Browse

RSS / Atom

Email me

textpattern

Creative Commons License
All Blog Articles, Data Models and Free Source Code by Simon Greener, The SpatialDB Advisor is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

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[[1]],
  61.            triangle$y[[1]],
  62.            triangle$z[[1]]
  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.

Creative Commons License

post this at del.icio.uspost this at Diggpost this at Technoratipost this at Redditpost this at Farkpost this at Yahoo! my webpost this at Windows Livepost this at Google Bookmarkspost this to Twitter

Comment