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

[Mike Leahy]

Hello again,

I figured out how to make the voronoi calculations work more
efficiently.  Rather than building rectangles based on the distance
between bisector lines and the neighbouring points, I rewrote the script
to expand the bisector lines into squares whose size is dependant on the
buffer threshold chosen, such that there are no gaps between the
rectangles associated with each neighbour point.  This worked great, as
it produced a single voronoi polygon for each point, and eliminated any
precision error reported by geos.  However, I then noticed tiny holes
between the polygons at their corners.  It dawned on me (somewhat
belatedly) that this was due to precision error based on the fact that
the buffer function only approximates curves, and therefore the
intersections between different buffers wouldn't provide precise
coordinates.

To fix this, I dropped the method of intersecting buffers.  Instead, I
just find the midpoint between each point and each of its
neighbours...calculate the endpoints of a line extending from that point
to bisect the space between the points (producing a more exact version
of the line obtained from intersected buffer rings), then I build
squares on the neighbour-side of these lines as before.  These
neighbour-squares are then unioned, and subtracted from a buffer around
the central point.  So far, I've had no precision errors from geos, and
I can't see any holes between the resulting polygons.  The processing is
also slightly faster.

I've attached a new version of the sql that does this to this email.
I'm fairly sure this script could be reworked to be much more
computationally efficient - I'm pretty sure my code is doing a great
deal of redundant calculations as it is now, but the results are fine
and are relatively fast for small datasets.  I'm sure that this would be
easily scriptable either as a standalone executable, or maybe even as a
function in PostGIS if someone with a better mind for that were to take
it on.

Regards,
Mike

