[postgis-users] Circle

[João Gonçalves]

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.