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