[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