[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