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

claudio cesar trevisani claudio.trevisani at pmmg.mg.gov.br
Sat Feb 4 11:36:44 PST 2006


Marco, 

I will wait for yours code.

Claudio Cesar Trevisani

Citando marco vieira <maovieira at gmail.com>:

> 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
> >
> >
> >
> 


-------------------------------------------------
This mail sent through IMP: http://horde.org/imp/




More information about the postgis-users mailing list