Mike Leahy wrote:
> 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
>>
>> ------------------------------------------------------------------------
>>
>> --Writtent by Michael G. Leahy (mgleahy {at} alumni.uwaterloo.ca), WLU/UW, April 4, 2006
>>
>> --this is a semi-automated procedure based on Mark Fenbers' script posted to the PostGIS-users mailing list
>> --    back in February will create voronoi polygons for an entire set of points in PostGIS...
>>
>> --requirements:  a source table (i.e., TABLE_NAME), with a unique identifier 'gid', and a point geometry 'the_geom'.
>>
>> --the data from the source table are copied to a new table named voronoi_output in this example.  
>> --Buffers are created for all points.
>> --two points from the instersection of the buffer outlines for each point's neighbour and its own buffer are 
>> --    converted into lines that bisect the space between the points, and stored in a new table with the bisecting 
>> --    line,  the point's gid field, and the neighouring point's gid
>> --Using Mark Fenbers' approach, these bisectors are converted into rectangles that cover the area between
>> --    the bisecting lines and the neighbour points, such that when all of these rectangles for a given point are
>> --    unioned, the space left around the central point is the area closest to that point.
>> --The unioned geometry of these bisector-based rectangles are then subtracted from the buffer around the central
>> --    point, which ideally leaves the voronoi polygon around the central point in the middle.  
>>
>> --Unfortunately, if a dataset  is irregularly distributed, it is likely that what is left from this operation for each point will 
>> --    be a multipolygon that contains the voronoi polygon that you're looking for, plus a bunch of other polgyons that resulted 
>> --    from many gaps between rectangles covering the central point's buffer.  Hence, my not-so-amazing solution to get the 
>> --    desired voronoi out of this multipolygon using the intersection between the centre point and incremeting GeometryN() 
>> --    results from the mutlipolygon.  If it were somehow possible to identify the index of the first intersecting polygon out of a 
>> --    multipolygon (e.g., GeometryN(multigeom,IntersectsN(multigeom,geom)), then this could be accomplished easliy with one
>> --    statement, rather than iteratively searching through the multipolygon.  Maybe someone can figure out how to automate this
>> --    better.  Another way to improve this might be to make the bisector-based rectangles much larger, so that there are fewer holes,
>> --    and perhaps fewer precision errors from Geos...
>>
>> DROP TABLE voronoi_output;
>> drop table voronoi_bisectors;
>> drop table voronoi_inverse;
>> select * into voronoi_output from {TABLE_NAME};
>> alter table voronoi_output add column id serial primary key;  --this lets me load the data in Quantum GIS...
>> SELECT AddGeometryColumn('public','voronoi_output','ring_geom',32718,'POLYGON',2); -- this will store the buffer around the point...
>> SELECT AddGeometryColumn('public','voronoi_output','voronoi_geom',32718,'POLYGON',2); -- this will store ultimately store the desired voronoi polygon.
>> SELECT AddGeometryColumn('public','voronoi_output','voronoi_multi_geom',32718,'MULTIPOLYGON',2); -- this will store the first remanants of the difference of each buffer from the bisector rectangles for all neighbouring points.
>>
>> -- compute buffer for the points (3km in this example, with 90 points per quarter-circle approximation...).  Using this approach,  bigger 
>> -- buffers means exponentially more points will be considered neighbours, so an appropriate point size is beneficial...
>> UPDATE voronoi_output SET ring_geom = buffer(the_geom, 3000, 90); 
>> CREATE INDEX voronoi_output_geom_idx ON voronoi_output USING GIST (the_geom);
>> CREATE INDEX voronoi_output_ring_idx ON voronoi_output USING GIST (ring_geom);
>>
>> -- build a query that gets the bisector between each point and any other points where their buffers touch each other...
>> select a.gid as agid, b.gid as bgid,  
>> LineFromMultiPoint(Intersection(ExteriorRing(a.ring_geom),ExteriorRing(b.ring_geom))) AS bisector_geom
>> into voronoi_bisectors
>> from voronoi_output a cross join voronoi_output b 
>> where a.ring_geom && b.ring_geom and intersects(a.ring_geom,b.ring_geom) and a.gid != b.gid;
>> alter table voronoi_bisectors add column id serial primary key; --again, this is so I can see the data in QGIS
>>
>> -- append the neighbouring point to each bisector linestring:
>> UPDATE voronoi_bisectors SET bisector_geom = AddPoint(bisector_geom,(SELECT the_geom FROM voronoi_output WHERE voronoi_output.gid = bgid),0);
>>
>> -- Now compute (from the 3 existing points) the remaining corner points of the rectangle associated with each bisector/neighbour triangle
>> UPDATE voronoi_bisectors SET bisector_geom = AddPoint(bisector_geom,MakePoint(2 * X(PointN(bisector_geom,1)) - X(PointN(bisector_geom,2)), 2 * Y(PointN(bisector_geom,1)) - Y(PointN(bisector_geom,2))));
>> UPDATE voronoi_bisectors SET bisector_geom = AddPoint(bisector_geom,MakePoint(2 * X(PointN(bisector_geom,1)) - X(PointN(bisector_geom,3)), 2 * Y(PointN(bisector_geom,1)) - Y(PointN(bisector_geom,3))));
>>
>> -- The convexHull() of these 5 points results in a rectangle covering an area adjacent to the centre point's buffer, with one side being the original bisector with endpoints intersecting the buffer.
>> SELECT AddGeometryColumn('public','voronoi_bisectors','bisector_rect_geom',32718,'POLYGON',2); -- this will store the buffer around the point...
>> UPDATE voronoi_bisectors SET bisector_rect_geom = ConvexHull((bisector_geom));
>>
>> select gid, (select geomunion(bisector_rect_geom) from voronoi_bisectors where voronoi_output.gid=agid) as inv_voronoi into voronoi_inverse from voronoi_output;
>> alter table voronoi_inverse add column id serial primary key; -- aslo so that I can see the data in QGIS
>>
>> --Subtract the unioned neighbour rectangles from each point's buffer.  As Mark Fenbers did in his script, the buffers and rectangles were 
>> --    snapped to a fine grid to avoid precision errors from disjoint lines calculated at 13 decimal places:
>> update voronoi_output set voronoi_multi_geom = multi(difference(snaptogrid(ring_geom,0.000001),snaptogrid(inv_voronoi,0.000001))) from voronoi_inverse where voronoi_inverse.gid = voronoi_output.gid;
>>
>> --add a counter field to the table, and set it to 1:
>> alter table voronoi_output add column i integer;
>> update voronoi_output set i = 1;
>>
>> --  Here's the downer: repeate the next three lines until the count returned from the third query reaches zero - the more points you have, and 
>> --    the bigger your buffers are, the more times you'll have to do this, but it should work - I just keep pasting them into the psql console until all 
>> --    voronoi polygons are found, but I'm sure this can be scripted relatively easily:
>> update voronoi_output set voronoi_geom = geometryn(voronoi_multi_geom,i) where voronoi_geom is null and intersects(the_geom,geometryn(voronoi_multi_geom,i));
>> update voronoi_output set i = i + 1 where voronoi_geom is null;
>> select count(*) from voronoi_output where voronoi_geom is null;
>>
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> 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: voronoi3.sql
Type: text/x-sql
Size: 6050 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20060405/a7c42536/attachment.bin>