2010-12-15 13:39:52 +00:00
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.refractions.net
--
2011-01-25 18:08:57 +00:00
-- Copyright (C) 2010, 2011 Sandro Santilli <strk@keybit.net>
2010-12-15 13:39:52 +00:00
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- Functions used to populate a topology
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2010-12-17 16:16:41 +00:00
--
-- Developed by Sandro Santilli <strk@keybit.net>
-- for Faunalia (http://www.faunalia.it) with funding from
-- Regione Toscana - Sistema Informativo per la Gestione del Territorio
-- e dell' Ambiente [RT-SIGTA].
-- For the project: "Sviluppo strumenti software per il trattamento di dati
-- geografici basati su QuantumGIS e Postgis (CIG 0494241492)"
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2010-12-15 13:39:52 +00:00
- - {
--
-- AddNode(atopology, point)
--
-- Add a node primitive to a topology and get it's identifier.
-- Returns an existing node at the same location, if any.
--
2010-12-17 16:16:41 +00:00
-- When adding a _new_ node it checks for the existance of any
-- edge crossing the given point, raising an exception if found.
--
-- The newly added nodes have no containing face.
2010-12-15 13:39:52 +00:00
--
--
2010-12-17 16:16:41 +00:00
CREATE OR REPLACE FUNCTION topology . AddNode ( varchar , geometry )
2010-12-15 13:39:52 +00:00
RETURNS int
AS
2010-12-17 16:16:41 +00:00
$ $
2010-12-15 13:39:52 +00:00
DECLARE
atopology ALIAS FOR $ 1 ;
apoint ALIAS FOR $ 2 ;
nodeid int ;
rec RECORD ;
BEGIN
--
-- Atopology and apoint are required
--
IF atopology IS NULL OR apoint IS NULL THEN
2010-12-17 16:16:41 +00:00
RAISE EXCEPTION ' Invalid null argument ' ;
2010-12-15 13:39:52 +00:00
END IF ;
--
-- Apoint must be a point
--
2010-12-17 16:16:41 +00:00
IF substring ( geometrytype ( apoint ) , 1 , 5 ) ! = ' POINT '
2010-12-15 13:39:52 +00:00
THEN
2010-12-17 16:16:41 +00:00
RAISE EXCEPTION ' Node geometry must be a point ' ;
2010-12-15 13:39:52 +00:00
END IF ;
--
-- Check if a coincident node already exists
--
-- We use index AND x/y equality
--
2010-12-17 16:16:41 +00:00
FOR rec IN EXECUTE ' SELECT node_id FROM '
| | quote_ident ( atopology ) | | ' .node ' | |
' WHERE geom && ' | | quote_literal ( apoint : : text ) | | ' ::geometry '
2011-01-10 17:27:22 +00:00
| | ' AND ST_X(geom) = ST_X( ' | | quote_literal ( apoint : : text ) | | ' ::geometry) '
| | ' AND ST_Y(geom) = ST_Y( ' | | quote_literal ( apoint : : text ) | | ' ::geometry) '
2010-12-15 13:39:52 +00:00
LOOP
RETURN rec . node_id ;
END LOOP ;
--
2010-12-17 16:16:41 +00:00
-- Check if any edge crosses this node
-- (endpoints are fine)
2010-12-15 13:39:52 +00:00
--
2010-12-17 16:16:41 +00:00
FOR rec IN EXECUTE ' SELECT edge_id FROM '
| | quote_ident ( atopology ) | | ' .edge '
2011-01-17 22:35:50 +00:00
| | ' WHERE ST_DWithin( '
| | quote_literal ( apoint : : text )
| | ' , geom, 0) AND NOT ST_Equals( '
| | quote_literal ( apoint : : text )
| | ' , ST_StartPoint(geom)) AND NOT ST_Equals( '
| | quote_literal ( apoint : : text )
| | ' , ST_EndPoint(geom)) '
2010-12-15 13:39:52 +00:00
LOOP
2011-01-17 22:35:50 +00:00
RAISE EXCEPTION ' An edge crosses the given node. ' ;
2010-12-15 13:39:52 +00:00
END LOOP ;
--
-- Get new node id from sequence
--
2011-02-04 16:31:34 +00:00
FOR rec IN EXECUTE ' SELECT nextval( ' | |
quote_literal (
quote_ident ( atopology ) | | ' .node_node_id_seq '
) | | ' ) '
2010-12-15 13:39:52 +00:00
LOOP
nodeid = rec . nextval ;
END LOOP ;
--
-- Insert the new row
--
2010-12-17 16:16:41 +00:00
EXECUTE ' INSERT INTO ' | | quote_ident ( atopology )
| | ' .node(node_id, geom)
VALUES ( ' ||nodeid|| ' , ' ||quote_literal(apoint::text)||
' ) ' ;
2010-12-15 13:39:52 +00:00
RETURN nodeid ;
END
2010-12-17 16:16:41 +00:00
$ $
2011-02-07 16:53:26 +00:00
LANGUAGE ' plpgsql ' VOLATILE ;
2010-12-15 13:39:52 +00:00
- - } AddNode
2010-12-17 16:16:41 +00:00
2010-12-17 17:16:00 +00:00
- - {
--
-- AddEdge(atopology, line)
--
-- Add an edge primitive to a topology and get it's identifier.
2010-12-17 17:19:07 +00:00
-- Edge endpoints will be added as nodes if missing.
2011-01-25 09:39:50 +00:00
-- Returns an existing edge at the same location, if any.
2010-12-17 17:16:00 +00:00
--
-- An exception is raised if the given line crosses an existing
-- node or interects with an existing edge on anything but endnodes.
--
-- The newly added edge has "universe" face on both sides
-- and links to itself.
-- Calling code is expected to do further linking.
--
--
CREATE OR REPLACE FUNCTION topology . AddEdge ( varchar , geometry )
RETURNS int
AS
$ $
DECLARE
atopology ALIAS FOR $ 1 ;
aline ALIAS FOR $ 2 ;
edgeid int ;
rec RECORD ;
2011-01-27 08:48:11 +00:00
ix geometry ;
2010-12-17 17:16:00 +00:00
BEGIN
--
-- Atopology and apoint are required
--
IF atopology IS NULL OR aline IS NULL THEN
RAISE EXCEPTION ' Invalid null argument ' ;
END IF ;
--
-- Aline must be a linestring
--
IF substring ( geometrytype ( aline ) , 1 , 4 ) ! = ' LINE '
THEN
RAISE EXCEPTION ' Edge geometry must be a linestring ' ;
END IF ;
2011-01-18 17:19:55 +00:00
--
-- Check there's no face registered in the topology
--
FOR rec IN EXECUTE ' SELECT count(face_id) FROM '
| | quote_ident ( atopology ) | | ' .face '
| | ' WHERE face_id != 0 LIMIT 1 '
LOOP
IF rec . count > 0 THEN
RAISE EXCEPTION ' AddEdge can only be used against topologies with no faces defined ' ;
END IF ;
END LOOP ;
2010-12-17 17:16:00 +00:00
--
-- Check if the edge crosses an existing node
--
FOR rec IN EXECUTE ' SELECT node_id FROM '
| | quote_ident ( atopology ) | | ' .node '
| | ' WHERE ST_Crosses( '
| | quote_literal ( aline : : text ) | | ' ::geometry, geom '
| | ' ) '
LOOP
RAISE EXCEPTION ' Edge crosses node % ' , rec . node_id ;
END LOOP ;
--
-- Check if the edge intersects an existing edge
-- on anything but endpoints
--
-- Following DE-9 Intersection Matrix represent
-- the only relation we accept.
--
2010-12-20 22:01:04 +00:00
-- F F 1
-- F * *
-- 1 * 2
2010-12-17 17:16:00 +00:00
--
2010-12-20 22:01:04 +00:00
-- Example1: linestrings touching at one endpoint
2011-01-14 16:40:01 +00:00
-- FF1 F00 102
-- FF1 F** 1*2 <-- our match
2010-12-20 22:01:04 +00:00
--
-- Example2: linestrings touching at both endpoints
2011-01-14 16:40:01 +00:00
-- FF1 F0F 1F2
-- FF1 F** 1*2 <-- our match
2010-12-20 22:01:04 +00:00
--
2011-01-14 16:40:01 +00:00
FOR rec IN EXECUTE ' SELECT edge_id, geom, ST_Relate( '
| | quote_literal ( aline : : text )
| | ' ::geometry, geom) as im '
| | ' FROM '
2010-12-17 17:16:00 +00:00
| | quote_ident ( atopology ) | | ' .edge '
| | ' WHERE '
| | quote_literal ( aline : : text ) | | ' ::geometry && geom '
LOOP
2011-01-14 16:40:01 +00:00
IF ST_RelateMatch ( rec . im , ' FF1F**1*2 ' ) THEN
CONTINUE ; -- no interior intersection
END IF ;
--
-- Closed lines have no boundary, so endpoint
-- intersection would be considered interior
-- See http://trac.osgeo.org/postgis/ticket/770
--
-- Possible relate patterns:
-- FF1 0F0 1F2 : first line is open, second is closed
-- F01 FFF 102 : first line is closed, second is open
-- 0F1 FFF 1F2 : both first and second line are closed
--
-- Note that the boundary of closed line never intersects
-- (_F_ _F_ _F_ for first, ___ FFF ___ for second) so we
-- can use that pattern to tell that a line is closed
-- (only exceptional case would be in presence of an empty
-- line operand, which we should deal with before anyway)
--
-- The problem here is that we have interior/interior (last case)
-- or interior/boundary (first 2 cases) intersection,
-- we can tell if it's puntual (dimension 0) but can't tell if
-- it is _only_ on an endpoint w/out
-- _computing_
-- the actual intersection and comparing.
--
-- For sure we know that if we are facing such a case, such
-- intersection will have dimension 0 and there would be NO
-- intersections on the closed line boundary
--
IF ST_RelateMatch ( rec . im , ' FF10F01F2 ' ) THEN
-- first line (aline) is open, second (rec.geom) is closed
-- first boundary has puntual intersection with second interior
--
-- compute intersection, check it equals second endpoint
IF ST_Equals ( ST_Intersection ( rec . geom , aline ) ,
ST_StartPoint ( rec . geom ) )
THEN
RAISE DEBUG ' Edge shares boundary with existing closed edge % ' ,
rec . edge_id ;
CONTINUE ;
END IF ;
END IF ;
IF ST_RelateMatch ( rec . im , ' F01FFF102 ' ) THEN
-- second line (rec.geom) is open, first (aline) is closed
-- second boundary has puntual intersection with first interior
--
-- compute intersection, check it equals first endpoint
IF ST_Equals ( ST_Intersection ( rec . geom , aline ) ,
ST_StartPoint ( aline ) )
THEN
RAISE DEBUG ' Closed edge shares boundary with existing edge % ' ,
rec . edge_id ;
CONTINUE ;
END IF ;
END IF ;
IF ST_RelateMatch ( rec . im , ' 0F1FFF1F2 ' ) THEN
-- both lines are closed (boundary intersects nothing)
-- they have puntual intersection between interiors
--
-- compute intersection, check it's a single point
2011-01-17 08:47:34 +00:00
-- and equals first's StartPoint _and_ second's StartPoint
2011-01-14 16:40:01 +00:00
IF ST_Equals ( ST_Intersection ( rec . geom , aline ) ,
2011-01-17 08:47:34 +00:00
ST_StartPoint ( aline ) ) AND
ST_Equals ( ST_StartPoint ( aline ) , ST_StartPoint ( rec . geom ) )
2011-01-14 16:40:01 +00:00
THEN
RAISE DEBUG
' Closed edge shares boundary with existing closed edge % ' ,
rec . edge_id ;
CONTINUE ;
END IF ;
END IF ;
2011-01-25 09:39:50 +00:00
-- Reuse an EQUAL edge (be it closed or not)
IF ST_RelateMatch ( rec . im , ' 1FFF*FFF2 ' ) THEN
RAISE DEBUG ' Edge already known as % ' , rec . edge_id ;
RETURN rec . edge_id ;
END IF ;
2011-01-14 16:40:01 +00:00
2011-01-27 08:48:11 +00:00
-- WARNING: the constructive operation might throw an exception
BEGIN
ix = ST_Intersection ( rec . geom , aline ) ;
EXCEPTION
WHEN OTHERS THEN
RAISE NOTICE ' Could not compute intersection between input edge (%) and edge % (%) ' , aline : : text , rec . edge_id , rec . geom : : text ;
END ;
2011-01-14 16:40:01 +00:00
2011-01-27 08:48:11 +00:00
RAISE EXCEPTION ' Edge intersects (not on endpoints) with existing edge % at or near point % ' , rec . edge_id , ST_AsText ( ST_PointOnSurface ( ix ) ) ;
END LOOP ;
2010-12-17 17:16:00 +00:00
--
-- Get new edge id from sequence
--
2011-02-04 16:31:34 +00:00
FOR rec IN EXECUTE ' SELECT nextval( ' | |
quote_literal (
quote_ident ( atopology ) | | ' .edge_data_edge_id_seq '
) | | ' ) '
2010-12-17 17:16:00 +00:00
LOOP
edgeid = rec . nextval ;
END LOOP ;
--
-- Insert the new row
--
EXECUTE ' INSERT INTO '
| | quote_ident ( atopology )
| | ' .edge(edge_id, start_node, end_node, '
| | ' next_left_edge, next_right_edge, '
| | ' left_face, right_face, '
| | ' geom) '
| | ' VALUES( '
-- edge_id
| | edgeid | | ' , '
-- start_node
| | ' topology.addNode( '
| | quote_literal ( atopology )
| | ' , ST_StartPoint( '
| | quote_literal ( aline : : text )
| | ' )) , '
-- end_node
| | ' topology.addNode( '
| | quote_literal ( atopology )
| | ' , ST_EndPoint( '
| | quote_literal ( aline : : text )
| | ' )) , '
-- next_left_edge
| | edgeid | | ' , '
-- next_right_edge
| | edgeid | | ' , '
-- left_face
| | ' 0, '
-- right_face
| | ' 0, '
-- geom
| | quote_literal ( aline : : text )
| | ' ) ' ;
RETURN edgeid ;
END
$ $
2011-02-07 16:53:26 +00:00
LANGUAGE ' plpgsql ' VOLATILE ;
2010-12-17 17:16:00 +00:00
- - } AddEdge
2010-12-22 17:27:02 +00:00
- - {
--
-- AddFace(atopology, poly)
--
-- Add a face primitive to a topology and get it's identifier.
-- Returns an existing face at the same location, if any.
--
-- For a newly added face, its edges will be appropriately
-- linked (marked as left-face or right-face).
--
-- The target topology is assumed to be valid (containing no
-- self-intersecting edges).
--
-- An exception is raised if:
-- o The polygon boundary is not fully defined by existing edges.
-- o The polygon overlaps an existing face.
--
--
CREATE OR REPLACE FUNCTION topology . AddFace ( varchar , geometry )
RETURNS int
AS
$ $
DECLARE
atopology ALIAS FOR $ 1 ;
apoly ALIAS FOR $ 2 ;
bounds geometry ;
symdif geometry ;
faceid int ;
rec RECORD ;
relate text ;
right_edges int [ ] ;
left_edges int [ ] ;
all_edges geometry ;
old_faceid int ;
old_edgeid 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 ;
--
-- Find all bounds edges, forcing right-hand-rule
-- to know what's left and what's right...
--
bounds = ST_Boundary ( ST_ForceRHR ( apoly ) ) ;
FOR rec IN EXECUTE ' SELECT geom, edge_id, '
| | ' left_face, right_face, '
| | ' (st_dump(st_sharedpaths(geom, '
| | quote_literal ( bounds : : text )
| | ' ))).path[1] as side FROM '
| | quote_ident ( atopology ) | | ' .edge '
| | ' WHERE '
| | quote_literal ( bounds : : text )
| | ' ::geometry && geom AND ST_Relate(geom, '
| | quote_literal ( bounds : : text )
| | ' , '' 1FF****** '' ) '
2011-01-15 23:14:30 +00:00
| | ' GROUP BY geom, edge_id, left_face, right_face, side '
2010-12-22 17:27:02 +00:00
LOOP
RAISE DEBUG ' Edge % (left:%, right:%) - side : % ' ,
rec . edge_id , rec . left_face , rec . right_face , rec . side ;
IF rec . side = 1 THEN
-- This face is on the right
right_edges : = array_append ( right_edges , rec . edge_id ) ;
old_faceid = rec . right_face ;
ELSE
-- This face is on the left
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 ;
IF all_edges IS NULL THEN
RAISE EXCEPTION ' Found no edges on the polygon boundary ' ;
END IF ;
RAISE DEBUG ' Left edges: % ' , left_edges ;
RAISE DEBUG ' Right edges: % ' , right_edges ;
--
-- Check that all edges found, taken togheter,
-- fully match the polygon boundary and nothing more
--
-- If the test fail either we need to add more edges
-- from the polygon boundary or we need to split
-- some of the existing ones.
--
IF NOT ST_isEmpty ( ST_SymDifference ( bounds , all_edges ) ) THEN
IF NOT ST_isEmpty ( ST_Difference ( bounds , all_edges ) ) THEN
2011-02-09 16:36:45 +00:00
RAISE EXCEPTION ' Polygon boundary is not fully defined by existing edges at or near point % ' , ST_AsText ( ST_PointOnSurface ( ST_Difference ( bounds , all_edges ) ) ) ;
2010-12-22 17:27:02 +00:00
ELSE
2011-02-09 16:36:45 +00:00
RAISE EXCEPTION ' Existing edges cover polygon boundary and more at or near point % (invalid topology?) ' , ST_AsText ( ST_PointOnSurface ( ST_Difference ( all_edges , bounds ) ) ) ;
2010-12-22 17:27:02 +00:00
END IF ;
END IF ;
-- EXECUTE 'SELECT ST_Collect(geom) FROM'
-- || quote_ident(atopology)
-- || '.edge_data '
-- || ' WHERE edge_id = ANY('
-- || quote_literal(array_append(left_edges, right_edges))
-- || ') ';
--
-- TODO:
-- Check that NO edge is contained in the face ?
--
RAISE WARNING ' Not checking if face contains any edge ' ;
IF faceid IS NOT NULL AND faceid ! = 0 THEN
RAISE DEBUG ' Face already known as % ' , faceid ;
RETURN faceid ;
END IF ;
--
-- Get new face id from sequence
--
2011-02-04 16:31:34 +00:00
FOR rec IN EXECUTE ' SELECT nextval( ' | |
quote_literal (
quote_ident ( atopology ) | | ' .face_face_id_seq '
) | | ' ) '
2010-12-22 17:27:02 +00:00
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
2011-01-28 08:55:22 +00:00
| | quote_literal ( ST_Envelope ( apoly ) : : text )
2010-12-22 17:27:02 +00:00
| | ' ) ' ;
--
-- Update all edges having this face on the left
--
2010-12-22 18:45:07 +00:00
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 ;
2010-12-22 17:27:02 +00:00
--
-- Update all edges having this face on the right
--
2010-12-22 18:45:07 +00:00
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 ;
2010-12-22 17:27:02 +00:00
2011-01-17 08:47:24 +00:00
--
-- TODO:
-- Set next_left_face and next_right_face !
-- These are required by the model, but not really used
-- by this implementation...
--
RAISE WARNING ' Not updating next_{left,right}_face fields of face boundary edges ' ;
2010-12-22 17:27:02 +00:00
RETURN faceid ;
END
$ $
2011-02-07 16:53:26 +00:00
LANGUAGE ' plpgsql ' VOLATILE ;
2010-12-22 17:27:02 +00:00
- - } AddFace