<div>Mark:</div>

<div>  Thank you. This functionality is obviously helpful to many people. I will try it and maybe I translate to a function to make it more generic (negative coordinates), eliminate the temporary table, do the loop and permit a destination table. I will post any improves I did.

</div>

<div>Thanks a lot again.<br><br>-- <br>Marco Vieira<br>+55 21 9499-6800<br>e-mail: <a href="mailto:maovieira@gmail.com">maovieira@gmail.com</a> <br> </div>

<div><span class="gmail_quote">2006/2/3, Mark Fenbers <<a href="mailto:Mark.Fenbers@noaa.gov">Mark.Fenbers@noaa.gov</a>>:</span>

<blockquote class="gmail_quote" style="PADDING-LEFT: 1ex; MARGIN: 0px 0px 0px 0.8ex; BORDER-LEFT: #ccc 1px solid">Attached (as promised yesterday) is an ASCII file giving the SQL to produce voronoi polygons given a set of points.

<br><br>Mark<br><br>marco vieira wrote:<br><br>> Hi list.<br>>   Someone writes a function to make a 2D Voronoi diagram from a set of<br>> points? Something like select voronoi(collect(geom)) from points<br>> I found many posts in 2004 about a Delaunay triangulation function and

> the last post from Erwan > <a href="http://postgis.refractions.net/pipermail/postgis-users/2004-October/005989.html">http://postgis.refractions.net/pipermail/postgis-users/2004-October/005989.html</a> > show me that something was implementing...

<br>><br>> --<br>> Marco Vieira<br>> +55 21 9499-6800<br>> e-mail: <a href="mailto:maovieira@gmail.com">maovieira@gmail.com</a> <mailto:<a href="mailto:maovieira@gmail.com">maovieira@gmail.com</a>><br>

><br>>------------------------------------------------------------------------<br>><br>>_______________________________________________<br>>postgis-users mailing list<br>><a href="mailto:postgis-users@postgis.refractions.net">

postgis-users@postgis.refractions.net</a><br>><a href="http://postgis.refractions.net/mailman/listinfo/postgis-users">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>><br>><br><br><br>-- written my Mark J Fenbers, OHRFC, June 2005

<br>-- needs an arg (:radar) from the command-line as in the following:<br>-- for radar in CLE BUF ILX CCX RLX LOT ILN DTX GRR IND LVX JKL OHX PAH FCX PBZ LSX MRX LWX IWX VWX ; do<br>--      psql -d hd_ob6tir -f ~/SQL/voronoi.sql -v radar="'$radar'";

<br>-- done<br><br>-- Description:  This SQL uses a table called radarloc which holds locations of radars<br>-- fields: radarID, lat, lon (positive westward), and a y/n flag whether the radar is in use<br>-- A new table is created (gisradarloc) which adds the lat/lon as a geometry and add a second

<br>-- geometry for the 125-nm ring around the point.  The original table is unaltered.<br>-- Results will go into gisRadarVoronoi.<br>-- Only one polygon ($radar) is added to the table per execution of this query.  This is why

-- this is typically called from a loop, such as the example above... DROP TABLE gisRadarLoc; SELECT * INTO gisradarloc FROM radarloc;  -- need radid, use_radar, and lat,lon fields from radaloc; other columns ignored...

<br>SELECT AddGeometryColumn('public','gisradarloc','the_geom',2792,'POINT',2); -- the center point geometry<br>SELECT AddGeometryColumn('public','gisradarloc','ring_geom',2792,'POLYGON',2); -- the "circle" geometry (radar umbrella)

<br><br>UPDATE gisradarloc SET the_geom = transform(SetSRID(makepoint(-1 * lon, lat), 4269), 2792);<br>UPDATE gisradarloc SET ring_geom = buffer(the_geom, 231500.0, 90); -- compute 125 nautical mile path around point<br>CREATE INDEX gisradarlocgix ON gisradarloc USING GIST (the_geom);

