[postgis-users] Voronoi / Dalaunay function with R

[Jan Hartmann]

As an alternative to Mike Leahy's and Mark Fenber's Voronoi postings:

There is an easy way to generate Voronoi polygons for PostGIS using R 
(wwww.r-project.org). In its simplest form you generate a text file with 
x-y coordinates from PostGIS, read these into R, compute the Voronoi 
diagram, and write the resulting polygons out as a text file with SQL 
insert commands:

- write out x-y-values to "points.txt"
   psql <database>
   \o points.txt
   \t
   select x,y from table
   \o

- The Voronoi package "deldir"
   (http://cran.r-project.org/src/contrib/Descriptions/deldir.html) is
   not part of the standard R distribution, so you have to install it
   with:
   "R CMD INSTALL deldir" (for Windows, install from the "packages" menu)

- R
   library(deldir)
   points = scan(file="points.txt",what=list(x=0.0,y=0.0),sep="|")
   voro = deldir(points$x,points$y)   # generate voronoi edges
   tiles = tile.list(voro)            # combine edges into polygons
   sink("voronoi.sql")                # redirect output to file
   for (i in 1:length(tiles)) {       # write out polygons
     tile = tiles[[i]]
     cat("insert into mytable (the_geom) values(geomfromtext('POLYGON((")
     for (j in 1:length(tile$x)) {
        cat (tile$x[[j]],' ',tile$y[[j]],",")
     }
     cat (tile$x[[1]],' ',tile$y[[1]]) #close polygon
     cat ("))',4326));\n")             # add SRID and newline
   }
   sink()                              # output to terminal

- In psql create a table "mytable" with a geometrycolumn "the_geom" and 
insert "voronoi.sql" (\i voronoi.sql). If you want additional 
information for each polygon, add variables to "points.txt" and adapt 
the "scan" statement accordingly. An additional variable (e.g. "id") can 
be accessed in the loop as points$id[[i]]

If you have PL/R installed, you can access everything from within 
PostgreSQL. In that case there is no need for intermediate ascii-files: 
  data can be read and written directly with "pg.spi.exec(SQL-statement)".

PostgreSQL data can also be read directly into R by means of the Rdbi 
and Rdbi.PgSQL packages 
(http://www.bioconductor.org/repository/release1.5/package/html/Rdbi.html 
and 
http://www.bioconductor.org/repository/release1.5/package/html/RdbiPgSQL.html)
They are only available for Unix and can be somewhat tricky to install. 
Alternatively, you can use the RODBC package.

Whatever method you choose, text-files, Pl/R or Rdbi/RODBC, the 
combination of PostgreSQL and R will give a huge usability boost to your 
PostGIS applications.

Jan

Jan Hartmann
Department of Geography
University of Amsterdam