[postgis-users] Voronoi / Dalaunay function (again)

[Mike Leahy]

Hello list,

I was hoping to get PostGIS to do voronoi polygons, and noticed this
posting from Mark Fenbers.  I thought I'd give it a shot, and took it
upon myself to rewrite the script to be a bit more generic, and so that
that it could do voronoi polygons for an entire dataset of points.  All
that is needed to use the sql in the attached file is to substitute your
source PostGIS table in 4th line of SQL (provided that it has a unique
'gid' field, and points stored in 'the_geom'), and maybe some
customization of values used in the later functions (e.g., srid, buffer
size, precision, etc...).

Mark Fenbers' method seems to work fairly well from tests with some of
my data, except for one problem for which I made a not-so-pretty, but
fairly easy workaround.  The issue I encountered was that with my point
datasets (which have inconsistent dispersion patterns), I needed to
calculate buffers at distances that were relatively large compared to
the density of the points in some areas.  Thus, I ended up producing
multipolygons that contained one voronoi polygon, and a bunch of
left-over polygons from gaps between rectangles calculated using Mark
Fenbers' approach, which intersected with the buffer around the central
point of the voronoi polygons.

To fix this, added a counter field 'i' to the table, and iteratively
checked to see if each point intersected with
geometryn(voronoi_multi_geom,i) in order to find the correct polygon
from the output.  If the the polygon at a given index intersected the
corresponding centre point, I saved it to the a single-polygon field.
Then I incremented the counter field in the table and checked the
polygon at the next index, and so on, until the voronoi polygon was
extracted for every point.  This method is described in a bit more
detail in the attached sql file.  Perhaps somebody can think of a better
way to automate this.  If there were some easy way to find out the index
of a shape in a multigeometry that intersects some other geometry, that
would resolve the issue.

At any rate, I hope that others will find this useful, and/or improve
upon it.  Thanks also to Mark Fenbers for sharing his script.

Mike

Mark Fenbers wrote:
> Attached (as promised yesterday) is an ASCII file giving the SQL to
> produce voronoi polygons given a set of points.
> 
> Mark
> 
> marco vieira wrote:
> 
>> Hi list.
>>   Someone writes a function to make a 2D Voronoi diagram from a set of
>> points? Something like select voronoi(collect(geom)) from points
>> I found many posts in 2004 about a Delaunay triangulation function and
>> the last post from Erwan
>> http://postgis.refractions.net/pipermail/postgis-users/2004-October/005989.html
>> show me that something was implementing...
>>
>> -- 
>> Marco Vieira
>> +55 21 9499-6800
>> e-mail: maovieira at gmail.com <mailto:maovieira at gmail.com>
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users at postgis.refractions.net
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>  
>>
> 
> ------------------------------------------------------------------------
> 
> -- written my Mark J Fenbers, OHRFC, June 2005
> -- needs an arg (:radar) from the command-line as in the following: 
> -- 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 
> -- 	psql -d hd_ob6tir -f ~/SQL/voronoi.sql -v radar="'$radar'"; 
> -- done
> 
> -- Description:  This SQL uses a table called radarloc which holds locations of radars 
> -- fields: radarID, lat, lon (positive westward), and a y/n flag whether the radar is in use
> -- A new table is created (gisradarloc) which adds the lat/lon as a geometry and add a second
> -- geometry for the 125-nm ring around the point.  The original table is unaltered.
> -- Results will go into gisRadarVoronoi.
> -- 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...
> SELECT AddGeometryColumn('public','gisradarloc','the_geom',2792,'POINT',2); -- the center point geometry
> SELECT AddGeometryColumn('public','gisradarloc','ring_geom',2792,'POLYGON',2); -- the "circle" geometry (radar umbrella)
> 
> UPDATE gisradarloc SET the_geom = transform(SetSRID(makepoint(-1 * lon, lat), 4269), 2792);
> UPDATE gisradarloc SET ring_geom = buffer(the_geom, 231500.0, 90); -- compute 125 nautical mile path around point
> CREATE INDEX gisradarlocgix ON gisradarloc USING GIST (the_geom);
> 
> DROP AGGREGATE MultiUnion (GEOMETRY);  -- if exists, drop the aggregate function
> CREATE AGGREGATE MultiUnion (	-- builds the Union of geometries in all rows
> 	BASETYPE = GEOMETRY,
> 	SFUNC = GeomUnion,
> 	STYPE = GEOMETRY
> );
> 
> -- next, build a TEMP table to store the circle around the radar site, 
> --   and the LINESTRINGs that bisect the intersecting regions of circle of adjacent radars
> SELECT radid, (SELECT ring_geom FROM gisRadarLoc WHERE radid = :radar) AS umb_geom,
> 	LineFromMultiPoint(
> 		Intersection(
> 			ExteriorRing(ring_geom),
> 			ExteriorRing((SELECT ring_geom FROM gisRadarLoc WHERE radid = :radar))
> 		)
> 	) AS x_geom
> 	INTO TEMP MJF FROM gisradarloc WHERE radid != :radar AND use_radar = 'T' 
> 	AND Intersects(ring_geom, (SELECT ring_geom FROM gisRadarLoc WHERE radid = :radar) );
> 
> -- Now use that LINESTRING an add the adjacent radar's center point	
> UPDATE MJF SET x_geom = AddPoint(x_geom,(SELECT the_geom FROM gisRadarLoc WHERE gisRadarLoc.radid = MJF.radid),0);
> 
> -- Now compute (from the 3 existing points) the remaining corner points of the rectangle
> 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))));
> 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))));
> 
> -- 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).
> -- Then subtract the union of the rectangles from the :radar's circle in one operation.
> SELECT MultiUnion(x_geom) AS Rectblob INTO TEMP WFOrects from MJF;
> 
> ------- Do next 3 lines (uncomment them) if you are creating/inserting the first record into gisRadarVoronoi ------
> --DROP TABLE gisRadarVoronoi;  -- will drop table if it already exists
> --CREATE TABLE gisRadarVoronoi (radid varchar(8)); -- create table if it does not exist...
> --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...
> 				(SELECT Rectblob FROM WFOrects)
> 			,0.000001)
> 		)
> 		,4269)
> 	)
> );
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
-------------- next part --------------
A non-text attachment was scrubbed...
Name: voronoi.sql
Type: text/x-sql
Size: 6912 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20060404/d06ac5d5/attachment.bin>