pggeodb.nancy.inra.fr/db_robot - db_robot on pggeodb.nancy.inra.fr

Function: addface

Schema

topology

Owner

albenard

Descriptions

args: toponame, apolygon, force_new=false - Registers a face primitive to a topology and gets its identifier.

Options

Definition

CREATE OR REPLACE FUNCTION topology.addface (
 atopology varchar,
 apoly public.geometry,
 force_new boolean = false
)
RETURNS integer AS
$span$
DECLARE
 bounds geometry;
 symdif geometry;
 faceid int;
 rec RECORD;
 rrec RECORD;
 relate text;
 right_edges int[];
 left_edges int[];
 all_edges geometry;
 old_faceid int;
 old_edgeid int;
 sql text;
 right_side bool;
 edgeseg geometry;
 p1 geometry;
 p2 geometry;
 p3 geometry;
 loc float8;
 segnum int;
 numsegs int;
BEGIN
 --
 -- Atopology and apoly are required
 --
 IF atopology IS NULL OR apoly IS NULL THEN
   RAISE EXCEPTION 'Invalid null argument';
 END IF;

 --
 -- Aline must be a polygon
 --
 IF substring(geometrytype(apoly), 1, 4) != 'POLY'
 THEN
   RAISE EXCEPTION 'Face geometry must be a polygon';
 END IF;

 for rrec IN SELECT (d).* FROM (
   SELECT ST_DumpRings(ST_ForceRHR(apoly)) d
 ) foo
 LOOP -- {
   --
   -- Find all bounds edges, forcing right-hand-rule
   -- to know what's left and what's right...
   --
   bounds = ST_Boundary(rrec.geom);

   sql := 'SELECT e.geom, e.edge_id, e.left_face, e.right_face FROM '
     || quote_ident(atopology)
     || '.edge e, (SELECT $1 as geom) r WHERE r.geom && e.geom'
   ;
   -- RAISE DEBUG 'SQL: %', sql;
   FOR rec IN EXECUTE sql USING bounds
   LOOP -- {
     --RAISE DEBUG 'Edge % has bounding box intersection', rec.edge_id;

     -- Find first non-empty segment of the edge
     numsegs = ST_NumPoints(rec.geom);
     segnum = 1;
     WHILE segnum < numsegs LOOP
       p1 = ST_PointN(rec.geom, segnum);
       p2 = ST_PointN(rec.geom, segnum+1);
       IF ST_Distance(p1, p2) > 0 THEN
         EXIT;
       END IF;
       segnum = segnum + 1;
     END LOOP;

     IF segnum = numsegs THEN
       RAISE WARNING 'Edge % is collapsed', rec.edge_id;
       CONTINUE; -- we don't want to spend time on it
     END IF;

     edgeseg = ST_MakeLine(p1, p2);

     -- Skip non-covered edges
     IF NOT ST_Equals(p2, ST_EndPoint(rec.geom)) THEN
       IF NOT ( _ST_Intersects(bounds, p1) AND _ST_Intersects(bounds, p2) )
       THEN
         --RAISE DEBUG 'Edge % has points % and % not intersecting with ring bounds', rec.edge_id, st_astext(p1), st_astext(p2);
         CONTINUE;
       END IF;
     ELSE
       -- must be a 2-points only edge, let's use Covers (more expensive)
       IF NOT _ST_Covers(bounds, edgeseg) THEN
         --RAISE DEBUG 'Edge % is not covered by ring', rec.edge_id;
         CONTINUE;
       END IF;
     END IF;

     p3 = ST_StartPoint(bounds);
     IF ST_DWithin(edgeseg, p3, 0) THEN
       -- Edge segment covers ring endpoint, See bug #874
       loc = ST_LineLocatePoint(edgeseg, p3);
       -- WARNING: this is as robust as length of edgeseg allows...
       IF loc > 0.9 THEN
         -- shift last point down
         p2 = ST_LineInterpolatePoint(edgeseg, loc - 0.1);
       ELSIF loc < 0.1 THEN
         -- shift first point up
         p1 = ST_LineInterpolatePoint(edgeseg, loc + 0.1);
       ELSE
         -- when ring start point is in between, we swap the points
         p3 = p1; p1 = p2; p2 = p3;
       END IF;
     END IF;

     right_side = ST_LineLocatePoint(bounds, p1) <
                  ST_LineLocatePoint(bounds, p2);


     IF right_side THEN
       right_edges := array_append(right_edges, rec.edge_id);
       old_faceid = rec.right_face;
     ELSE
       left_edges := array_append(left_edges, rec.edge_id);
       old_faceid = rec.left_face;
     END IF;

     IF faceid IS NULL OR faceid = 0 THEN
       faceid = old_faceid;
       old_edgeid = rec.edge_id;
     ELSIF faceid != old_faceid THEN
       RAISE EXCEPTION 'Edge % has face % registered on the side of this face, while edge % has face % on the same side', rec.edge_id, old_faceid, old_edgeid, faceid;
     END IF;

     -- Collect all edges for final full coverage check
     all_edges = ST_Collect(all_edges, rec.geom);

   END LOOP; -- }
 END LOOP; -- }

 IF all_edges IS NULL THEN
   RAISE EXCEPTION 'Found no edges on the polygon boundary';
 END IF;


 --
 -- Check that all edges found, taken togheter,
 -- fully match the ring boundary and nothing more
 --
 -- If the test fail either we need to add more edges
 -- from the polygon ring or we need to split
 -- some of the existing ones.
 --
 bounds = ST_Boundary(apoly);
 IF NOT ST_isEmpty(ST_SymDifference(bounds, all_edges)) THEN
   IF NOT ST_isEmpty(ST_Difference(bounds, all_edges)) THEN
     RAISE EXCEPTION 'Polygon boundary is not fully defined by existing edges at or near point %', ST_AsText(ST_PointOnSurface(ST_Difference(bounds, all_edges)));
   ELSE
     RAISE EXCEPTION 'Existing edges cover polygon boundary and more at or near point % (invalid topology?)', ST_AsText(ST_PointOnSurface(ST_Difference(all_edges, bounds)));
   END IF;
 END IF;

 IF faceid IS NOT NULL AND faceid != 0 THEN
   IF NOT force_new THEN
     RETURN faceid;
   ELSE
   END IF;
 END IF;

 --
 -- Get new face id from sequence
 --
 FOR rec IN EXECUTE 'SELECT nextval(' ||
   quote_literal(
     quote_ident(atopology) || '.face_face_id_seq'
   ) || ')'
 LOOP
   faceid = rec.nextval;
 END LOOP;

 --
 -- Insert new face
 --
 EXECUTE 'INSERT INTO '
   || quote_ident(atopology)
   || '.face(face_id, mbr) VALUES('
   -- face_id
   || faceid || ','
   -- minimum bounding rectangle
   || '$1)'
   USING ST_Envelope(apoly);

 --
 -- Update all edges having this face on the left
 --
 IF left_edges IS NOT NULL THEN
   EXECUTE 'UPDATE '
   || quote_ident(atopology)
   || '.edge_data SET left_face = '
   || quote_literal(faceid)
   || ' WHERE edge_id = ANY('
   || quote_literal(left_edges)
   || ') ';
 END IF;

 --
 -- Update all edges having this face on the right
 --
 IF right_edges IS NOT NULL THEN
   EXECUTE 'UPDATE '
   || quote_ident(atopology)
   || '.edge_data SET right_face = '
   || quote_literal(faceid)
   || ' WHERE edge_id = ANY('
   || quote_literal(right_edges)
   || ') ';
 END IF;


 --
 -- Set left_face/right_face of any contained edge
 --
 EXECUTE 'UPDATE '
   || quote_ident(atopology)
   || '.edge_data SET right_face = '
   || quote_literal(faceid)
   || ', left_face = '
   || quote_literal(faceid)
   || ' WHERE ST_Contains($1, geom)'
   USING apoly;

 --
 -- Set containing_face of any contained node
 --
 EXECUTE 'UPDATE '
   || quote_ident(atopology)
   || '.node SET containing_face = '
   || quote_literal(faceid)
   || ' WHERE containing_face IS NOT NULL AND ST_Contains($1, geom)'
   USING apoly;

 RETURN faceid;

END
$span$
LANGUAGE 'plpgsql'
VOLATILE
CALLED ON NULL INPUT
SECURITY INVOKER
COST 100;

COMMENT ON FUNCTION topology.addface(atopology varchar, apoly public.geometry, force_new boolean)
IS 'args: toponame, apolygon, force_new=false - Registers a face primitive to a topology and gets its identifier.';

This file was generated with SQL Manager for PostgreSQL (www.pgsqlmanager.com) at 26/02/2014 11:51