[postgis-users] Buffers as kms

[Andreas Forø Tollefsen]

Hi all,

I need to buffer a set of point features by using a radius column which
indicate the buffer radius for each point in kilometers.

I tried to first create points in lat/long WGS 84, then use the st_transform
and buffer using ECKERT VI which should be equal area meter as far as i
know.
However, it does not turn out correctly.
Any idea what i do wrong in the script below?

I first inserted the SRID reference:
INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext)
values ( 953010, 'esri', 53010, '+proj=eck6 +lon_0=0 +x_0=0 +y_0=0
+a=6371000 +b=6371000 +units=m +no_defs ',
'PROJCS["Sphere_Eckert_VI",GEOGCS["GCS_Sphere",DATUM["Not_specified_based_on_Authalic_Sphere",SPHEROID["Sphere",6371000,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Eckert_VI"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],UNIT["Meter",1],AUTHORITY["EPSG","53010"]]');

My script for buffering (ignore the python syntax, it is done using
PyGresql):

import pg
db1 = pg.connect('postgis', 'localhost', 5432, None, None, 'postgres',
'G83rec')

# Copy csv file into database
dropafr = "DROP TABLE IF EXISTS afr;"
db1.query(dropafr)
createtable = "CREATE TABLE afr (id int, year int, lat decimal, long
decimal, radius decimal, confsite varchar, confterr varchar, version
varchar);"
db1.query(createtable)
copyfromcsv = "COPY afr FROM 'c:/prio_grid/source/confsite/africa.csv'
DELIMITER ';';"
db1.query(copyfromcsv)
addgeomfield = "SELECT addgeometrycolumn('public', 'afr', 'centroid', 4326,
'POINT', 2);"
db1.query(addgeomfield)
fillgeomfield = "UPDATE afr SET centroid=ST_SetSRID(ST_MakePoint(long, lat),
4326) WHERE long != -99 OR lat != -99;"
db1.query(fillgeomfield)
addprimkey = "ALTER TABLE afr ADD column gid serial primary key;"
db1.query(addprimkey)

# Create list for for loop
yearlist = range(1989, 2009, 1)
for x in yearlist:
    dropsplit = "DROP TABLE IF EXISTS afr"+str(x)+";"
    db1.query(dropsplit)
    splittable = "SELECT * INTO afr"+str(x)+" FROM afr WHERE year =
"+str(x)+";"
    db1.query(splittable)
    addgeom = "SELECT addgeometrycolumn('public', 'afr"+str(x)+"', 'buffer',
953010, 'POLYGON', 2);"
    db1.query(addgeom)
    updatebuffer = "UPDATE afr"+str(x)+" SET buffer =
ST_Buffer(ST_Transform(centroid, 953010),radius) WHERE radius > 0;"
    db1.query(updatebuffer)

Andreas
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20100927/a15d27c9/attachment.html>