[postgis-users] Voronoi / Dalaunay function with R
Jan Hartmann
j.l.h.hartmann at uva.nl
Thu Apr 6 12:49:09 PDT 2006
Hi Mike,
You can change the boundary behavior of deldir in two ways:
- voro = deldir(points$x,points$y,rw=c(minx,maxx,miny,maxy))
changes the boundary rectangle of the voronoi diagram. By default,
this rectangle is 10% + the actual data range. If you set
minx, miny,maxx,maxx to large values, all data points will fall within
the diagram
- voro = deldir(x=points$x,y=points$y,dpl=list(ndx=2,ndy=2))
will add four dummy points at the edges of the bounding rectangle.
Note that this will result in four extra polygons. These are added at
the end of the polygon list.
- voro =
deldir(x=points$x,y=points$y,dpl=list(ndx=2,ndy=2),rw=c(minx,maxx,miny,maxy))
combines both methods. Experiment for good-looking results
The number of digits can be changed by the parameter "digits=" of
deldir. This does not seem to work however, probably a bug. You can
output numbers with an explicit C-style format like this:
for (j in 1:length(tile$x)) {
cat (sprintf("%.6f",tile$x[[j]]),' ',sprintf("%.6f",tile$y[[j]]),",")
}
Jan
Mike Leahy wrote:
> Jan,
>
> Thanks - that works pretty well. It's certainly way faster than my
> script. However, the polygons aren't quite right around the outlying
> area when I tried it with one of my datasets. In the image attached,
> you can see an example of the voronoi polygons from my script (the
> bigger darker polygons), and from the R-generated polygons using your
> script below (the smaller lighter polygons). You'll notice that in one
> case, there's a point that isn't even inside its own voronoi polygon
> because the output from R got clipped off at too short of a distance.
> I'm pretty new to R (i.e., this is the first time I've ever done
> anything besides install it), so I'm not sure how this could be
> resolved. Do you have any suggestions? I'm thinking maybe I could add
> four outlying corner-points located well-beyond the extent of the
> dataset to force the R calculations to fill-in the entire dataset, then
> clip the result to a given buffer distance.
>
> Also, there is a slight variation in the accuracy of the vertices - from
> what I can tell, R is rounding off the coordinates of the vertices to
> one decimal point, whereas the sql in the script I posted does no
> rounding/snapping at all. Would anyone know of a way to get around
> that? I guess that's really a question for the R-project list. Looks
> like I have a new learning curve to face here...
>
> Thanks again,
> Mike
>
> Jan Hartmann wrote:
>
>>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
>>
>>
>>
>>_______________________________________________
>>postgis-users mailing list
>>postgis-users at postgis.refractions.net
>>http://postgis.refractions.net/mailman/listinfo/postgis-users
>>
>>
>>
>> ------------------------------------------------------------------------
>>
>>
>>------------------------------------------------------------------------
>>
>>_______________________________________________
>>postgis-users mailing list
>>postgis-users at postgis.refractions.net
>>http://postgis.refractions.net/mailman/listinfo/postgis-users
More information about the postgis-users
mailing list