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 ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html