[postgis-users] Circle
João Gonçalves
joaofgo at gmail.com
Sat Jun 14 21:20:53 PDT 2008
Hi! I'm testing different solution to this problem using PL/R and
pgirmess R library. This library provides a mehod to approximate a
polygonal geometry to a circular one using n (user defined) vertexes.
Drawing up some equations to get a comparison ratio between the circle
area and it's polygonal approximation area I've found that about 200
vertexes suffice to provide an area difference of just 0.00001% :: [100
- ((polygon area / circle area) * 100)].
The implementation uses a set of points and draws a polycircle for each
one. By convention uses gid as PKey identifier in both point dataset (as
circle centers) and output table; output and points tables geometry
field is always named the_geom.
DROP FUNCTION IF EXISTS circle2polygon(TEXT, TEXT, FLOAT8, INT4, INT4)
CASCADE;
CREATE OR REPLACE FUNCTION circle2polygon(TEXT, TEXT, FLOAT8, INT4,
INT4) RETURNS TEXT AS '
# Uses a set of point features to generate a polygonal approximation
# to a circular geometry centered in the inputed points
# Parameters
# arg1: Input table name (points)
# arg2: Output table name (as polygon)
# arg3: Circle radius
# arg4: SRID value
# arg5: Number of vertexes (for the polygonal approximation)
library(pgirmess)
gids<-pg.spi.exec(sprintf("SELECT gid FROM %1$s;",arg1))
for (i in 1:length(gids$gid)){
# Get circle center points
points<-pg.spi.exec(sprintf("SELECT st_x(the_geom) as x,
st_y(the_geom) as y FROM %1$s WHERE gid = %2$i", arg1, gids$gid[[i]]))
# Get polycircle vertexes coordinates
coords<-polycirc(arg3, pts = c(points$x, points$y), nbr = arg5)
# Open the geometry string
geom_txt = "GeomFromText(''POLYGON(("
# Concatenate polygon vertexes
for(j in 1:length(coords[,1])){
# Vertexes
geom_txt = sprintf("%s %.6f %.6f,", geom_txt, coords[j,1], coords[j,2])
}
# Close geometry string with the start vertex
geom_txt = sprintf("%s %.6f %.6f))'', %i)", geom_txt, coords[1,1],
coords[1,2], arg4)
# Insert into the database
pg.spi.exec(sprintf("INSERT INTO %1$s (gid, the_geom) VALUES (%2$i,
%3$s)", arg2, gids$gid[[i]], geom_txt))
}
return("OK")
' LANGUAGE 'plr';
DROP TABLE IF EXISTS test_polycircle;
CREATE TABLE "public"."test_polycircle" (
"gid" INTEGER NOT NULL,
"the_geom" GEOMETRY,
PRIMARY KEY("gid")
) WITH OIDS;
-- Test example:
SELECT circle2polygon('point_table_name', 'test_polycircle', 10.785, -1,
250);
With some adaptions this can be used for other applications besides this
one.
I'm hoping that someone can help me to improve this, probably making it
more robust.
More information about the postgis-users
mailing list