<br><br>DROP AGGREGATE MultiUnion (GEOMETRY);  -- if exists, drop the aggregate function<br>CREATE AGGREGATE MultiUnion (   -- builds the Union of geometries in all rows<br>       BASETYPE = GEOMETRY,<br>       SFUNC = GeomUnion,

<br>       STYPE = GEOMETRY<br>);<br><br>-- next, build a TEMP table to store the circle around the radar site,<br>--   and the LINESTRINGs that bisect the intersecting regions of circle of adjacent radars<br>SELECT radid, (SELECT ring_geom FROM gisRadarLoc WHERE radid = :radar) AS umb_geom,

<br>       LineFromMultiPoint(<br>               Intersection(<br>                       ExteriorRing(ring_geom),<br>                       ExteriorRing((SELECT ring_geom FROM gisRadarLoc WHERE radid = :radar))<br>               )

<br>       ) AS x_geom<br>       INTO TEMP MJF FROM gisradarloc WHERE radid != :radar AND use_radar = 'T'<br>       AND Intersects(ring_geom, (SELECT ring_geom FROM gisRadarLoc WHERE radid = :radar) );<br><br>-- Now use that LINESTRING an add the adjacent radar's center point

<br>UPDATE MJF SET x_geom = AddPoint(x_geom,(SELECT the_geom FROM gisRadarLoc WHERE gisRadarLoc.radid = MJF.radid),0);<br><br>-- Now compute (from the 3 existing points) the remaining corner points of the rectangle<br>UPDATE MJF SET x_geom = AddPoint(x_geom,MakePoint(2 * X(PointN(x_geom,1)) - X(PointN(x_geom,2)), 2 * Y(PointN(x_geom,1)) - Y(PointN(x_geom,2))));

<br>UPDATE MJF SET x_geom = AddPoint(x_geom,MakePoint(2 * X(PointN(x_geom,1)) - X(PointN(x_geom,3)), 2 * Y(PointN(x_geom,1)) - Y(PointN(x_geom,3))));<br><br>-- The convexHull() of these 5 points results in a rectangle.  2 adjacent corners intersect the circle of :radar

UPDATE MJF SET x_geom = ConvexHull((x_geom)); -- Now what is needed is to subtract the rectangles from the circle to make a voronoi polygon.  But since I can't guarantee the order that -- an aggregate function will do the subtraction, I use and aggregate to build a union() of all the rectangles (from all the adjacent radars).

<br>-- Then subtract the union of the rectangles from the :radar's circle in one operation.<br>SELECT MultiUnion(x_geom) AS Rectblob INTO TEMP WFOrects from MJF;<br><br>------- Do next 3 lines (uncomment them) if you are creating/inserting the first record into gisRadarVoronoi ------

<br>--DROP TABLE gisRadarVoronoi;  -- will drop table if it already exists<br>--CREATE TABLE gisRadarVoronoi (radid varchar(8)); -- create table if it does not exist...<br>--SELECT AddGeometryColumn('public','gisradarvoronoi','the_geom',4269,'POLYGON',2);

-- The insert will add one voronoi polygon geometry for :radar to the table... INSERT INTO gisRadarVoronoi values(        :radar,        (SELECT transform(                Difference(                        SnapToGrid( -- avoid disjoint lines due to precision errors at 13 decimal places...

(SELECT ring_geom FROM gisradarloc WHERE radid = :radar)                        ,0.000001),                        SnapToGrid( -- avoid disjoint lines due to precision errors at 13 decimal places...

<br>                               (SELECT Rectblob FROM WFOrects)<br>                       ,0.000001)<br>               )<br>               ,4269)<br>       )<br>);<br><br><br>_______________________________________________

<br>postgis-users mailing list<br><a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br><a href="http://postgis.refractions.net/mailman/listinfo/postgis-users">http://postgis.refractions.net/mailman/listinfo/postgis-users

</a><br><br><br></blockquote></div><br><br clear="all">