[R] rounding of voronoi vertices using deldir()

Mike Leahy mgleahy at alumni.uwaterloo.ca
Thu Apr 6 05:44:58 CEST 2006


Hello list,

I'm just getting started with using R - I have been trying over the past
day or so to work out a method for generating voronoi polygons for
PostGIS using SQL.  I was able to put together a procedure which works
relatively well, but is somewhat inefficient.  Someone on the PostGIS
list pointed me to the deldir() function in R, for which I can export a
text file with x/y coordinates from a PostGIS table, and write an SQL
script with insert statements containing PostGIS ewkt-compatible
geometries representing the voronoi polygons generated by deldir().

The script I'm using is below.  I've added a very small frac parameter
to ensure none of the points are dropped from my dataset (which was
actually happening for me with some fairly dispersed clusters of
points), and a rectangular window that is much larger than the dataset
itself to ensure that the voronoi polygons extend far enough to actually
cover all of the points in my dataset.

The problem I'm having is that the vertices at the corners of these
polygons seem to have their coordinates rounded to one decimal place.
Based on the documentation, the 'digits' option should allow me to
change this, but I'm not having getting anything different no matter
what I set it to.  Is there something obvious that I've missed?  I
realize that in most practical applications it's not a significant
issue, but I'd prefer more accuracy as I may be using voronoi polygons
to evaluate some statistical methods.

Below are some sample coordinates from my points.txt file, and the
script I'm running in R to write out the voronoi polygons.  If anyone
can see what's going, I'd be happy to hear it, as the R calculations are
much faster than the method I'm using within PostGIS/SQL:

points.txt:
 266662.462042329 | 8665540.49905558
 270171.836250559 |  8667802.6446983
 268895.572741816 | 8674257.75469324
 270054.378262961 | 8666483.37597101
 268402.641255299 | 8664853.87941629
 265707.056272354 | 8665434.09025432
 269985.118229025 | 8667743.14071004
 269282.034045422 | 8665403.39312076
....


R
  library(deldir)
  points = scan(file="points.txt",what=list(x=0.0,y=0.0),sep="|")
  voro =
deldir(points$x,points$y,digits=10,frac=0.000000000000001,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))))
  # 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 ("))',32718));\n")             # add SRID and newline
  }
  sink()                              # output to terminal
q()
n



Regards,
Mike




More information about the R-help mailing list