[postgis-users] Voronoi / Dalaunay function (again)
marco vieira
maovieira at gmail.com
Sat Feb 4 08:15:27 PST 2006
Mark:
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.
Thanks a lot again.
--
Marco Vieira
+55 21 9499-6800
e-mail: maovieira at gmail.com
2006/2/3, Mark Fenbers <Mark.Fenbers at noaa.gov>:
>
> 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 --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20060204/85cdaeda/attachment.html>
More information about the postgis-users
mailing list