-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2010, 2011 Sandro Santilli <strk@kbt.io>
-- Copyright (C) 2005 Refractions Research Inc.
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- Author: Sandro Santilli <strk@kbt.io>
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- STATUS:
--
--  All objects are created in the 'topology' schema.
--
--  We have PostGIS-specific objects and SQL/MM objects.
--  PostGIS-specific objects have no prefix, SQL/MM ones
--  have the ``ST_' prefix.
--
-- [PostGIS-specific]
--
-- TABLE topology
--  Table storing topology info (name, srid, precision)
--
-- TYPE TopoGeometry
--  Complex type storing topology_id, layer_id, geometry type
--  and topogeometry id.
--
-- DOMAIN TopoElement
--  An array of two elements: element_id and element_type.
--  In fact, an array of integers.
--
-- DOMAIN TopoElementArray
--  An array of element_id,element_type values.
--  In fact, a bidimensional array of integers:
--  '{{id,type}, {id,type}, ...}'
--
-- FUNCTION CreateTopology(name, [srid], [precision])
--  Initialize a new topology (creating schema with
--  edge,face,node,relation) and add a record into
--  the topology.topology table.
--  TODO: add triggers (or rules, or whatever) enforcing
--        precision to the edge and node tables.
--
-- FUNCTION DropTopology(name)
--    Delete a topology removing reference from the
--    topology.topology table
--
-- FUNCTION GetTopologyId(name)
-- FUNCTION GetTopologySRID(name)
-- FUNCTION GetTopologyName(id)
--    Return info about a Topology
--
-- FUNCTION AddTopoGeometryColumn(toponame, schema, table, column, geomtype)
--    Add a TopoGeometry column to a table, making it a topology layer.
--    Returns created layer id.
--
-- FUNCTION DropTopoGeometryColumn(schema, table, column)
--    Drop a TopoGeometry column, unregister the associated layer,
--    cleanup the relation table.
--
-- FUNCTION CreateTopoGeom(toponame, geomtype, layer_id, topo_objects)
--    Create a TopoGeometry object from existing Topology elements.
--    The "topo_objects" parameter is of TopoElementArray type.
--
-- FUNCTION GetTopoGeomElementArray(toponame, layer_id, topogeom_id)
-- FUNCTION GetTopoGeomElementArray(TopoGeometry)
--    Returns a TopoElementArray object containing the topological
--    elements of the given TopoGeometry.
--
-- FUNCTION GetTopoGeomElements(toponame, layer_id, topogeom_id)
-- FUNCTION GetTopoGeomElements(TopoGeometry)
--    Returns a set of TopoElement objects containing the
--    topological elements of the given TopoGeometry (primitive
--    elements)
--
-- FUNCTION ValidateTopology(toponame)
--    Run validity checks on the topology, returning, for each
--    detected error, a 3-columns row containing error string
--    and references to involved topo elements: error, id1, id2
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- Overloaded functions for TopoGeometry inputs
--
-- FUNCTION intersects(TopoGeometry, TopoGeometry)
-- FUNCTION equals(TopoGeometry, TopoGeometry)
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- FUNCTION TopoGeo_AddPoint(toponame, point)
--    Add a Point geometry to the topology
--    TODO: accept a topology/layer id
--          rework to use existing node if existent
--
-- FUNCTION TopoGeo_AddLinestring(toponame, line)
--    Add a LineString geometry to the topology
--    TODO: accept a topology/layer id
--          rework to use existing nodes/edges
--          splitting them if required
--
-- FUNCTION TopoGeo_AddPolygon(toponame, polygon)
--    Add a Polygon geometry to the topology
--    TODO: implement
--
-- TYPE GetFaceEdges_ReturnType
--    Complex type used to return tuples from ST_GetFaceEdges
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- [SQL/MM]
--
--   ST_InitTopoGeo
--    Done, can be modified to include explicit sequences or
--    more constraints. Very noisy due to implicit index creations
--    for primary keys and sequences for serial fields...
--
--   ST_CreateTopoGeo
--    Complete
--
--   ST_AddIsoNode
--    Complete
--
--   ST_RemoveIsoNode
--    Complete
--
--   ST_MoveIsoNode
--    Complete
--
--   ST_AddIsoEdge
--    Complete
--
--   ST_RemoveIsoEdge
--    Complete, exceptions untested
--
--   ST_ChangeEdgeGeom
--    Complete
--
--   ST_NewEdgesSplit
--    Complete
--    Also updates the Relation table
--
--   ST_ModEdgeSplit
--    Complete
--    Also updates the Relation table
--
--   ST_AddEdgeNewFaces
--    Complete
--    Also updates the Relation table
--
--   ST_AddEdgeModFace
--    Complete
--    Also updates the Relation table
--
--   ST_GetFaceEdges
--    Complete
--
--   ST_ModEdgeHeal
--    Complete
--    Also updates the Relation table
--
--   ST_NewEdgeHeal
--    Complete
--    Also updates the Relation table
--
--   ST_GetFaceGeometry
--    Implemented using ST_BuildArea()
--
--   ST_RemEdgeNewFace
--    Complete
--    Also updates the Relation table
--
--   ST_RemEdgeModFace
--    Complete
--    Also updates the Relation table
--
--   ST_ValidateTopoGeo
--    Unimplemented (probably a wrapper around ValidateTopology)
--
--

-- Uninstalling previous installation isn't really a good habit ...
-- Let people decide about that
-- DROP SCHEMA topology CASCADE;
















CREATE SCHEMA topology;
COMMENT ON SCHEMA topology IS 'PostGIS Topology schema';

-- Doing everything outside of a transaction helps
-- upgrading in the best case.
BEGIN;

--={ ----------------------------------------------------------------
--  POSTGIS-SPECIFIC block
--
--  This part contains function NOT in the SQL/MM specification
--
---------------------------------------------------------------------

--
-- Topology table.
-- Stores id,name,precision and SRID of topologies.
--
CREATE TABLE topology.topology (
  id SERIAL NOT NULL PRIMARY KEY,
  name VARCHAR NOT NULL UNIQUE,
  SRID INTEGER NOT NULL,
  precision FLOAT8 NOT NULL,
  hasz BOOLEAN NOT NULL DEFAULT false
);

--{ LayerTrigger()
--
-- Layer integrity trigger
--
CREATE OR REPLACE FUNCTION topology.LayerTrigger()
  RETURNS trigger
AS
$$
DECLARE
  rec RECORD;
  ok BOOL;
  toponame varchar;
  query TEXT;
BEGIN

  --RAISE NOTICE 'LayerTrigger called % % at % level', TG_WHEN, TG_OP, TG_LEVEL;

  IF TG_OP = 'INSERT' THEN
    RAISE EXCEPTION 'LayerTrigger not meant to be called on INSERT';
  ELSIF TG_OP = 'UPDATE' THEN
    RAISE EXCEPTION 'The topology.layer table cannot be updated';
  END IF;

  -- Check for existance of any feature column referencing
  -- this layer
  FOR rec IN SELECT * FROM pg_namespace n, pg_class c, pg_attribute a
    WHERE text(n.nspname) = OLD.schema_name
    AND c.relnamespace = n.oid
    AND text(c.relname) = OLD.table_name
    AND a.attrelid = c.oid
    AND text(a.attname) = OLD.feature_column
  LOOP
    query = 'SELECT * '
         ' FROM ' || quote_ident(OLD.schema_name)
      || '.' || quote_ident(OLD.table_name)
      || ' WHERE layer_id('
      || quote_ident(OLD.feature_column)||') '
         '=' || OLD.layer_id
      || ' LIMIT 1';
    --RAISE NOTICE '%', query;
    FOR rec IN EXECUTE query
    LOOP
      RAISE NOTICE 'A feature referencing layer % of topology % still exists in %.%.%', OLD.layer_id, OLD.topology_id, OLD.schema_name, OLD.table_name, OLD.feature_column;
      RETURN NULL;
    END LOOP;
  END LOOP;

  -- Get topology name
  SELECT name FROM topology.topology INTO toponame
    WHERE id = OLD.topology_id;

  IF toponame IS NULL THEN
    RAISE NOTICE 'Could not find name of topology with id %',
      OLD.layer_id;
  END IF;

  -- Check if any record in the relation table references this layer
  FOR rec IN SELECT c.oid FROM pg_namespace n, pg_class c
    WHERE text(n.nspname) = toponame AND c.relnamespace = n.oid
          AND c.relname = 'relation'
  LOOP
    query = 'SELECT * '
         ' FROM ' || quote_ident(toponame)
      || '.relation '
         ' WHERE layer_id = '|| OLD.layer_id
      || ' LIMIT 1';
    --RAISE NOTICE '%', query;
    FOR rec IN EXECUTE query
    LOOP
      RAISE NOTICE 'A record in %.relation still references layer %', toponame, OLD.layer_id;
      RETURN NULL;
    END LOOP;
  END LOOP;

  RETURN OLD;
END;
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
--} LayerTrigger()

--{
-- Layer table.
-- Stores topology layer informations
--
CREATE TABLE topology.layer (
  topology_id INTEGER NOT NULL
    REFERENCES topology.topology(id),
  layer_id integer NOT NULL,
  schema_name VARCHAR NOT NULL,
  table_name VARCHAR NOT NULL,
  feature_column VARCHAR NOT NULL,
  feature_type integer NOT NULL,
  level INTEGER NOT NULL DEFAULT 0,
  child_id INTEGER DEFAULT NULL,
  UNIQUE(schema_name, table_name, feature_column),
  PRIMARY KEY(topology_id, layer_id)
);

CREATE TRIGGER layer_integrity_checks BEFORE UPDATE OR DELETE
ON topology.layer FOR EACH ROW EXECUTE PROCEDURE topology.LayerTrigger();

--} Layer table.

--
-- TopoGeometry type
--
CREATE TYPE topology.TopoGeometry AS (
  topology_id integer,
  layer_id integer,
  id integer,
  type integer -- 1: [multi]point, 2: [multi]line,
               -- 3: [multi]polygon, 4: collection
);

--
-- TopoElement domain
--
-- This is an array of two elements: element_id and element_type.
--
-- When used to define _simple_ TopoGeometries,
-- element_type can be:
--  0: a node
--  1: an edge
--  2: a face
-- and element_id will be the node, edge or face identifier
--
-- When used to define _hierarchical_ TopoGeometries,
-- element_type will be the child layer identifier and
-- element_id will be composing TopoGeometry identifier
--
CREATE DOMAIN topology.TopoElement AS integer[]
  CONSTRAINT DIMENSIONS CHECK (
    array_upper(VALUE, 2) IS NULL
    AND array_upper(VALUE, 1) = 2
  );
ALTER DOMAIN topology.TopoElement ADD
        CONSTRAINT lower_dimension CHECK (
    array_lower(VALUE, 1) = 1
  );
ALTER DOMAIN topology.TopoElement DROP CONSTRAINT
  IF EXISTS
  type_range;
ALTER DOMAIN topology.TopoElement ADD
        CONSTRAINT type_range CHECK (
    VALUE[2] > 0
  );

--
-- TopoElementArray domain
--
-- Changed: 3.1.0 - get rid of IS NOT NULL check
CREATE DOMAIN topology.TopoElementArray AS integer[][]
  CONSTRAINT DIMENSIONS CHECK (
    array_upper(VALUE, 2) = 2
    AND array_upper(VALUE, 3) IS NULL
  );

-- Changed: 3.1.0
ALTER DOMAIN topology.TopoElementArray DROP CONSTRAINT
  IF EXISTS
  DIMENSIONS;

-- Changed: 3.1.0 - get rid of IS NOT NULL check
ALTER DOMAIN topology.TopoElementArray ADD
        CONSTRAINT type_range CHECK (
    array_upper(VALUE, 2) = 2
    AND array_upper(VALUE, 3) IS NULL
  );

--{ RelationTrigger()
--
-- Relation integrity trigger
--
CREATE OR REPLACE FUNCTION topology.RelationTrigger()
  RETURNS trigger
AS
$$
DECLARE
  toponame varchar;
  topoid integer;
  plyr RECORD; -- parent layer
  rec RECORD;
  ok BOOL;

BEGIN
  IF TG_NARGS != 2 THEN
    RAISE EXCEPTION 'RelationTrigger called with wrong number of arguments';
  END IF;

  topoid = TG_ARGV[0];
  toponame = TG_ARGV[1];

  --RAISE NOTICE 'RelationTrigger called % % on %.relation for a %', TG_WHEN, TG_OP, toponame, TG_LEVEL;

  IF TG_OP = 'DELETE' THEN
    RAISE EXCEPTION 'RelationTrigger not meant to be called on DELETE';
  END IF;

  -- Get layer info (and verify it exists)
  ok = false;
  FOR plyr IN EXECUTE 'SELECT * FROM topology.layer '
       'WHERE '
       ' topology_id = ' || topoid
    || ' AND'
       ' layer_id = ' || NEW.layer_id
  LOOP
    ok = true;
    EXIT;
  END LOOP;
  IF NOT ok THEN
    RAISE EXCEPTION 'Layer % does not exist in topology %',
      NEW.layer_id, topoid;
    RETURN NULL;
  END IF;

  IF plyr.level > 0 THEN -- this is hierarchical layer

    -- ElementType must be the layer child id
    IF NEW.element_type != plyr.child_id THEN
      RAISE EXCEPTION 'Type of elements in layer % must be set to its child layer id %', plyr.layer_id, plyr.child_id;
      RETURN NULL;
    END IF;

    -- ElementId must be an existent TopoGeometry in child layer
    ok = false;
    FOR rec IN EXECUTE 'SELECT topogeo_id FROM '
      || quote_ident(toponame) || '.relation '
         ' WHERE layer_id = ' || plyr.child_id
      || ' AND topogeo_id = ' || NEW.element_id
    LOOP
      ok = true;
      EXIT;
    END LOOP;
    IF NOT ok THEN
      RAISE EXCEPTION 'TopoGeometry % does not exist in the child layer %', NEW.element_id, plyr.child_id;
      RETURN NULL;
    END IF;

  ELSE -- this is a basic layer

    -- ElementType must be compatible with layer type
    IF plyr.feature_type != 4
      AND plyr.feature_type != NEW.element_type
    THEN
      RAISE EXCEPTION 'Element of type % is not compatible with layer of type %', NEW.element_type, plyr.feature_type;
      RETURN NULL;
    END IF;

    --
    -- Now lets see if the element is consistent, which
    -- is it exists in the topology tables.
    --

    --
    -- Element is a Node
    --
    IF NEW.element_type = 1
    THEN
      ok = false;
      FOR rec IN EXECUTE 'SELECT node_id FROM '
        || quote_ident(toponame) || '.node '
           ' WHERE node_id = ' || NEW.element_id
      LOOP
        ok = true;
        EXIT;
      END LOOP;
      IF NOT ok THEN
        RAISE EXCEPTION 'Node % does not exist in topology %', NEW.element_id, toponame;
        RETURN NULL;
      END IF;

    --
    -- Element is an Edge
    --
    ELSIF NEW.element_type = 2
    THEN
      ok = false;
      FOR rec IN EXECUTE 'SELECT edge_id FROM '
        || quote_ident(toponame) || '.edge_data '
           ' WHERE edge_id = ' || abs(NEW.element_id)
      LOOP
        ok = true;
        EXIT;
      END LOOP;
      IF NOT ok THEN
        RAISE EXCEPTION 'Edge % does not exist in topology %', NEW.element_id, toponame;
        RETURN NULL;
      END IF;

    --
    -- Element is a Face
    --
    ELSIF NEW.element_type = 3
    THEN
      IF NEW.element_id = 0 THEN
        RAISE EXCEPTION 'Face % cannot be associated with any feature', NEW.element_id;
        RETURN NULL;
      END IF;
      ok = false;
      FOR rec IN EXECUTE 'SELECT face_id FROM '
        || quote_ident(toponame) || '.face '
           ' WHERE face_id = ' || NEW.element_id
      LOOP
        ok = true;
        EXIT;
      END LOOP;
      IF NOT ok THEN
        RAISE EXCEPTION 'Face % does not exist in topology %', NEW.element_id, toponame;
        RETURN NULL;
      END IF;
    END IF;

  END IF;

  RETURN NEW;
END;
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
--} RelationTrigger()

--{
--  DropTopoGeometryColumn(schema, table, colum)
--
--  Drop a TopoGeometry column, unregister the associated layer,
--  cleanup the relation table.
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.DropTopoGeometryColumn(schema varchar, tbl varchar, col varchar)
  RETURNS text
AS
$BODY$
DECLARE
  rec RECORD;
  lyrinfo RECORD;
  ok BOOL;
  result text;
  sql TEXT;
BEGIN

  -- Get layer and topology info

  sql := $$
    SELECT t.name as toponame, l.*
    FROM topology.topology t, topology.layer l
    WHERE l.topology_id = t.id
    AND l.schema_name = $1
    AND l.table_name = $2
    AND l.feature_column = $3
  $$;

  ok := false;
  FOR rec IN EXECUTE sql USING schema, tbl, col
  LOOP
    ok := true;
    lyrinfo := rec;
  END LOOP;

  -- Layer not found
  IF NOT ok THEN
    RAISE EXCEPTION 'No layer registered on %.%.%',
      schema,tbl,col;
  END IF;

  -- Cleanup the relation table (if it exists)
  BEGIN
    sql := format(
      'DELETE FROM %I.relation WHERE layer_id = $1',
      lyrinfo.toponame
    );
    EXECUTE sql USING lyrinfo.layer_id;
  EXCEPTION
    WHEN UNDEFINED_TABLE THEN
      RAISE NOTICE '%', SQLERRM;
    WHEN OTHERS THEN
      RAISE EXCEPTION 'Got % (%)', SQLERRM, SQLSTATE;
  END;

  -- Drop the sequence for topogeoms in this layer
  sql := format(
    'DROP SEQUENCE IF EXISTS %I.topogeo_s_%s',
    lyrinfo.toponame,
    lyrinfo.layer_id
  );
  EXECUTE sql;

  ok = false;
  FOR rec IN SELECT * FROM pg_namespace n, pg_class c, pg_attribute a
    WHERE text(n.nspname) = schema
    AND c.relnamespace = n.oid
    AND text(c.relname) = tbl
    AND a.attrelid = c.oid
    AND text(a.attname) = col
  LOOP
    ok = true;
    EXIT;
  END LOOP;

  IF ok THEN
    -- Drop the layer column
    sql := format(
      'ALTER TABLE %I.%I DROP %I CASCADE',
      schema, tbl, col
    );
    EXECUTE sql;
  END IF;

  -- Delete the layer record
  sql := $$
    DELETE FROM topology.layer
    WHERE topology_id = $1
    AND layer_id = $2
  $$;
  EXECUTE sql USING lyrinfo.topology_id, lyrinfo.layer_id;


  result := format(
    'Layer %s (%I.%I.%I) dropped',
    lyrinfo.layer_id, schema, tbl, col
  );

  RETURN result;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE;
--
--} DropTopoGeometryColumn

-- {
--
-- populate_topology_layer
--
-- Register missing layers into topology.topology, looking at
-- their constraints.
--
-- The function doesn't attempt to determine if a layer is
-- hierarchical or primitive, but always assumes primitive.
--
-- }{
DROP FUNCTION IF EXISTS topology.populate_topology_layer();
CREATE OR REPLACE FUNCTION topology.populate_topology_layer()
	RETURNS TABLE(schema_name text, table_name text, feature_column text)
AS
$$
  INSERT INTO topology.layer
  WITH checks AS (
  SELECT
    n.nspname sch, r.relname tab,
    replace(c.conname, 'check_topogeom_', '') col,
    --c.consrc src,
    regexp_matches(c.consrc,
      '\.topology_id = (\d+).*\.layer_id = (\d+).*\.type = (\d+)') inf
  FROM (SELECT conname, connamespace, conrelid, conkey, pg_get_constraintdef(oid) As consrc
		    FROM pg_constraint) AS c, pg_class r, pg_namespace n
  WHERE c.conname LIKE 'check_topogeom_%'
    AND r.oid = c.conrelid
    AND n.oid = r.relnamespace
  ), newrows AS (
    SELECT inf[1]::int as topology_id,
           inf[2]::int as layer_id,
          sch, tab, col, inf[3]::int as feature_type --, src
    FROM checks c
    WHERE NOT EXISTS (
      SELECT * FROM topology.layer l
      WHERE l.schema_name = c.sch
        AND l.table_name = c.tab
        AND l.feature_column = c.col
    )
  )
  SELECT topology_id, layer_id, sch,
         tab, col, feature_type,
         0, NULL
  FROM newrows RETURNING schema_name,table_name,feature_column;
$$
LANGUAGE 'sql' VOLATILE;

--{
-- CreateTopoGeom(topology_name, topogeom_type, layer_id, elements)
--
-- Create a TopoGeometry object from Topology elements.
-- The elements parameter is a two-dimensional array.
-- Every element of the array is either a Topology element represented by
-- (id, type) or a TopoGeometry element represented by (id, layer).
-- The actual semantic depends on the TopoGeometry layer, either at
-- level 0 (elements are topological primitives) or higer (elements
-- are TopoGeoms from child layer).
--
-- @param toponame Topology name
--
-- @param tg_type Spatial type of geometry
--          1:[multi]point (puntal)
--          2:[multi]line  (lineal)
--          3:[multi]poly  (areal)
--          4:collection   (mixed)
--
-- @param layer_id Layer identifier
--
-- @param tg_objs Array of components
--
-- Return a topology.TopoGeometry object.
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.CreateTopoGeom(toponame varchar, tg_type integer, layer_id integer, tg_objs topology.TopoElementArray)
  RETURNS topology.TopoGeometry
AS
$$
DECLARE
  i integer;
  dims varchar;
  outerdims varchar;
  innerdims varchar;
  obj_type integer;
  obj_id integer;
  ret topology.TopoGeometry;
  rec RECORD;
  layertype integer;
  layerlevel integer;
  layerchild integer;
BEGIN

  IF tg_type < 1 OR tg_type > 4 THEN
    RAISE EXCEPTION 'Invalid TopoGeometry type % (must be in the range 1..4)', tg_type;
  END IF;

  -- Get topology id into return TopoGeometry
  SELECT id INTO ret.topology_id
    FROM topology.topology WHERE name = toponame;

  IF NOT FOUND THEN
    RAISE EXCEPTION 'Topology % does not exist', quote_literal(toponame);
  END IF;

  --
  -- Get layer info
  --
  layertype := NULL;
  FOR rec IN EXECUTE 'SELECT * FROM topology.layer'
       ' WHERE topology_id = ' || ret.topology_id
    || ' AND layer_id = ' || layer_id
  LOOP
    layertype = rec.feature_type;
    layerlevel = rec.level;
    layerchild = rec.child_id;
  END LOOP;

  -- Check for existence of given layer id
  IF layertype IS NULL THEN
    RAISE EXCEPTION 'No layer with id % is registered with topology %', layer_id, toponame;
  END IF;

  -- Verify compatibility between layer geometry type and
  -- TopoGeom requested geometry type
  IF layertype != 4 and layertype != tg_type THEN
    RAISE EXCEPTION 'A Layer of type % cannot contain a TopoGeometry of type %', layertype, tg_type;
  END IF;

  -- Set layer id and type in return object
  ret.layer_id = layer_id;
  ret.type = tg_type;

  --
  -- Get new TopoGeo id from sequence
  --
  FOR rec IN EXECUTE 'SELECT nextval(' ||
    quote_literal(
      quote_ident(toponame) || '.topogeo_s_' || layer_id
    ) || ')'
  LOOP
    ret.id = rec.nextval;
  END LOOP;

  -- Loop over outer dimension
  i = array_lower(tg_objs, 1);
  LOOP
    obj_id = tg_objs[i][1];
    obj_type = tg_objs[i][2];

    -- Elements of type 0 represent emptiness, just skip them
    IF obj_type = 0 THEN
      IF obj_id != 0 THEN
        RAISE EXCEPTION 'Malformed empty topo element {0,%} -- id must be 0 as well', obj_id;
      END IF;
    ELSE
      IF layerlevel = 0 THEN -- array specifies lower-level objects
        IF tg_type != 4 and tg_type != obj_type THEN
          RAISE EXCEPTION 'A TopoGeometry of type % cannot contain topology elements of type %', tg_type, obj_type;
        END IF;
      ELSE -- array specifies lower-level topogeometries
        IF obj_type != layerchild THEN
          RAISE EXCEPTION 'TopoGeom element layer do not match TopoGeom child layer';
        END IF;
        -- TODO: verify that the referred TopoGeometry really
        -- exists in the relation table ?
      END IF;

      --RAISE NOTICE 'obj:% type:% id:%', i, obj_type, obj_id;

      --
      -- Insert record into the Relation table
      --
      EXECUTE 'INSERT INTO '||quote_ident(toponame)
        || '.relation(topogeo_id, layer_id, '
           'element_id,element_type) '
           ' VALUES ('||ret.id
        ||','||ret.layer_id
        || ',' || obj_id || ',' || obj_type || ');';
    END IF;

    i = i+1;
    IF i > array_upper(tg_objs, 1) THEN
      EXIT;
    END IF;
  END LOOP;

  RETURN ret;

END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
--} CreateTopoGeom(toponame,topogeom_type, layer_id, TopoElementArray)

--{
-- CreateTopoGeom(topology_name, topogeom_type, layer_id) - creates the empty topogeom
-- Availability: 1.1.0
CREATE OR REPLACE FUNCTION topology.CreateTopoGeom(toponame varchar, tg_type integer, layer_id integer)
  RETURNS topology.TopoGeometry
AS
$$
  SELECT topology.CreateTopoGeom($1,$2,$3,'{{0,0}}');
$$ LANGUAGE 'sql' VOLATILE STRICT;
--} CreateTopoGeom(toponame, topogeom_type, layer_id)

--{
-- GetTopologyName(topology_id)
--
-- TODO: rewrite in SQL, as a wrapper around findTopology
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.GetTopologyName(topoid integer)
  RETURNS varchar
AS
$$
DECLARE
  ret varchar;
BEGIN
        SELECT name FROM topology.topology into ret
                WHERE id = topoid;
  RETURN ret;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;
--} GetTopologyName(topoid)

--{
-- GetTopologyId(toponame)
--
-- TODO: rewrite in SQL, as a wrapper around findTopology
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.GetTopologyId(toponame varchar)
  RETURNS integer
AS
$$
DECLARE
  ret integer;
BEGIN
  SELECT id INTO ret
    FROM topology.topology WHERE name = toponame;

  IF NOT FOUND THEN
    RAISE EXCEPTION 'Topology % does not exist', quote_literal(toponame);
  END IF;

  RETURN ret;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;
--} GetTopologyId(toponame)

--{
-- GetTopologySRID(toponame)
--
CREATE OR REPLACE FUNCTION topology.GetTopologySRID(toponame varchar)
  RETURNS integer
AS $$
  SELECT SRID FROM topology.topology WHERE name = $1;
$$ LANGUAGE 'sql' STABLE STRICT;
--} GetTopologySRID(toponame)

--{
-- GetTopoGeomElementArray(toponame, layer_id, topogeom_id)
-- GetTopoGeomElementArray(TopoGeometry)
--
-- Returns a set of element_id,element_type
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.GetTopoGeomElementArray(toponame varchar, layer_id integer, tgid integer)
  RETURNS topology.TopoElementArray
AS
$$
DECLARE
  rec RECORD;
  tg_objs varchar := '{';
  i integer;
  query text;
BEGIN

  query = 'SELECT * FROM topology.GetTopoGeomElements('
    || quote_literal(toponame) || ','
    || quote_literal(layer_id) || ','
    || quote_literal(tgid)
    || ') as obj ORDER BY obj';


  -- TODO: why not using array_agg here ?

  i = 1;
  FOR rec IN EXECUTE query
  LOOP
    IF i > 1 THEN
      tg_objs = tg_objs || ',';
    END IF;
    tg_objs = tg_objs || '{'
      || rec.obj[1] || ',' || rec.obj[2]
      || '}';
    i = i+1;
  END LOOP;

  tg_objs = tg_objs || '}';

  RETURN tg_objs;
END;
$$
LANGUAGE 'plpgsql' STABLE STRICT;

CREATE OR REPLACE FUNCTION topology.GetTopoGeomElementArray(tg topology.TopoGeometry)
  RETURNS topology.TopoElementArray
AS
$$
DECLARE
  toponame varchar;
BEGIN
  toponame = topology.GetTopologyName(tg.topology_id);
  RETURN topology.GetTopoGeomElementArray(toponame, tg.layer_id, tg.id);
END;
$$
LANGUAGE 'plpgsql' STABLE STRICT;

--} GetTopoGeomElementArray()

--{
-- GetTopoGeomElements(toponame, layer_id, topogeom_id)
-- GetTopoGeomElements(TopoGeometry)
--
-- Returns a set of element_id,element_type
--
CREATE OR REPLACE FUNCTION topology.GetTopoGeomElements(toponame varchar, layerid integer, tgid integer)
  RETURNS SETOF topology.TopoElement
AS
$$
DECLARE
  ret topology.TopoElement;
  rec RECORD;
  rec2 RECORD;
  query text;
  query2 text;
  lyr RECORD;
  ok bool;
  topoid INTEGER;
BEGIN

  -- Get topology id
  SELECT id INTO topoid
    FROM topology.topology WHERE name = toponame;
  IF NOT FOUND THEN
    RAISE EXCEPTION 'Topology % does not exist', quote_literal(toponame);
  END IF;

  -- Get layer info
  ok = false;
  FOR rec IN EXECUTE 'SELECT * FROM topology.layer '
       ' WHERE layer_id = $1 AND topology_id = $2'
       USING layerid, topoid
  LOOP
    lyr = rec;
    ok = true;
  END LOOP;

  IF NOT ok THEN
    RAISE EXCEPTION 'Layer % does not exist', layerid;
  END IF;

  query = 'SELECT abs(element_id) as element_id, element_type FROM '
    || quote_ident(toponame) || '.relation WHERE '
       ' layer_id = ' || layerid
    || ' AND topogeo_id = ' || quote_literal(tgid)
    || ' ORDER BY element_type, element_id';

  --RAISE NOTICE 'Query: %', query;

  FOR rec IN EXECUTE query
  LOOP
    IF lyr.level > 0 THEN
      query2 = 'SELECT * from topology.GetTopoGeomElements('
        || quote_literal(toponame) || ','
        || rec.element_type
        || ','
        || rec.element_id
        || ') as ret;';
      --RAISE NOTICE 'Query2: %', query2;
      FOR rec2 IN EXECUTE query2
      LOOP
        RETURN NEXT rec2.ret;
      END LOOP;
    ELSE
      ret = '{' || rec.element_id || ',' || rec.element_type || '}';
      RETURN NEXT ret;
    END IF;

  END LOOP;

  RETURN;
END;
$$
LANGUAGE 'plpgsql' STABLE STRICT;

CREATE OR REPLACE FUNCTION topology.GetTopoGeomElements(tg topology.TopoGeometry)
  RETURNS SETOF topology.TopoElement
AS
$$
DECLARE
  toponame varchar;
  rec RECORD;
BEGIN
  toponame = topology.GetTopologyName(tg.topology_id);
  FOR rec IN SELECT * FROM topology.GetTopoGeomElements(toponame,
    tg.layer_id,tg.id) as ret
  LOOP
    RETURN NEXT rec.ret;
  END LOOP;
  RETURN;
END;
$$
LANGUAGE 'plpgsql' STABLE STRICT;

--} GetTopoGeomElements()

--{
-- Geometry(TopoGeometry)
--
-- Construct a Geometry from a TopoGeometry.
--
-- }{
CREATE OR REPLACE FUNCTION topology.Geometry(topogeom topology.TopoGeometry)
  RETURNS Geometry
AS $BODY$
DECLARE
  toponame varchar;
  toposrid INT;
  geom geometry;
  elements_count INT;
  rec RECORD;
  plyr RECORD;
  clyr RECORD;
  sql TEXT;
BEGIN

  -- Get topology name
  SELECT name, srid FROM topology.topology
  WHERE id = topogeom.topology_id
  INTO toponame, toposrid;
  IF toponame IS NULL THEN
    RAISE EXCEPTION 'Invalid TopoGeometry (unexistent topology id %)', topogeom.topology_id;
  END IF;

  -- Get layer info
  SELECT * FROM topology.layer
    WHERE topology_id = topogeom.topology_id
    AND layer_id = topogeom.layer_id
    INTO plyr;
  IF plyr IS NULL THEN
    RAISE EXCEPTION 'Could not find TopoGeometry layer % in topology %', topogeom.layer_id, topogeom.topology_id;
  END IF;

  --
  -- If this feature layer is on any level > 0 we will
  -- compute the topological union of all child features
  -- in fact recursing.
  --
  IF plyr.level > 0 THEN -- {

    -- Get child layer info
    SELECT * FROM topology.layer WHERE layer_id = plyr.child_id
      AND topology_id = topogeom.topology_id
      INTO clyr;
    IF clyr IS NULL THEN
      RAISE EXCEPTION 'Invalid layer % in topology % (unexistent child layer %)', topogeom.layer_id, topogeom.topology_id, plyr.child_id;
    END IF;

    sql := 'SELECT st_multi(st_union(topology.Geometry('
      || quote_ident(clyr.feature_column)
      || '))) as geom FROM '
      || quote_ident(clyr.schema_name) || '.'
      || quote_ident(clyr.table_name)
      || ', ' || quote_ident(toponame) || '.relation pr'
         ' WHERE '
         ' pr.topogeo_id = ' || topogeom.id
      || ' AND '
         ' pr.layer_id = ' || topogeom.layer_id
      || ' AND '
         ' id('||quote_ident(clyr.feature_column)
      || ') = pr.element_id '
         ' AND '
         'layer_id('||quote_ident(clyr.feature_column)
      || ') = pr.element_type ';
    --RAISE DEBUG '%', query;
    EXECUTE sql INTO geom;

  ELSIF topogeom.type = 3 THEN -- [multi]polygon -- }{

    sql := format(
      $$
SELECT
  count(element_id),
  ST_Multi(
    ST_Union(
      topology.ST_GetFaceGeometry(
        %1$L,
        element_id
      )
    )
  ) as g
FROM %1$I.relation
WHERE topogeo_id = %2$L
AND layer_id = %3$L
AND element_type = 3
      $$,
      toponame,
      topogeom.id,
      topogeom.layer_id
    );
    EXECUTE sql INTO elements_count, geom;


  ELSIF topogeom.type = 2 THEN -- [multi]line -- }{

    sql := format(
      $$
SELECT
  st_multi(
    ST_LineMerge(
      ST_Collect(
        CASE
          WHEN r.element_id > 0 THEN
            e.geom
          ELSE
            ST_Reverse(e.geom)
        END
      )
    )
  ) as g
FROM %1$I.edge e, %1$I.relation r
WHERE r.topogeo_id = id($1)
AND r.layer_id = layer_id($1)
AND r.element_type = 2
AND abs(r.element_id) = e.edge_id
      $$,
      toponame
    );
    EXECUTE sql USING topogeom INTO geom;

  ELSIF topogeom.type = 1 THEN -- [multi]point -- }{

    sql :=
      'SELECT st_multi(st_union(n.geom)) as g FROM '
      || quote_ident(toponame) || '.node n, '
      || quote_ident(toponame) || '.relation r '
         ' WHERE r.topogeo_id = ' || topogeom.id
      || ' AND r.layer_id = ' || topogeom.layer_id
      || ' AND r.element_type = 1 '
         ' AND r.element_id = n.node_id';
    EXECUTE sql INTO geom;

  ELSIF topogeom.type = 4 THEN -- mixed collection -- }{

    sql := 'WITH areas AS ( SELECT ST_Union('
         'topology.ST_GetFaceGeometry('
      || quote_literal(toponame) || ','
      || 'element_id)) as g FROM '
      || quote_ident(toponame)
      || '.relation WHERE topogeo_id = '
      || topogeom.id || ' AND layer_id = '
      || topogeom.layer_id || ' AND element_type = 3), '
         'lines AS ( SELECT ST_LineMerge(ST_Collect(e.geom)) as g FROM '
      || quote_ident(toponame) || '.edge e, '
      || quote_ident(toponame) || '.relation r '
         ' WHERE r.topogeo_id = ' || topogeom.id
      || ' AND r.layer_id = ' || topogeom.layer_id
      || ' AND r.element_type = 2 '
         ' AND abs(r.element_id) = e.edge_id ), '
         ' points as ( SELECT st_union(n.geom) as g FROM '
      || quote_ident(toponame) || '.node n, '
      || quote_ident(toponame) || '.relation r '
         ' WHERE r.topogeo_id = ' || topogeom.id
      || ' AND r.layer_id = ' || topogeom.layer_id
      || ' AND r.element_type = 1 '
         ' AND r.element_id = n.node_id ), '
         ' un as ( SELECT g FROM areas UNION ALL SELECT g FROM lines '
         '          UNION ALL SELECT g FROM points ) '
         'SELECT ST_Multi(ST_Collect(g)) FROM un';
    EXECUTE sql INTO geom;

  ELSE -- }{

    RAISE EXCEPTION 'Invalid TopoGeometries (unknown type %)', topogeom.type;

  END IF; -- }

  IF geom IS NULL THEN
    IF topogeom.type = 3 THEN -- [multi]polygon
      geom := 'MULTIPOLYGON EMPTY';
    ELSIF topogeom.type = 2 THEN -- [multi]line
      geom := 'MULTILINESTRING EMPTY';
    ELSIF topogeom.type = 1 THEN -- [multi]point
      geom := 'MULTIPOINT EMPTY';
    ELSE
      geom := 'GEOMETRYCOLLECTION EMPTY';
    END IF;
    geom := ST_SetSRID(geom, toposrid);
  END IF;

  RETURN geom;
END
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;
--} Geometry(TopoGeometry)

-- 7.3+ explicit cast
CREATE CAST (topology.TopoGeometry AS Geometry) WITH FUNCTION topology.Geometry(topology.TopoGeometry) AS IMPLICIT;

--{
--  DropTopology(name)
--
-- Drops a topology schema getting rid of every dependent object.
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.DropTopology(atopology varchar)
RETURNS text
AS
$$
DECLARE
  topoid integer;
  rec RECORD;
  sql TEXT;
  toposchema REGNAMESPACE;
  deferred_constraints TEXT[];
BEGIN
  -- Get topology id
  SELECT id INTO topoid
    FROM topology.topology WHERE name = atopology;

  IF NOT FOUND THEN
    RAISE EXCEPTION 'Topology % does not exist', quote_literal(atopology);
  END IF;

  RAISE NOTICE 'Dropping all layers from topology % (%)',
    quote_literal(atopology), topoid;

  -- Drop all layers in the topology
  sql := 'SELECT * FROM topology.layer WHERE topology_id = $1';
  FOR rec IN EXECUTE sql USING topoid
  LOOP

    sql := format(
      'SELECT topology.DropTopoGeometryColumn(%L, %L, %L)',
      rec.schema_name, rec.table_name, rec.feature_column
    );
    EXECUTE sql;
  END LOOP;


  -- Delete record from topology.topology
  sql := 'DELETE FROM topology.topology WHERE id = $1';
  EXECUTE sql USING topoid;


  -- Drop the schema (if it exists)
  SELECT oid FROM pg_namespace WHERE text(nspname) = atopology
  INTO toposchema;

  IF toposchema IS NOT NULL THEN --{

    -- Give immediate execution to pending constraints
    -- in the topology schema.
    --
    -- See https://trac.osgeo.org/postgis/ticket/5115
    SELECT array_agg(format('%I.%I', atopology, conname))
    FROM pg_constraint c
    WHERE connamespace = toposchema AND condeferred
    INTO deferred_constraints;

    IF deferred_constraints IS NOT NULL THEN --{
      sql := format(
        'SET constraints %s IMMEDIATE',
        array_to_string(deferred_constraints, ',')
      );
      EXECUTE sql;
    END IF; --}

    sql := format('DROP SCHEMA %I CASCADE', atopology);
    EXECUTE sql;
  END IF; --}

  RETURN format('Topology %L dropped', atopology);
END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
--} DropTopology

-- Spatial predicates
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011-2012 Sandro Santilli <strk@kbt.io>
-- Copyright (C) 2005 Refractions Research Inc.
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- Author: Sandro Santilli <strk@kbt.io>
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--  Overloaded spatial predicates for TopoGeometry inputs
--
--   FUNCTION intersects(TopoGeometry, TopoGeometry)
--   FUNCTION equals(TopoGeometry, TopoGeometry)
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
-- intersects(TopoGeometry, TopoGeometry)
--
CREATE OR REPLACE FUNCTION topology.intersects(tg1 topology.TopoGeometry, tg2 topology.TopoGeometry)
  RETURNS bool
AS
$$
DECLARE
  tgbuf topology.TopoGeometry;
  rec RECORD;
  toponame varchar;
  query text;
BEGIN
  IF tg1.topology_id != tg2.topology_id THEN
    -- TODO: revert to ::geometry instead ?
    RAISE EXCEPTION 'Cannot compute intersection between TopoGeometries from different topologies';
  END IF;

  -- Order TopoGeometries so that tg1 has less-or-same
  -- dimensionality of tg1 (point,line,polygon,collection)
  IF tg1.type > tg2.type THEN
    tgbuf := tg2;
    tg2 := tg1;
    tg1 := tgbuf;
  END IF;

  --RAISE NOTICE 'tg1.id:% tg2.id:%', tg1.id, tg2.id;
  -- Geometry collection are not currently supported
  IF tg2.type = 4 THEN
    RAISE EXCEPTION 'GeometryCollection are not supported by intersects()';
  END IF;

        -- Get topology name
        SELECT name FROM topology.topology into toponame
                WHERE id = tg1.topology_id;

  -- Hierarchical TopoGeometries are not currently supported
  query = 'SELECT level FROM topology.layer'
    || ' WHERE '
    || ' topology_id = ' || tg1.topology_id
    || ' AND '
    || '( layer_id = ' || tg1.layer_id
    || ' OR layer_id = ' || tg2.layer_id
    || ' ) '
    || ' AND level > 0 ';

  --RAISE NOTICE '%', query;

  FOR rec IN EXECUTE query
  LOOP
    -- TODO: revert to ::geometry instead ?
    RAISE EXCEPTION 'Hierarchical TopoGeometries are not currently supported by intersects()';
  END LOOP;

  IF tg1.type = 1 THEN -- [multi]point

    IF tg2.type = 1 THEN -- point/point
  ---------------------------------------------------------
  --
  --  Two [multi]point features intersect if they share
  --  any Node
  --
  --
  --
      query =
        'SELECT a.topogeo_id FROM '
        || quote_ident(toponame) ||
        '.relation a, '
        || quote_ident(toponame) ||
        '.relation b '
        || 'WHERE a.layer_id = ' || tg1.layer_id
        || ' AND b.layer_id = ' || tg2.layer_id
        || ' AND a.topogeo_id = ' || tg1.id
        || ' AND b.topogeo_id = ' || tg2.id
        || ' AND a.element_id = b.element_id '
        || ' LIMIT 1';
      --RAISE NOTICE '%', query;
      FOR rec IN EXECUTE query
      LOOP
        RETURN TRUE; -- they share an element
      END LOOP;
      RETURN FALSE; -- no elements shared
  --
  ---------------------------------------------------------

    ELSIF tg2.type = 2 THEN -- point/line
  ---------------------------------------------------------
  --
  --  A [multi]point intersects a [multi]line if they share
  --  any Node.
  --
  --
  --
      query =
        'SELECT a.topogeo_id FROM '
        || quote_ident(toponame) ||
        '.relation a, '
        || quote_ident(toponame) ||
        '.relation b, '
        || quote_ident(toponame) ||
        '.edge_data e '
        || 'WHERE a.layer_id = ' || tg1.layer_id
        || ' AND b.layer_id = ' || tg2.layer_id
        || ' AND a.topogeo_id = ' || tg1.id
        || ' AND b.topogeo_id = ' || tg2.id
        || ' AND abs(b.element_id) = e.edge_id '
        || ' AND ( '
          || ' e.start_node = a.element_id '
          || ' OR '
          || ' e.end_node = a.element_id '
        || ' )'
        || ' LIMIT 1';
      --RAISE NOTICE '%', query;
      FOR rec IN EXECUTE query
      LOOP
        RETURN TRUE; -- they share an element
      END LOOP;
      RETURN FALSE; -- no elements shared
  --
  ---------------------------------------------------------

    ELSIF tg2.type = 3 THEN -- point/polygon
  ---------------------------------------------------------
  --
  --  A [multi]point intersects a [multi]polygon if any
  --  Node of the point is contained in any face of the
  --  polygon OR ( is end_node or start_node of any edge
  --  of any polygon face ).
  --
  --  We assume the Node-in-Face check is faster becasue
  --  there will be less Faces then Edges in any polygon.
  --
  --
  --
  --
      -- Check if any node is contained in a face
      query =
        'SELECT n.node_id as id FROM '
        || quote_ident(toponame) ||
        '.relation r1, '
        || quote_ident(toponame) ||
        '.relation r2, '
        || quote_ident(toponame) ||
        '.node n '
        || 'WHERE r1.layer_id = ' || tg1.layer_id
        || ' AND r2.layer_id = ' || tg2.layer_id
        || ' AND r1.topogeo_id = ' || tg1.id
        || ' AND r2.topogeo_id = ' || tg2.id
        || ' AND n.node_id = r1.element_id '
        || ' AND r2.element_id = n.containing_face '
        || ' LIMIT 1';
      --RAISE NOTICE '%', query;
      FOR rec IN EXECUTE query
      LOOP
        --RAISE NOTICE 'Node % in polygon face', rec.id;
        RETURN TRUE; -- one (or more) nodes are
                     -- contained in a polygon face
      END LOOP;

      -- Check if any node is start or end of any polygon
      -- face edge
      query =
        'SELECT n.node_id as nid, e.edge_id as eid '
        || ' FROM '
        || quote_ident(toponame) ||
        '.relation r1, '
        || quote_ident(toponame) ||
        '.relation r2, '
        || quote_ident(toponame) ||
        '.edge_data e, '
        || quote_ident(toponame) ||
        '.node n '
        || 'WHERE r1.layer_id = ' || tg1.layer_id
        || ' AND r2.layer_id = ' || tg2.layer_id
        || ' AND r1.topogeo_id = ' || tg1.id
        || ' AND r2.topogeo_id = ' || tg2.id
        || ' AND n.node_id = r1.element_id '
        || ' AND ( '
        || ' e.left_face = r2.element_id '
        || ' OR '
        || ' e.right_face = r2.element_id '
        || ' ) '
        || ' AND ( '
        || ' e.start_node = r1.element_id '
        || ' OR '
        || ' e.end_node = r1.element_id '
        || ' ) '
        || ' LIMIT 1';
      --RAISE NOTICE '%', query;
      FOR rec IN EXECUTE query
      LOOP
        --RAISE NOTICE 'Node % on edge % bound', rec.nid, rec.eid;
        RETURN TRUE; -- one node is start or end
                     -- of a face edge
      END LOOP;

      RETURN FALSE; -- no intersection
  --
  ---------------------------------------------------------

    ELSIF tg2.type = 4 THEN -- point/collection
      RAISE EXCEPTION 'Intersection point/collection not implemented yet';

    ELSE
      RAISE EXCEPTION 'Invalid TopoGeometry type %', tg2.type;
    END IF;

  ELSIF tg1.type = 2 THEN -- [multi]line
    IF tg2.type = 2 THEN -- line/line
  ---------------------------------------------------------
  --
  --  A [multi]line intersects a [multi]line if they share
  --  any Node.
  --
  --
  --
      query =
        'SELECT e1.start_node FROM '
        || quote_ident(toponame) ||
        '.relation r1, '
        || quote_ident(toponame) ||
        '.relation r2, '
        || quote_ident(toponame) ||
        '.edge_data e1, '
        || quote_ident(toponame) ||
        '.edge_data e2 '
        || 'WHERE r1.layer_id = ' || tg1.layer_id
        || ' AND r2.layer_id = ' || tg2.layer_id
        || ' AND r1.topogeo_id = ' || tg1.id
        || ' AND r2.topogeo_id = ' || tg2.id
        || ' AND abs(r1.element_id) = e1.edge_id '
        || ' AND abs(r2.element_id) = e2.edge_id '
        || ' AND ( '
        || ' e1.start_node = e2.start_node '
        || ' OR '
        || ' e1.start_node = e2.end_node '
        || ' OR '
        || ' e1.end_node = e2.start_node '
        || ' OR '
        || ' e1.end_node = e2.end_node '
        || ' )'
        || ' LIMIT 1';
      --RAISE NOTICE '%', query;
      FOR rec IN EXECUTE query
      LOOP
        RETURN TRUE; -- they share an element
      END LOOP;
      RETURN FALSE; -- no elements shared
  --
  ---------------------------------------------------------

    ELSIF tg2.type = 3 THEN -- line/polygon
  ---------------------------------------------------------
  --
  -- A [multi]line intersects a [multi]polygon if they share
  -- any Node (touch-only case), or if any line edge has any
  -- polygon face on the left or right (full-containment case
  -- + edge crossing case).
  --
  --
      -- E1 are line edges, E2 are polygon edges
      -- R1 are line relations.
      -- R2 are polygon relations.
      -- R2.element_id are FACE ids
      query =
        'SELECT e1.edge_id'
        || ' FROM '
        || quote_ident(toponame) ||
        '.relation r1, '
        || quote_ident(toponame) ||
        '.relation r2, '
        || quote_ident(toponame) ||
        '.edge_data e1, '
        || quote_ident(toponame) ||
        '.edge_data e2 '
        || 'WHERE r1.layer_id = ' || tg1.layer_id
        || ' AND r2.layer_id = ' || tg2.layer_id
        || ' AND r1.topogeo_id = ' || tg1.id
        || ' AND r2.topogeo_id = ' || tg2.id

        -- E1 are line edges
        || ' AND e1.edge_id = abs(r1.element_id) '

        -- E2 are face edges
        || ' AND ( e2.left_face = r2.element_id '
        || '   OR e2.right_face = r2.element_id ) '

        || ' AND ( '

        -- Check if E1 have left-or-right face
        -- being part of R2.element_id
        || ' e1.left_face = r2.element_id '
        || ' OR '
        || ' e1.right_face = r2.element_id '

        -- Check if E1 share start-or-end node
        -- with any E2.
        || ' OR '
        || ' e1.start_node = e2.start_node '
        || ' OR '
        || ' e1.start_node = e2.end_node '
        || ' OR '
        || ' e1.end_node = e2.start_node '
        || ' OR '
        || ' e1.end_node = e2.end_node '

        || ' ) '

        || ' LIMIT 1';
      --RAISE NOTICE '%', query;
      FOR rec IN EXECUTE query
      LOOP
        RETURN TRUE; -- either common node
                     -- or edge-in-face
      END LOOP;

      RETURN FALSE; -- no intersection
  --
  ---------------------------------------------------------

    ELSIF tg2.type = 4 THEN -- line/collection
      RAISE EXCEPTION 'Intersection line/collection not implemented yet';

    ELSE
      RAISE EXCEPTION 'Invalid TopoGeometry type %', tg2.type;
    END IF;

  ELSIF tg1.type = 3 THEN -- [multi]polygon

    IF tg2.type = 3 THEN -- polygon/polygon
  ---------------------------------------------------------
  --
  -- A [multi]polygon intersects a [multi]polygon if they share
  -- any Node (touch-only case), or if any face edge has any of the
  -- other polygon face on the left or right (full-containment case
  -- + edge crossing case).
  --
  --
      -- E1 are poly1 edges.
      -- E2 are poly2 edges
      -- R1 are poly1 relations.
      -- R2 are poly2 relations.
      -- R1.element_id are poly1 FACE ids
      -- R2.element_id are poly2 FACE ids
      query =
        'SELECT e1.edge_id'
        || ' FROM '
        || quote_ident(toponame) ||
        '.relation r1, '
        || quote_ident(toponame) ||
        '.relation r2, '
        || quote_ident(toponame) ||
        '.edge_data e1, '
        || quote_ident(toponame) ||
        '.edge_data e2 '
        || 'WHERE r1.layer_id = ' || tg1.layer_id
        || ' AND r2.layer_id = ' || tg2.layer_id
        || ' AND r1.topogeo_id = ' || tg1.id
        || ' AND r2.topogeo_id = ' || tg2.id

        -- E1 are poly1 edges
        || ' AND ( e1.left_face = r1.element_id '
        || '   OR e1.right_face = r1.element_id ) '

        -- E2 are poly2 edges
        || ' AND ( e2.left_face = r2.element_id '
        || '   OR e2.right_face = r2.element_id ) '

        || ' AND ( '

        -- Check if any edge from a polygon face
        -- has any of the other polygon face
        -- on the left or right
        || ' e1.left_face = r2.element_id '
        || ' OR '
        || ' e1.right_face = r2.element_id '
        || ' OR '
        || ' e2.left_face = r1.element_id '
        || ' OR '
        || ' e2.right_face = r1.element_id '

        -- Check if E1 share start-or-end node
        -- with any E2.
        || ' OR '
        || ' e1.start_node = e2.start_node '
        || ' OR '
        || ' e1.start_node = e2.end_node '
        || ' OR '
        || ' e1.end_node = e2.start_node '
        || ' OR '
        || ' e1.end_node = e2.end_node '

        || ' ) '

        || ' LIMIT 1';
      --RAISE NOTICE '%', query;
      FOR rec IN EXECUTE query
      LOOP
        RETURN TRUE; -- either common node
                     -- or edge-in-face
      END LOOP;

      RETURN FALSE; -- no intersection
  --
  ---------------------------------------------------------

    ELSIF tg2.type = 4 THEN -- polygon/collection
      RAISE EXCEPTION 'Intersection poly/collection not implemented yet';

    ELSE
      RAISE EXCEPTION 'Invalid TopoGeometry type %', tg2.type;
    END IF;

  ELSIF tg1.type = 4 THEN -- collection
    IF tg2.type = 4 THEN -- collection/collection
      RAISE EXCEPTION 'Intersection collection/collection not implemented yet';
    ELSE
      RAISE EXCEPTION 'Invalid TopoGeometry type %', tg2.type;
    END IF;

  ELSE
    RAISE EXCEPTION 'Invalid TopoGeometry type %', tg1.type;
  END IF;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;
--} intersects(TopoGeometry, TopoGeometry)

--{
--  equals(TopoGeometry, TopoGeometry)
--
CREATE OR REPLACE FUNCTION topology.equals(tg1 topology.TopoGeometry, tg2 topology.TopoGeometry)
  RETURNS bool
AS
$$
DECLARE
  rec RECORD;
  toponame varchar;
  query text;
BEGIN

  IF tg1.topology_id != tg2.topology_id THEN
    -- TODO: revert to ::geometry instead ?
    RAISE EXCEPTION 'Cannot compare TopoGeometries from different topologies';
  END IF;

  -- Not the same type, not equal
  IF tg1.type != tg2.type THEN
    RETURN FALSE;
  END IF;

  -- Geometry collection are not currently supported
  IF tg2.type = 4 THEN
    RAISE EXCEPTION 'GeometryCollection are not supported by equals()';
  END IF;

        -- Get topology name
        SELECT name FROM topology.topology into toponame
                WHERE id = tg1.topology_id;

  -- Two geometries are equal if they are composed by
  -- the same TopoElements
  FOR rec IN EXECUTE 'SELECT * FROM '
    || ' topology.GetTopoGeomElements('
    || quote_literal(toponame) || ', '
    || tg1.layer_id || ',' || tg1.id || ') '
    || ' EXCEPT SELECT * FROM '
    || ' topology.GetTopogeomElements('
    || quote_literal(toponame) || ', '
    || tg2.layer_id || ',' || tg2.id || ');'
  LOOP
    RETURN FALSE;
  END LOOP;

  FOR rec IN EXECUTE 'SELECT * FROM '
    || ' topology.GetTopoGeomElements('
    || quote_literal(toponame) || ', '
    || tg2.layer_id || ',' || tg2.id || ')'
    || ' EXCEPT SELECT * FROM '
    || ' topology.GetTopogeomElements('
    || quote_literal(toponame) || ', '
    || tg1.layer_id || ',' || tg1.id || '); '
  LOOP
    RETURN FALSE;
  END LOOP;
  RETURN TRUE;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;
--} equals(TopoGeometry, TopoGeometry)


--  Querying
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011 Andrea Peri <aperi2007@gmail.com>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-- {
--
-- Andrea Peri (19 Jan 2011) creation
-- Andrea Peri (14 Feb 2011) minor issues
--
-- getnodebypoint(atopology, point, tolerance)
--
-- Retrieve a Node ID given a POINT and a tolerance
-- tolerance = 0 mean exactly intersection
--
-- Returns the integer ID if there is a Node on the Point.
--
-- If there isn't any node in the Point, GetNodeByPoint return 0.
--
-- if near the point there are two or more nodes it throw an exception.
--
CREATE OR REPLACE FUNCTION topology.GetNodeByPoint(atopology varchar, apoint geometry, tol1 float8)
	RETURNS int AS
	'$libdir/postgis_topology-3', 'GetNodeByPoint'
	LANGUAGE 'c' STABLE STRICT;
--} GetNodeByPoint

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011 Andrea Peri <aperi2007@gmail.com>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--
-- Andrea Peri (19 Jan 2011) creation
-- Andrea Peri (14 Feb 2011) minor issues
--
-- GetEdgeByPoint(atopology, point, tol)
--
-- Retrieve an Edge ID given a POINT and a tolerance
-- tolerance = 0 mean exactly intersection
--
-- Returns return the integer ID if there is an edge on the Point.
--
-- When the Point is even a Node it raise an exception.
-- This case is testable with the GetNodeByPoint(atopology, apoint, tol)
--
-- If there isn't any edge in the Point, GetEdgeByPoint return 0.
--
-- if near the point there are two or more edges it throw an exception.
--
CREATE OR REPLACE FUNCTION topology.GetEdgeByPoint(atopology varchar, apoint geometry, tol1 float8)
	RETURNS int AS
	'$libdir/postgis_topology-3', 'GetEdgeByPoint'
	LANGUAGE 'c' STABLE STRICT;
--} GetEdgeByPoint
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2021 Sandro Santilli <strk@kbt.io>
-- Copyright (C) 2011 Andrea Peri <aperi2007@gmail.com>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--
-- GetFaceByPoint(atopology, point, tol)
--
-- Retrieve the ID of a face unambiguously referenced by the given POINT.
--
-- A face is unambiguously referenced by a point if a circle of radius
-- tol1 around the point intersects a single face of the topology.
--
-- In case of ambiguity (multiple faces intersect the circle)
-- this function will raise an exception
--
-- If the circle is completely outside any face, 0 is returned
-- (the circle is contained in the Universe Face).
--
-- Initial implementation by Andrea Peri (27 Feb 2011)
-- Rewrite by Sandro Santilli (29 Jun 2021)
--
-- }{
CREATE OR REPLACE FUNCTION topology.GetFaceByPoint(atopology varchar, apoint geometry, tol1 float8)
	RETURNS int AS
$BODY$ -- }{
DECLARE
  rec RECORD;
  sql TEXT;
  sideFaces INT[];
BEGIN

  -- Check if any edge intersects the query circle
  sql := format(
    $$
      WITH edges_in_circle AS (
        SELECT
          left_face,
          right_face
        FROM
          %1$I.edge
        WHERE
          ST_DWithin(geom, $1, $2)
      ), side_faces AS (
        SELECT left_face f FROM edges_in_circle
          UNION
        SELECT right_face FROM edges_in_circle
      )
      SELECT array_agg(f ORDER BY f) FROM side_faces;
    $$,
    atopology
  );
  EXECUTE sql
  USING apoint, tol1
  INTO sideFaces;

  RAISE DEBUG 'Side faces: %', sideFaces;

  IF array_upper(sideFaces, 1) = 1
  THEN
    -- Edges intersecting the circle
    -- have a single side-face, our circle
    -- is surely in that face
    --
    -- NOTE: this may also be the universe face
    --
    RETURN sideFaces[1];
  END IF;

  IF array_upper(sideFaces, 1) = 2
  THEN
    IF sideFaces[1] = 0
    THEN
      -- Edges intersecting the circle
      -- have a single real side-face,
      -- we'll consider our query to be in that face
      RETURN sideFaces[2];
    ELSE
      -- Edges have multiple real side-faces
      RAISE EXCEPTION 'Two or more faces found';
    END IF;
  END IF;

  IF array_upper(sideFaces, 1) > 2
  THEN
      RAISE EXCEPTION 'Two or more faces found';
  END IF;

  -- No edge intersects the circle, check for containment
  RETURN topology.GetFaceContainingPoint(atopology, apoint);
END;
$BODY$ LANGUAGE 'plpgsql' STABLE STRICT;
--} GetFaceByPoint

-- This function finds the id of a face whose geometry
-- would properly contain the given query point.
--
-- The implementation of this function relies on a properly built
-- topology, so results against an invalid topology are unpredictable.
--
-- If the point falls on the boundary of any face an exception will be
-- raised. Dangling edges will not be considered as face boundaries,
-- which is consistent with ST_GetFaceGeometry.
--
-- The following query should always return true
-- for non-boundary query points (:point)
--
--  ST_Contains(
--    ST_GetFaceGeometry(
--       :topo,
--       ST_GetFaceContainingPoint(:topo, :point)
--    ),
--    :point
--  )
--
--
-- Availability: 3.2.0
CREATE OR REPLACE FUNCTION topology.GetFaceContainingPoint(
  atopology text,
  apoint geometry
)
RETURNS INT AS
	'$libdir/postgis_topology-3', 'GetFaceContainingPoint'
LANGUAGE 'c' STABLE;

--  Populating
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2010-2012 Sandro Santilli <strk@kbt.io>
--
-- 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
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -



-- {
--  Compute the max size of the double-precision floating point grid
--  cell required to cover the given geometry
--
--
-- A pragmatic test conducted using algoritm shown here:
-- http://stackoverflow.com/questions/7408407/generate-next-largest-or-smallest-representable-floating-point-number-without-bi
-- showed the "tolerance" growing by an order of magnitude proportionally
-- with the order of magnitude of the input, starting with something like
-- 3.5527136788005009294e-15 for the starting value of 9.0
--
-- }{
CREATE OR REPLACE FUNCTION topology._st_mintolerance(ageom Geometry)
  RETURNS float8
AS $$
    SELECT 3.6 * power(10,  - ( 15 - log(coalesce(
      nullif(
        greatest(abs(ST_xmin($1)), abs(ST_ymin($1)),
                 abs(ST_xmax($1)), abs(ST_ymax($1))),
        0),
      1)) ));
$$ LANGUAGE 'sql' IMMUTABLE STRICT;
-- }

-- {
-- Get tolerance for a given topology
-- and if zero the minimum for the given geometry
--
-- }{
CREATE OR REPLACE FUNCTION topology._st_mintolerance(atopology varchar, ageom Geometry)
  RETURNS float8
AS $$
DECLARE
  ret FLOAT8;
BEGIN
  SELECT COALESCE(
    NULLIF(precision, 0),
    topology._st_mintolerance($2))
  FROM topology.topology
  WHERE name = $1 INTO ret;
  IF NOT FOUND THEN
    RAISE EXCEPTION
      'No topology with name "%" in topology.topology', atopology;
  END IF;
  return ret;
END;
$$ LANGUAGE 'plpgsql' STABLE STRICT;
-- }

--{
--
-- AddNode(atopology, point, allowEdgeSplitting, setContainingFace)
--
-- Add a node primitive to a topology and get its identifier.
-- Returns an existing node at the same location, if any.
--
-- 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.
--
-- Developed by Sandro Santilli <strk@kbt.io>
-- 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)"
--
-- }{
CREATE OR REPLACE FUNCTION topology.AddNode(atopology varchar, apoint geometry, allowEdgeSplitting boolean DEFAULT false, setContainingFace boolean DEFAULT false)
	RETURNS int
AS
$$
DECLARE
	nodeid int;
	rec RECORD;
  containing_face int;
BEGIN
	--
	-- Atopology and apoint are required
	--
	IF atopology IS NULL OR apoint IS NULL THEN
		RAISE EXCEPTION 'Invalid null argument';
	END IF;

	--
	-- Apoint must be a point
	--
	IF substring(geometrytype(apoint), 1, 5) != 'POINT'
	THEN
		RAISE EXCEPTION 'Node geometry must be a point';
	END IF;

	--
	-- Check if a coincident node already exists
	--
	-- We use index AND x/y equality
	--
	FOR rec IN EXECUTE 'SELECT node_id FROM '
		|| quote_ident(atopology) || '.node ' ||
		'WHERE geom && $1 AND ST_X(geom) = ST_X($1) AND ST_Y(geom) = ST_Y($1)'
    USING apoint
	LOOP
		RETURN  rec.node_id;
	END LOOP;

	--
	-- Check if any edge crosses this node
	-- (endpoints are fine)
	--
	FOR rec IN EXECUTE 'SELECT edge_id FROM '
		|| quote_ident(atopology) || '.edge '
		|| 'WHERE ST_DWithin($1, geom, 0) AND '
    || 'NOT ST_Equals($1, ST_StartPoint(geom)) AND '
    || 'NOT ST_Equals($1, ST_EndPoint(geom))'
    USING apoint
	LOOP
    IF allowEdgeSplitting THEN
      RETURN topology.ST_ModEdgeSplit(atopology, rec.edge_id, apoint);
    ELSE
		  RAISE EXCEPTION 'An edge crosses the given node.';
    END IF;
	END LOOP;

  IF setContainingFace THEN
    containing_face := topology.GetFaceByPoint(atopology, apoint, 0);
  ELSE
    containing_face := NULL;
  END IF;

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

	--
	-- Insert the new row
	--
	EXECUTE 'INSERT INTO ' || quote_ident(atopology)
		|| '.node(node_id, containing_face, geom)
		VALUES(' || nodeid || ',' || coalesce(containing_face::text, 'NULL')
    || ',$1)' USING apoint;

	RETURN nodeid;

END
$$
LANGUAGE 'plpgsql' VOLATILE;
--} AddNode

--{
--
-- AddEdge(atopology, line)
--
-- Add an edge primitive to a topology and get its identifier.
-- Edge endpoints will be added as nodes if missing.
-- Returns an existing edge at the same location, if any.
--
-- 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 as per next left/right edge.
-- Calling code is expected to do further linking.
--
-- Developed by Sandro Santilli <strk@kbt.io>
-- 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)"
--
CREATE OR REPLACE FUNCTION topology.AddEdge(atopology varchar, aline geometry)
	RETURNS int
AS
$$
DECLARE
	edgeid int;
	rec RECORD;
  ix geometry;
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;

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

	--
	-- Check if the edge crosses an existing node
	--
	FOR rec IN EXECUTE 'SELECT node_id FROM '
		|| quote_ident(atopology) || '.node '
		|| 'WHERE ST_Crosses($1, geom)'
    USING aline
	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.
	--
	--    F F 1
	--    F * *
	--    1 * 2
	--
	-- Example1: linestrings touching at one endpoint
	--    FF1 F00 102
	--    FF1 F** 1*2 <-- our match
	--
	-- Example2: linestrings touching at both endpoints
	--    FF1 F0F 1F2
	--    FF1 F** 1*2 <-- our match
	--
	FOR rec IN EXECUTE 'SELECT edge_id, geom, ST_Relate($1, geom, 2) as im FROM '
		|| quote_ident(atopology) || '.edge WHERE $1 && geom'
    USING aline
	LOOP

	  IF ST_RelateMatch(rec.im, 'FF1F**1*2') THEN
	    CONTINUE; -- no interior intersection
	  END IF;

	  -- Reuse an EQUAL edge (be it closed or not)
	  IF ST_RelateMatch(rec.im, '1FFF*FFF2') THEN
	      RETURN rec.edge_id;
	  END IF;

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

    -- Find a point on the intersection which
    -- is NOT an endpoint of "aline"
    IF ST_Dimension(ix) = 0
    THEN
      WITH SharedBounds AS (
        (
          SELECT ST_Force2D(ST_StartPoint(rec.geom)) g
          UNION
          SELECT ST_Force2D(ST_EndPoint(rec.geom))
        )
        INTERSECT
        (
          SELECT ST_Force2D(ST_StartPoint(aline))
          UNION
          SELECT ST_Force2D(ST_EndPoint(aline))
        )
      )
      SELECT d.geom
      FROM ST_DumpPoints(ix) d
      WHERE ST_Force2D(geom) NOT IN ( SELECT g FROM SharedBounds )
      ORDER BY d.path
      LIMIT 1
      INTO STRICT ix;
    ELSE
      -- for linear intersection we pick
      -- an internal point.
      ix := ST_PointOnSurface(ix);
    END IF;

    RAISE EXCEPTION
      'Edge intersects (not on endpoints)'
        ' with existing edge % at or near point %',
      rec.edge_id,
      ST_AsText(ix);

	END LOOP;

	--
	-- Get new edge id from sequence
	--
	FOR rec IN EXECUTE 'SELECT nextval(' ||
		quote_literal(
			quote_ident(atopology) || '.edge_data_edge_id_seq'
		) || ')'
	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($1)), '

		-- end_node
		|| 'topology.addNode('
		|| quote_literal(atopology)
		|| ', ST_EndPoint($1)), '

		-- next_left_edge
		|| -edgeid ||','

		-- next_right_edge
		|| edgeid ||','

		-- left_face
		|| '0,'

		-- right_face
		|| '0,'

		-- geom
		|| '$1)'
    USING aline;

	RETURN edgeid;

END
$$
LANGUAGE 'plpgsql' VOLATILE;
--} AddEdge

--{
--
-- AddFace(atopology, poly, [<force_new>=true])
--
-- Add a face primitive to a topology and get its identifier.
-- Returns an existing face at the same location, if any, unless
-- true is passed as the force_new argument
--
-- For a newly added face, its edges will be appropriately
-- linked (marked as left-face or right-face), and any contained
-- edges and nodes would also be marked as such.
--
-- When forcing re-registration of an existing face, no action will be
-- taken to deal with the face being substituted. Which means
-- a record about the old face and any record in the relation table
-- referencing the existing face will remain untouched, effectively
-- leaving the topology in a possibly invalid state.
-- It is up to the caller to deal with that.
--
-- 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.
--
-- Developed by Sandro Santilli <strk@kbt.io>
-- 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)"
--
CREATE OR REPLACE FUNCTION topology.AddFace(atopology varchar, apoly geometry, force_new boolean DEFAULT FALSE)
	RETURNS int
AS
$$
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
$$
LANGUAGE 'plpgsql' VOLATILE;
--} AddFace

-- ----------------------------------------------------------------------------
--
-- Functions to incrementally populate a topology
--
-- ----------------------------------------------------------------------------

--{
--  TopoGeo_AddPoint(toponame, pointgeom, tolerance)
--
--  Add a Point into a topology, with an optional tolerance
--
CREATE OR REPLACE FUNCTION topology.TopoGeo_AddPoint(atopology varchar, apoint geometry, tolerance float8 DEFAULT 0)
	RETURNS int AS
	'$libdir/postgis_topology-3', 'TopoGeo_AddPoint'
  LANGUAGE 'c' VOLATILE;
--} TopoGeo_AddPoint

--{
--  TopoGeo_addLinestring(toponame, linegeom, tolerance)
--
--  Add a LineString into a topology
--
-- }{
CREATE OR REPLACE FUNCTION topology.TopoGeo_addLinestring(atopology varchar, aline geometry, tolerance float8 DEFAULT 0)
	RETURNS SETOF int AS
	'$libdir/postgis_topology-3', 'TopoGeo_AddLinestring'
  LANGUAGE 'c' VOLATILE;
--} TopoGeo_addLinestring

--{
--  TopoGeo_AddPolygon(toponame, polygeom, tolerance)
--
--  Add a Polygon into a topology
--
-- }{
CREATE OR REPLACE FUNCTION topology.TopoGeo_AddPolygon(atopology varchar, apoly geometry, tolerance float8 DEFAULT 0)
	RETURNS SETOF int AS
	'$libdir/postgis_topology-3', 'TopoGeo_AddPolygon'
  LANGUAGE 'c' VOLATILE;
--} TopoGeo_AddPolygon

--{
--  TopoGeo_AddGeometry(toponame, geom, tolerance)
--
--  Add a Geometry into a topology
--
CREATE OR REPLACE FUNCTION topology.TopoGeo_AddGeometry(atopology varchar, ageom geometry, tolerance float8 DEFAULT 0)
	RETURNS void AS
$$
DECLARE
BEGIN
	RAISE EXCEPTION 'TopoGeo_AddGeometry not implemented yet';
END
$$
LANGUAGE 'plpgsql';
--} TopoGeo_AddGeometry
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011 Sandro Santilli <strk@kbt.io>
--
-- 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 polygonize topology edges
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--
-- int Polygonize(toponame)
--
-- TODO: allow restricting polygonization to a bounding box
--
-- }{
CREATE OR REPLACE FUNCTION topology.polygonize(toponame varchar)
  RETURNS text
AS
$$
DECLARE
  sql text;
  rec RECORD;
  faces int;
BEGIN

  sql := 'SELECT (st_dump(st_polygonize(geom))).geom from '
         || quote_ident(toponame) || '.edge_data';

  faces = 0;
  FOR rec in EXECUTE sql LOOP
    BEGIN
      PERFORM topology.AddFace(toponame, rec.geom);
      faces = faces + 1;
    EXCEPTION
      WHEN OTHERS THEN
        RAISE WARNING 'Error registering face % (%)', rec.geom, SQLERRM;
    END;
  END LOOP;
  RETURN faces || ' faces registered';
END
$$
LANGUAGE 'plpgsql';
--} Polygonize(toponame)

--  TopoElement
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2010, 2011 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- TopoElement management functions
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- Developed by Sandro Santilli <strk@kbt.io>
-- 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)"
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--
-- TopoElementArray TopoElementArray_append(<TopoElement>)
--
-- Append a TopoElement to a TopoElementArray
--
CREATE OR REPLACE FUNCTION topology.TopoElementArray_append(topology.TopoElementArray, topology.TopoElement)
	RETURNS topology.TopoElementArray
AS
$$
	SELECT CASE
		WHEN $1 IS NULL THEN
			topology.TopoElementArray('{' || $2::text || '}')
		ELSE
			topology.TopoElementArray($1::int[][]||$2::int[])
		END;
$$
LANGUAGE 'sql' IMMUTABLE;
--} TopoElementArray_append

--{
--
-- TopoElementArray TopoElementArray_agg(<setof TopoElement>)
--
-- Aggregates a set of TopoElement values into a TopoElementArray
--
-- Availability: 2.0.0
DROP AGGREGATE IF EXISTS topology.TopoElementArray_agg(topology.TopoElement);
CREATE AGGREGATE topology.TopoElementArray_agg(
	sfunc = topology.TopoElementArray_append,
	basetype = topology.TopoElement,
	stype = topology.TopoElementArray
	);
--} TopoElementArray_agg

--  TopoGeometry
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-- {
--  Override geometrytype() for topogeometry objects
--
--  Note: For performance reasons, this function always assumes
--        TopoGeometry are of the MULTI type. This may not always
--        be the case if you convert the TopoGeometry to an actual
--        Geometry.
--
-- }{
CREATE OR REPLACE FUNCTION topology.GeometryType(tg topology.TopoGeometry)
	RETURNS text
AS
$$
	SELECT CASE
		WHEN type($1) = 1 THEN 'MULTIPOINT'
		WHEN type($1) = 2 THEN 'MULTILINESTRING'
		WHEN type($1) = 3 THEN 'MULTIPOLYGON'
		WHEN type($1) = 4 THEN 'GEOMETRYCOLLECTION'
		ELSE 'UNEXPECTED'
		END;
$$
LANGUAGE 'sql' STABLE STRICT;
-- }

-- {
--  Override st_geometrytype() for topogeometry objects
--
--  Note: For performance reasons, this function always assumes
--        TopoGeometry are of the MULTI type. This may not always
--        be the case if you convert the TopoGeometry to an actual
--        Geometry.
--
-- }{
CREATE OR REPLACE FUNCTION topology.ST_GeometryType(tg topology.TopoGeometry)
	RETURNS text
AS
$$
	SELECT CASE
		WHEN type($1) = 1 THEN 'ST_MultiPoint'
		WHEN type($1) = 2 THEN 'ST_MultiLinestring'
		WHEN type($1) = 3 THEN 'ST_MultiPolygon'
		WHEN type($1) = 4 THEN 'ST_GeometryCollection'
		ELSE 'ST_Unexpected'
		END;
$$
LANGUAGE 'sql' STABLE STRICT;
-- }
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2021-2022 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-- {
-- Override ST_Srid() for topogeometry objects
--
-- Availability: 3.2.0
--
-- }{
CREATE OR REPLACE FUNCTION topology.ST_Srid(tg topology.TopoGeometry)
	RETURNS INT
AS
$$
	SELECT srid FROM topology.topology
  WHERE id = topology_id(tg);
$$
LANGUAGE 'sql' VOLATILE STRICT;
-- }

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2012 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-- {
--  Clear the contents of a TopoGeometry
--
-- }{
CREATE OR REPLACE FUNCTION topology.clearTopoGeom(tg topology.TopoGeometry)
  RETURNS topology.TopoGeometry
AS
$$
DECLARE
  topology_info RECORD;
  sql TEXT;
BEGIN

  -- Get topology information
  SELECT id, name FROM topology.topology
    INTO topology_info
    WHERE id = topology_id(tg);
  IF NOT FOUND THEN
      RAISE EXCEPTION 'No topology with id "%" in topology.topology', topology_id(tg);
  END IF;

  -- Clear the TopoGeometry contents
  sql := 'DELETE FROM ' || quote_ident(topology_info.name)
        || '.relation WHERE layer_id = '
        || layer_id(tg)
        || ' AND topogeo_id = '
        || id(tg);
  EXECUTE sql;

  RETURN tg;

END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-- }

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2012 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-- {
--  Get a simplified geometry version from a TopoGeometry
--
--  Performs Douglas Peucker algorithm on each edge composing
--  the given TopoGeometry
--
-- }{
CREATE OR REPLACE FUNCTION topology.ST_Simplify(tg topology.TopoGeometry, tolerance float8)
  RETURNS geometry
AS
$$
DECLARE
  topology_info RECORD;
  layer_info RECORD;
  child_layer_info RECORD;
  geom geometry;
  sql TEXT;
BEGIN

  -- Get topology information
  SELECT id, name FROM topology.topology
    INTO topology_info
    WHERE id = tg.topology_id;
  IF NOT FOUND THEN
      RAISE EXCEPTION 'No topology with id "%" in topology.topology', tg.topology_id;
  END IF;

  -- Get layer info
  SELECT * FROM topology.layer
    WHERE topology_id = tg.topology_id
    AND layer_id = tg.layer_id
    INTO layer_info;
  IF NOT FOUND THEN
    RAISE EXCEPTION 'Could not find TopoGeometry layer % in topology %', tg.layer_id, tg.topology_id;
  END IF;

  --
  -- If this feature layer is on any level > 0 we will
  -- compute the topological union of all simplified child
  -- features in fact recursing.
  --
  IF layer_info.level > 0 THEN -- {

    -- Get child layer info
    SELECT * FROM topology.layer WHERE layer_id = layer_info.child_id
      AND topology_id = tg.topology_id
      INTO child_layer_info;
    IF NOT FOUND THEN
      RAISE EXCEPTION 'Invalid layer % in topology % (unexistent child layer %)', tg.layer_id, tg.topology_id, layer_info.child_id;
    END IF;

    sql := 'SELECT st_multi(st_union(topology.ST_Simplify('
      || quote_ident(child_layer_info.feature_column)
      || ',' || tolerance || '))) as geom FROM '
      || quote_ident(child_layer_info.schema_name) || '.'
      || quote_ident(child_layer_info.table_name)
      || ', ' || quote_ident(topology_info.name) || '.relation pr'
      || ' WHERE '
      || ' pr.topogeo_id = ' || tg.id
      || ' AND '
      || ' pr.layer_id = ' || tg.layer_id
      || ' AND '
      || ' id('||quote_ident(child_layer_info.feature_column)
      || ') = pr.element_id '
      || ' AND '
      || 'layer_id('||quote_ident(child_layer_info.feature_column)
      || ') = pr.element_type ';
    RAISE DEBUG '%', sql;
    EXECUTE sql INTO geom;

  ELSIF tg.type = 3 THEN -- [multi]polygon -- }{

    -- TODO: use ST_GetFaceEdges
    -- TODO: is st_unaryunion needed?
    sql := 'SELECT st_multi(st_unaryunion(ST_BuildArea(ST_Node(ST_Collect(ST_Simplify(geom, '
      || tolerance || ')))))) as geom FROM '
      || quote_ident(topology_info.name)
      || '.edge_data e, '
      || quote_ident(topology_info.name)
      || '.relation r WHERE ( e.left_face = r.element_id'
      || ' OR e.right_face = r.element_id )'
      || ' AND r.topogeo_id = ' || tg.id
      || ' AND r.layer_id = ' || tg.layer_id
      || ' AND element_type = 3 ';
    RAISE DEBUG '%', sql;
    EXECUTE sql INTO geom;

  ELSIF tg.type = 2 THEN -- [multi]line -- }{

    sql :=
      'SELECT st_multi(ST_LineMerge(ST_Node(ST_Collect(ST_Simplify(e.geom,'
      || tolerance || '))))) as g FROM '
      || quote_ident(topology_info.name) || '.edge e, '
      || quote_ident(topology_info.name) || '.relation r '
      || ' WHERE r.topogeo_id = ' || tg.id
      || ' AND r.layer_id = ' || tg.layer_id
      || ' AND r.element_type = 2 '
      || ' AND abs(r.element_id) = e.edge_id';
    EXECUTE sql INTO geom;

  ELSIF tg.type = 1 THEN -- [multi]point -- }{

    -- Can't simplify points...
    geom := topology.Geometry(tg);

  ELSIF tg.type = 4 THEN -- mixed collection -- }{

   sql := 'WITH areas AS ( '
      || 'SELECT st_multi(st_union(ST_BuildArea(ST_Node(ST_Collect(ST_Simplify(geom, '
      || tolerance || ')))) as geom FROM '
      || quote_ident(topology_info.name)
      || '.edge_data e, '
      || quote_ident(topology_info.name)
      || '.relation r WHERE ( e.left_face = r.element_id'
      || ' OR e.right_face = r.element_id )'
      || ' AND r.topogeo_id = ' || tg.id
      || ' AND r.layer_id = ' || tg.layer_id
      || ' AND element_type = 3 ), '
      || 'lines AS ( '
      || 'SELECT st_multi(ST_LineMerge(ST_Collect(ST_Simplify(e.geom,'
      || tolerance || ')))) as g FROM '
      || quote_ident(topology_info.name) || '.edge e, '
      || quote_ident(topology_info.name) || '.relation r '
      || ' WHERE r.topogeo_id = ' || tg.id
      || ' AND r.layer_id = ' || tg.layer_id
      || ' AND r.element_type = 2 '
      || ' AND abs(r.element_id) = e.edge_id ), '
      || ' points as ( SELECT st_union(n.geom) as g FROM '
      || quote_ident(topology_info.name) || '.node n, '
      || quote_ident(topology_info.name) || '.relation r '
      || ' WHERE r.topogeo_id = ' || tg.id
      || ' AND r.layer_id = ' || tg.layer_id
      || ' AND r.element_type = 1 '
      || ' AND r.element_id = n.node_id ), '
      || ' un as ( SELECT g FROM areas UNION ALL SELECT g FROM lines '
      || '          UNION ALL SELECT g FROM points ) '
      || 'SELECT ST_Multi(ST_Collect(g)) FROM un';
    EXECUTE sql INTO geom;

  ELSE -- }{

    RAISE EXCEPTION 'Invalid TopoGeometries (unknown type %)', tg.type;

  END IF; -- }

  RETURN geom;

END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-- }

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011-2012 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -



-- {
--  Convert a simple geometry to a topologically-defined one
--
--  See http://trac.osgeo.org/postgis/ticket/1017
--
-- }{
CREATE OR REPLACE FUNCTION topology.toTopoGeom(ageom Geometry, atopology varchar, alayer int, atolerance float8 DEFAULT 0)
  RETURNS topology.TopoGeometry
AS
$$
DECLARE
  layer_info RECORD;
  topology_info RECORD;
  tg topology.TopoGeometry;
  typ TEXT;
BEGIN

  -- Get topology information
  BEGIN
    SELECT *
    FROM topology.topology
      INTO STRICT topology_info WHERE name = atopology;
  EXCEPTION
    WHEN NO_DATA_FOUND THEN
      RAISE EXCEPTION 'No topology with name "%" in topology.topology',
        atopology;
  END;

  -- Get layer information
  BEGIN
    SELECT *, CASE
      WHEN feature_type = 1 THEN 'puntal'
      WHEN feature_type = 2 THEN 'lineal'
      WHEN feature_type = 3 THEN 'areal'
      WHEN feature_type = 4 THEN 'mixed'
      ELSE 'unexpected_'||feature_type
      END as typename
    FROM topology.layer l
      INTO STRICT layer_info
      WHERE l.layer_id = alayer
      AND l.topology_id = topology_info.id;
  EXCEPTION
    WHEN NO_DATA_FOUND THEN
      RAISE EXCEPTION 'No layer with id "%" in topology "%"',
        alayer, atopology;
  END;

  -- Can't convert to a hierarchical topogeometry
  IF layer_info.level > 0 THEN
      RAISE EXCEPTION 'Layer "%" of topology "%" is hierarchical, cannot convert to it.',
        alayer, atopology;
  END IF;

  --
  -- Check type compatibility and create empty TopoGeometry
  -- 1:puntal, 2:lineal, 3:areal, 4:collection
  --
  typ = geometrytype(ageom);
  IF typ = 'GEOMETRYCOLLECTION' THEN
    --  A collection can only go collection layer
    IF layer_info.feature_type != 4 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold a collection feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg := topology.CreateTopoGeom(atopology, 4, alayer);
  ELSIF typ = 'POINT' OR typ = 'MULTIPOINT' THEN -- puntal
    --  A point can go in puntal or collection layer
    IF layer_info.feature_type != 4 and layer_info.feature_type != 1 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold a puntal feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg := topology.CreateTopoGeom(atopology, 1, alayer);
  ELSIF typ = 'LINESTRING' or typ = 'MULTILINESTRING' THEN -- lineal
    --  A line can go in lineal or collection layer
    IF layer_info.feature_type != 4 and layer_info.feature_type != 2 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold a lineal feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg := topology.CreateTopoGeom(atopology, 2, alayer);
  ELSIF typ = 'POLYGON' OR typ = 'MULTIPOLYGON' THEN -- areal
    --  An area can go in areal or collection layer
    IF layer_info.feature_type != 4 and layer_info.feature_type != 3 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold an areal feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg := topology.CreateTopoGeom(atopology, 3, alayer);
  ELSE
      -- Should never happen
      RAISE EXCEPTION
        'Unsupported feature type %', typ;
  END IF;

  tg := topology.toTopoGeom(ageom, tg, atolerance);

  RETURN tg;

END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-- }

-- {
--  Convert a simple geometry to a topologically-defined one
--  adding its components to a pre-existing TopoGeometry
--
-- }{
CREATE OR REPLACE FUNCTION topology.toTopoGeom(ageom Geometry, tg topology.TopoGeometry, atolerance float8 DEFAULT 0)
  RETURNS topology.TopoGeometry
AS
$$
DECLARE
  layer_info RECORD;
  topology_info RECORD;
  rec RECORD;
  rec2 RECORD;
  elem TEXT;
  elems TEXT[];
  sql TEXT;
  typ TEXT;
  tolerance FLOAT8;
  alayer INT;
  atopology TEXT;
BEGIN


  -- Get topology information
  SELECT id, name FROM topology.topology
    INTO topology_info
    WHERE id = topology_id(tg);
  IF NOT FOUND THEN
    RAISE EXCEPTION 'No topology with id "%" in topology.topology',
                    topology_id(tg);
  END IF;

  alayer := layer_id(tg);
  atopology := topology_info.name;

  -- Get tolerance, if 0 was given
  tolerance := COALESCE( NULLIF(atolerance, 0), topology._st_mintolerance(topology_info.name, ageom) );

  -- Get layer information
  BEGIN
    SELECT *, CASE
      WHEN feature_type = 1 THEN 'puntal'
      WHEN feature_type = 2 THEN 'lineal'
      WHEN feature_type = 3 THEN 'areal'
      WHEN feature_type = 4 THEN 'mixed'
      ELSE 'unexpected_'||feature_type
      END as typename
    FROM topology.layer l
      INTO STRICT layer_info
      WHERE l.layer_id = layer_id(tg)
      AND l.topology_id = topology_info.id;
  EXCEPTION
    WHEN NO_DATA_FOUND THEN
      RAISE EXCEPTION 'No layer with id "%" in topology "%"',
        alayer, atopology;
  END;

  -- Can't convert to a hierarchical topogeometry
  IF layer_info.level > 0 THEN
      RAISE EXCEPTION 'Layer "%" of topology "%" is hierarchical, cannot convert a simple geometry to it.',
        alayer, atopology;
  END IF;

  --
  -- Check type compatibility and set TopoGeometry type
  -- 1:puntal, 2:lineal, 3:areal, 4:collection
  --
  typ = geometrytype(ageom);
  IF typ = 'GEOMETRYCOLLECTION' THEN
    --  A collection can only go to collection layer
    IF layer_info.feature_type != 4 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold a collection feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg.type := 4;
  ELSIF typ = 'POINT' OR typ = 'MULTIPOINT' THEN -- puntal
    --  A point can go in puntal or collection layer
    IF layer_info.feature_type != 4 and layer_info.feature_type != 1 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold a puntal feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg.type := CASE WHEN tg.type = 1 THEN 1 ELSE 4 END;
  ELSIF typ = 'LINESTRING' or typ = 'MULTILINESTRING' THEN -- lineal
    --  A line can go in lineal or collection layer
    IF layer_info.feature_type != 4 and layer_info.feature_type != 2 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold a lineal feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg.type := CASE WHEN tg.type = 2 THEN 2 ELSE 4 END;
  ELSIF typ = 'POLYGON' OR typ = 'MULTIPOLYGON' THEN -- areal
    --  An area can go in areal or collection layer
    IF layer_info.feature_type != 4 and layer_info.feature_type != 3 THEN
      RAISE EXCEPTION
        'Layer "%" of topology "%" is %, cannot hold an areal feature.',
        layer_info.layer_id, topology_info.name, layer_info.typename;
    END IF;
    tg.type := CASE WHEN tg.type = 3 THEN 3 ELSE 4 END;
  ELSE
      -- Should never happen
      RAISE EXCEPTION
        'Unexpected feature dimension %', ST_Dimension(ageom);
  END IF;

  -- Now that we have an empty topogeometry, we loop over distinct components
  -- and add them to the definition of it. We add them as soon
  -- as possible so that each element can further edit the
  -- definition by splitting
  FOR rec IN SELECT id(tg), alayer as lyr,
    geom, ST_Dimension(gd.geom) as dims
    FROM ST_Dump(ageom) AS gd
    WHERE NOT ST_IsEmpty(gd.geom)
  LOOP
    -- NOTE: Switched from using case to this because of PG 10 behavior change
    -- Using a UNION ALL only one will be processed because of the WHERE
    -- Since the WHERE clause will be processed first
    FOR rec2 IN SELECT primitive
          FROM
            (
              SELECT topology.topogeo_addPoint(atopology, rec.geom, tolerance)
                WHERE rec.dims = 0
              UNION ALL
              SELECT topology.topogeo_addLineString(atopology, rec.geom, tolerance)
                WHERE rec.dims = 1
              UNION ALL
              SELECT topology.topogeo_addPolygon(atopology, rec.geom, tolerance)
                WHERE rec.dims = 2
            ) AS f(primitive)
    LOOP
      elem := ARRAY[rec.dims+1, rec2.primitive]::text;
      IF elems @> ARRAY[elem] THEN
      ELSE
        elems := elems || elem;
        -- TODO: consider use a single INSERT statement for the whole thing
        sql := 'INSERT INTO ' || quote_ident(atopology)
            || '.relation(topogeo_id, layer_id, element_type, element_id) VALUES ('
            || rec.id || ',' || rec.lyr || ',' || rec.dims+1
            || ',' || rec2.primitive || ')'
            -- NOTE: we're avoiding duplicated rows here
            || ' EXCEPT SELECT ' || rec.id || ', ' || rec.lyr
            || ', element_type, element_id FROM '
            || quote_ident(topology_info.name)
            || '.relation WHERE layer_id = ' || rec.lyr
            || ' AND topogeo_id = ' || rec.id;
        EXECUTE sql;
      END IF;
    END LOOP;
  END LOOP;

  RETURN tg;

END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-- }
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2015 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -



-- {
--  Add an element to a TopoGeometry definition
--
-- }{
CREATE OR REPLACE FUNCTION topology.TopoGeom_addElement(tg topology.TopoGeometry, el topology.TopoElement)
  RETURNS topology.TopoGeometry
AS
$$
DECLARE
  toponame TEXT;
  sql TEXT;
BEGIN

  -- Get topology name
  BEGIN
    SELECT name
    FROM topology.topology
      INTO STRICT toponame WHERE id = topology_id(tg);
  EXCEPTION
    WHEN NO_DATA_FOUND THEN
      RAISE EXCEPTION 'No topology with name "%" in topology.topology',
        atopology;
  END;

  -- Insert new element
  sql := format('INSERT INTO %s.relation'
         '(topogeo_id,layer_id,element_id,element_type)'
         ' VALUES($1,$2,$3,$4)', quote_ident(toponame));
  BEGIN
    EXECUTE sql USING id(tg), layer_id(tg), el[1], el[2];
  EXCEPTION
    WHEN unique_violation THEN
      -- already present, let go
    WHEN OTHERS THEN
      RAISE EXCEPTION 'Got % (%)', SQLERRM, SQLSTATE;
  END;

  RETURN tg;

END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-- }

-- {
--  Remove an element from a TopoGeometry definition
--
-- }{
CREATE OR REPLACE FUNCTION topology.TopoGeom_remElement(tg topology.TopoGeometry, el topology.TopoElement)
  RETURNS topology.TopoGeometry
AS
$$
DECLARE
  toponame TEXT;
  sql TEXT;
BEGIN

  -- Get topology name
  BEGIN
    SELECT name
    FROM topology.topology
      INTO STRICT toponame WHERE id = topology_id(tg);
  EXCEPTION
    WHEN NO_DATA_FOUND THEN
      RAISE EXCEPTION 'No topology with name "%" in topology.topology',
        atopology;
  END;

  -- Delete the element
  sql := format('DELETE FROM %s.relation WHERE '
         'topogeo_id = $1 AND layer_id = $2 AND '
         'element_id = $3 AND element_type = $4',
         quote_ident(toponame));
  EXECUTE sql USING id(tg), layer_id(tg), el[1], el[2];

  RETURN tg;

END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-- }


-- {
-- Add the component of a TopoGeometry to the definition of
-- another TopoGeometry.
--
-- The two TopoGeometry objects need to be defined on the *same*
-- topology and need to be compatible (both simple or built over
-- the same child layer, and the target TopoGeometry needs to allow
-- for holding components of the type found in the source TopoGeometry)
--
-- }{
CREATE OR REPLACE FUNCTION topology.TopoGeom_addTopoGeom(tgt topology.TopoGeometry, src topology.TopoGeometry)
  RETURNS topology.TopoGeometry
AS
$BODY$
DECLARE
  sql TEXT;
  topo topology.topology;
  srcElementTypes int[];
  srcLayer topology.layer;
  tgtLayer topology.layer;
  maxElemType int;
BEGIN

  -- Get topology information
  topo := topology.FindTopology(topology_id(src));


  IF topology_id(src) != topology_id(tgt) THEN
    RAISE EXCEPTION 'Source and target TopoGeometry objects need be defined on the same topology';
  END IF;


  SELECT * FROM topology.layer
  WHERE topology_id = topo.id
    AND layer_id = layer_id(src)
  INTO srcLayer;

  SELECT * FROM topology.layer
  WHERE topology_id = topo.id
    AND layer_id = layer_id(tgt)
  INTO tgtLayer;

  -- Check simple/hierarchical compatibility
  IF srcLayer.child_id IS NULL THEN
    IF srcLayer.child_id IS NOT NULL THEN
      RAISE EXCEPTION 'Cannot add components of hierarchical TopoGeometry to a non-hierarchical TopoGeometry';
    END IF;
  ELSIF tgtLayer.child_id IS NULL THEN
      RAISE EXCEPTION 'Cannot add components of non-hierarchical TopoGeometry to a hierarchical TopoGeometry';
  ELSIF tgtLayer.child_id != srcLayer.childId THEN
      RAISE EXCEPTION 'Cannot add components of hierarchical TopoGeometry to a hierarchical TopoGeometry based on different layer';
  END IF;

  -- Add every element of the source TopoGeometry to
  -- the definition of the target TopoGeometry
  sql := format($$
WITH inserted AS (
  INSERT INTO %1$I.relation(
    topogeo_id,
    layer_id,
    element_id,
    element_type
  )
  SELECT %2$s, %3$s, element_id, element_type
  FROM %1$I.relation
  WHERE topogeo_id = %4$L
  AND layer_id = %5$L
  EXCEPT
  SELECT %2$s, %3$s, element_id, element_type
  FROM %1$I.relation
  WHERE topogeo_id = %2$L
  AND layer_id = %3$L
  RETURNING element_type
)
SELECT array_agg(DISTINCT element_type) FROM inserted
    $$,
    topo.name,      -- %1
    id(tgt),        -- %2
    layer_id(tgt),  -- %3
    id(src),        -- %4
    layer_id(src)   -- %5
  );

  RAISE DEBUG 'SQL: %', sql;

  EXECUTE sql INTO srcElementTypes;

  -- TODO: Check layer's feature_type compatibility ?
  -- or let the relationTrigger take care of it ?
--  IF tgtLayer.feature_type != 4 THEN -- 'mixed' typed target can accept anything
--    IF srcLayer.feature_type != tgtLayer.feature_type THEN
--    END IF;
--  END IF;

  RAISE DEBUG 'Target type: %', type(tgt);
  RAISE DEBUG 'Detected source element types: %', srcElementTypes;

  -- Check if target TopoGeometry type needs be changed
  IF type(tgt) != 4 -- collection TopoGeometry accept anything
  THEN
    IF array_upper(srcElementTypes, 1) > 1
    OR srcElementTypes[1] != tgt.type
    THEN
      -- source is mixed-typed or typed differently from
      -- target, so we turn target type to collection
      RAISE DEBUG 'Changing target element type to collection';
      tgt.type = 4;
    END IF;
  END IF;




  RETURN tgt;
END
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;
-- }

--  Exports
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2010, 2011 Sandro Santilli <strk@kbt.io>
--
-- 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 for topology GML output
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- Developed by Sandro Santilli <strk@kbt.io>
-- 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)"
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--
-- INTERNAL FUNCTION
-- text _AsGMLNode(id, point, nsprefix, precision, options, idprefix, gmlver)
--
-- }{
CREATE OR REPLACE FUNCTION topology._AsGMLNode(id int, point geometry, nsprefix_in text, prec int, options int, idprefix text, gmlver int)
  RETURNS text
AS
$$
DECLARE
  nsprefix text;
  gml text;
BEGIN

  nsprefix := 'gml:';
  IF NOT nsprefix_in IS NULL THEN
    IF nsprefix_in = '' THEN
      nsprefix = nsprefix_in;
    ELSE
      nsprefix = nsprefix_in || ':';
    END IF;
  END IF;

  gml := '<' || nsprefix || 'Node ' || nsprefix
    || 'id="' || idprefix || 'N' || id || '"';
  IF point IS NOT NULL THEN
    gml = gml || '>'
              || '<' || nsprefix || 'pointProperty>'
              || ST_AsGML(gmlver, point, prec, options, nsprefix_in)
              || '</' || nsprefix || 'pointProperty>'
              || '</' || nsprefix || 'Node>';
  ELSE
    gml = gml || '/>';
  END IF;
  RETURN gml;
END
$$
LANGUAGE 'plpgsql' IMMUTABLE;
--} _AsGMLNode(id, point, nsprefix, precision, options, idprefix, gmlVersion)

--{
--
-- INTERNAL FUNCTION
-- text _AsGMLEdge(edge_id, start_node, end_node, line, visitedTable,
--                 nsprefix, precision, options, idprefix, gmlVersion)
--
-- }{
CREATE OR REPLACE FUNCTION topology._AsGMLEdge(edge_id int, start_node int,end_node int, line geometry, visitedTable regclass, nsprefix_in text, prec int, options int, idprefix text, gmlver int)
  RETURNS text
AS
$$
DECLARE
  visited bool;
  nsprefix text;
  gml text;
BEGIN

  nsprefix := 'gml:';
  IF nsprefix_in IS NOT NULL THEN
    IF nsprefix_in = '' THEN
      nsprefix = nsprefix_in;
    ELSE
      nsprefix = nsprefix_in || ':';
    END IF;
  END IF;

  gml := '<' || nsprefix || 'Edge ' || nsprefix
    || 'id="' || idprefix || 'E' || edge_id || '">';

  -- Start node
  gml = gml || '<' || nsprefix || 'directedNode orientation="-"';
  -- Do visited bookkeeping if visitedTable was given
  visited = NULL;
  IF visitedTable IS NOT NULL THEN
    EXECUTE 'SELECT true FROM '
            || visitedTable::text
            || ' WHERE element_type = 1 AND element_id = '
            || start_node LIMIT 1 INTO visited;
    IF visited IS NOT NULL THEN
      gml = gml || ' xlink:href="#' || idprefix || 'N' || start_node || '" />';
    ELSE
      -- Mark as visited
      EXECUTE 'INSERT INTO ' || visitedTable::text
        || '(element_type, element_id) VALUES (1, '
        || start_node || ')';
    END IF;
  END IF;
  IF visited IS NULL THEN
    gml = gml || '>';
    gml = gml || topology._AsGMLNode(start_node, NULL, nsprefix_in,
                                     prec, options, idprefix, gmlver);
    gml = gml || '</' || nsprefix || 'directedNode>';
  END IF;

  -- End node
  gml = gml || '<' || nsprefix || 'directedNode';
  -- Do visited bookkeeping if visitedTable was given
  visited = NULL;
  IF visitedTable IS NOT NULL THEN
    EXECUTE 'SELECT true FROM '
            || visitedTable::text
            || ' WHERE element_type = 1 AND element_id = '
            || end_node LIMIT 1 INTO visited;
    IF visited IS NOT NULL THEN
      gml = gml || ' xlink:href="#' || idprefix || 'N' || end_node || '" />';
    ELSE
      -- Mark as visited
      EXECUTE 'INSERT INTO ' || visitedTable::text
        || '(element_type, element_id) VALUES (1, '
        || end_node || ')';
    END IF;
  END IF;
  IF visited IS NULL THEN
    gml = gml || '>';
    gml = gml || topology._AsGMLNode(end_node, NULL, nsprefix_in,
                                     prec, options, idprefix, gmlver);
    gml = gml || '</' || nsprefix || 'directedNode>';
  END IF;

  IF line IS NOT NULL THEN
    gml = gml || '<' || nsprefix || 'curveProperty>'
              || ST_AsGML(gmlver, line, prec, options, nsprefix_in)
              || '</' || nsprefix || 'curveProperty>';
  END IF;

  gml = gml || '</' || nsprefix || 'Edge>';

  RETURN gml;
END
$$
LANGUAGE 'plpgsql' VOLATILE; -- writes into visitedTable
--} _AsGMLEdge(id, start_node, end_node, line, visitedTable, nsprefix, precision, options, idprefix, gmlver)

--{
--
-- INTERNAL FUNCTION
-- text _AsGMLFace(toponame, face_id, visitedTable,
--                 nsprefix, precision, options, idprefix, gmlVersion)
--
-- }{
CREATE OR REPLACE FUNCTION topology._AsGMLFace(toponame text, face_id int, visitedTable regclass, nsprefix_in text, prec int, options int, idprefix text, gmlver int)
  RETURNS text
AS
$$
DECLARE
  visited bool;
  nsprefix text;
  gml text;
  rec RECORD;
  rec2 RECORD;
  bounds geometry;
BEGIN

  nsprefix := 'gml:';
  IF nsprefix_in IS NOT NULL THEN
    IF nsprefix_in = '' THEN
      nsprefix = nsprefix_in;
    ELSE
      nsprefix = nsprefix_in || ':';
    END IF;
  END IF;

  gml := '<' || nsprefix || 'Face ' || nsprefix
    || 'id="' || idprefix || 'F' || face_id || '">';

  -- Construct the face geometry, then for each polygon:
  FOR rec IN SELECT (ST_DumpRings((ST_Dump(ST_ForceRHR(
    topology.ST_GetFaceGeometry(toponame, face_id)))).geom)).geom
  LOOP

      -- Contents of a directed face are the list of edges
      -- that cover the specific ring
      bounds = ST_Boundary(rec.geom);

      FOR rec2 IN EXECUTE
        'SELECT e.*, ST_LineLocatePoint($1'
        || ', ST_LineInterpolatePoint(e.geom, 0.2)) as pos'
        || ', ST_LineLocatePoint($1'
        || ', ST_LineInterpolatePoint(e.geom, 0.8)) as pos2 FROM '
        || quote_ident(toponame)
        || '.edge e WHERE ( e.left_face = $2'
        || ' OR e.right_face = $2'
        || ') AND ST_Covers($1'
        || ', e.geom) ORDER BY pos'
        USING bounds, face_id
      LOOP

        gml = gml || '<' || nsprefix || 'directedEdge';

        -- if this edge goes in same direction to the
        --       ring bounds, make it with negative orientation
        IF rec2.pos2 > rec2.pos THEN -- edge goes in same direction
          gml = gml || ' orientation="-"';
        END IF;

        -- Do visited bookkeeping if visitedTable was given
        IF visitedTable IS NOT NULL THEN

          EXECUTE 'SELECT true FROM '
            || visitedTable::text
            || ' WHERE element_type = 2 AND element_id = '
            || rec2.edge_id LIMIT 1 INTO visited;
          IF visited THEN
            -- Use xlink:href if visited
            gml = gml || ' xlink:href="#' || idprefix || 'E'
                      || rec2.edge_id || '" />';
            CONTINUE;
          ELSE
            -- Mark as visited otherwise
            EXECUTE 'INSERT INTO ' || visitedTable::text
              || '(element_type, element_id) VALUES (2, '
              || rec2.edge_id || ')';
          END IF;

        END IF;

        gml = gml || '>';

        gml = gml || topology._AsGMLEdge(rec2.edge_id, rec2.start_node,
                                        rec2.end_node, rec2.geom,
                                        visitedTable, nsprefix_in,
                                        prec, options, idprefix, gmlver);
        gml = gml || '</' || nsprefix || 'directedEdge>';

      END LOOP;
    END LOOP;

  gml = gml || '</' || nsprefix || 'Face>';

  RETURN gml;
END
$$
LANGUAGE 'plpgsql' VOLATILE; -- writes into visited table
--} _AsGMLFace(toponame, id, visitedTable, nsprefix, precision, options, idprefix, gmlver)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry, nsprefix, precision, options, visitedTable, idprefix, gmlver)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry, nsprefix_in text, precision_in int, options_in int, visitedTable regclass, idprefix text, gmlver int)
  RETURNS text
AS
$$
DECLARE
  nsprefix text;
  precision int;
  options int;
  visited bool;
  toponame text;
  gml text;
  sql text;
  rec RECORD;
  rec2 RECORD;
BEGIN

  nsprefix := 'gml:';
  IF nsprefix_in IS NOT NULL THEN
    IF nsprefix_in = '' THEN
      nsprefix = nsprefix_in;
    ELSE
      nsprefix = nsprefix_in || ':';
    END IF;
  END IF;

  precision := 15;
  IF precision_in IS NOT NULL THEN
    precision = precision_in;
  END IF;

  options := 1;
  IF options_in IS NOT NULL THEN
    options = options_in;
  END IF;

  -- Get topology name (for subsequent queries)
  SELECT name FROM topology.topology into toponame
              WHERE id = tg.topology_id;

  -- Puntual TopoGeometry
  IF tg.type = 1 THEN
    gml = '<' || nsprefix || 'TopoPoint>';
    -- For each defining node, print a directedNode
    FOR rec IN  EXECUTE 'SELECT r.element_id, n.geom from '
      || quote_ident(toponame) || '.relation r LEFT JOIN '
      || quote_ident(toponame) || '.node n ON (r.element_id = n.node_id)'
      || ' WHERE r.layer_id = ' || tg.layer_id
      || ' AND r.topogeo_id = ' || tg.id
    LOOP
      gml = gml || '<' || nsprefix || 'directedNode';
      -- Do visited bookkeeping if visitedTable was given
      IF visitedTable IS NOT NULL THEN
        EXECUTE 'SELECT true FROM '
                || visitedTable::text
                || ' WHERE element_type = 1 AND element_id = '
                || rec.element_id LIMIT 1 INTO visited;
        IF visited IS NOT NULL THEN
          gml = gml || ' xlink:href="#' || idprefix || 'N' || rec.element_id || '" />';
          CONTINUE;
        ELSE
          -- Mark as visited
          EXECUTE 'INSERT INTO ' || visitedTable::text
            || '(element_type, element_id) VALUES (1, '
            || rec.element_id || ')';
        END IF;
      END IF;
      gml = gml || '>';
      gml = gml || topology._AsGMLNode(rec.element_id, rec.geom, nsprefix_in, precision, options, idprefix, gmlver);
      gml = gml || '</' || nsprefix || 'directedNode>';
    END LOOP;
    gml = gml || '</' || nsprefix || 'TopoPoint>';
    RETURN gml;

  ELSIF tg.type = 2 THEN -- lineal
    gml = '<' || nsprefix || 'TopoCurve>';

    FOR rec IN SELECT (ST_Dump(topology.Geometry(tg))).geom
    LOOP
      FOR rec2 IN EXECUTE
        'SELECT e.*, ST_LineLocatePoint($1'
        || ', ST_LineInterpolatePoint(e.geom, 0.2)) as pos'
        || ', ST_LineLocatePoint($1'
        || ', ST_LineInterpolatePoint(e.geom, 0.8)) as pos2 FROM '
        || quote_ident(toponame)
        || '.edge e WHERE ST_Covers($1'
        || ', e.geom) ORDER BY pos'
        -- TODO: add relation to the conditional, to reduce load ?
        USING rec.geom
      LOOP

        gml = gml || '<' || nsprefix || 'directedEdge';

        -- if this edge goes in opposite direction to the
        --       line, make it with negative orientation
        IF rec2.pos2 < rec2.pos THEN -- edge goes in opposite direction
          gml = gml || ' orientation="-"';
        END IF;

        -- Do visited bookkeeping if visitedTable was given
        IF visitedTable IS NOT NULL THEN

          EXECUTE 'SELECT true FROM '
            || visitedTable::text
            || ' WHERE element_type = 2 AND element_id = '
            || rec2.edge_id LIMIT 1 INTO visited;
          IF visited THEN
            -- Use xlink:href if visited
            gml = gml || ' xlink:href="#' || idprefix || 'E' || rec2.edge_id || '" />';
            CONTINUE;
          ELSE
            -- Mark as visited otherwise
            EXECUTE 'INSERT INTO ' || visitedTable::text
              || '(element_type, element_id) VALUES (2, '
              || rec2.edge_id || ')';
          END IF;

        END IF;

        gml = gml || '>';

        gml = gml || topology._AsGMLEdge(rec2.edge_id,
                                        rec2.start_node,
                                        rec2.end_node, rec2.geom,
                                        visitedTable,
                                        nsprefix_in, precision,
                                        options, idprefix, gmlver);

        gml = gml || '</' || nsprefix || 'directedEdge>';
      END LOOP;
    END LOOP;

    gml = gml || '</' || nsprefix || 'TopoCurve>';
    return gml;

  ELSIF tg.type = 3 THEN -- areal
    gml = '<' || nsprefix || 'TopoSurface>';

    -- For each defining face, print a directedFace
    FOR rec IN  EXECUTE 'SELECT f.face_id from '
      || quote_ident(toponame) || '.relation r LEFT JOIN '
      || quote_ident(toponame) || '.face f ON (r.element_id = f.face_id)'
      || ' WHERE r.layer_id = ' || tg.layer_id
      || ' AND r.topogeo_id = ' || tg.id
    LOOP
      gml = gml || '<' || nsprefix || 'directedFace';
      -- Do visited bookkeeping if visitedTable was given
      IF visitedTable IS NOT NULL THEN
        EXECUTE 'SELECT true FROM '
                || visitedTable::text
                || ' WHERE element_type = 3 AND element_id = '
                || rec.face_id LIMIT 1 INTO visited;
        IF visited IS NOT NULL THEN
          gml = gml || ' xlink:href="#' || idprefix || 'F' || rec.face_id || '" />';
          CONTINUE;
        ELSE
          -- Mark as visited
          EXECUTE 'INSERT INTO ' || visitedTable::text
            || '(element_type, element_id) VALUES (3, '
            || rec.face_id || ')';
        END IF;
      END IF;
      gml = gml || '>';
      gml = gml || topology._AsGMLFace(toponame, rec.face_id, visitedTable,
                                       nsprefix_in, precision,
                                       options, idprefix, gmlver);
      gml = gml || '</' || nsprefix || 'directedFace>';
    END LOOP;
    gml = gml || '</' || nsprefix || 'TopoSurface>';
    RETURN gml;

  ELSIF tg.type = 4 THEN -- collection
    RAISE EXCEPTION 'Collection TopoGeometries are not supported by AsGML';

  END IF;

  RETURN gml;

END
$$
LANGUAGE 'plpgsql' VOLATILE; -- writes into visited table
--} AsGML(TopoGeometry, nsprefix, precision, options, visitedTable, idprefix, gmlver)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry, nsprefix, precision, options, visitedTable,
--            idprefix)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry,nsprefix text, prec int, options int, visitedTable regclass, idprefix text)
  RETURNS text
AS
$$
 SELECT topology.AsGML($1, $2, $3, $4, $5, $6, 3);
$$
LANGUAGE 'sql' VOLATILE; -- writes into visited table
--} AsGML(TopoGeometry, nsprefix, precision, options, visitedTable, idprefix)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry, nsprefix, precision, options, visitedTable)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry, nsprefix text, prec int, options int, vis regclass)
  RETURNS text AS
$$
 SELECT topology.AsGML($1, $2, $3, $4, $5, '');
$$ LANGUAGE 'sql' VOLATILE; -- writes into visited table
-- } AsGML(TopoGeometry, nsprefix, precision, options)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry, nsprefix, precision, options)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry, nsprefix text, prec int, opts int)
  RETURNS text AS
$$
 SELECT topology.AsGML($1, $2, $3, $4, NULL);
$$ LANGUAGE 'sql' STABLE; -- does NOT write into visited table
-- } AsGML(TopoGeometry, nsprefix, precision, options)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry, nsprefix)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry, nsprefix text)
  RETURNS text AS
$$
 SELECT topology.AsGML($1, $2, 15, 1, NULL);
$$ LANGUAGE 'sql' STABLE; -- does NOT write into visited table
-- } AsGML(TopoGeometry, nsprefix)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry, visited_table)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry, visitedTable regclass)
  RETURNS text AS
$$
 SELECT topology.AsGML($1, 'gml', 15, 1, $2);
$$ LANGUAGE 'sql' VOLATILE; -- writes into visited table
-- } AsGML(TopoGeometry, visited_table)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry, visited_table, nsprefix)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry, visitedTable regclass, nsprefix text)
  RETURNS text AS
$$
 SELECT topology.AsGML($1, $3, 15, 1, $2);
$$ LANGUAGE 'sql' VOLATILE; -- writes into visited table
-- } AsGML(TopoGeometry, visited_table, nsprefix)

--{
--
-- API FUNCTION
--
-- text AsGML(TopoGeometry)
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsGML(tg topology.TopoGeometry)
  RETURNS text AS
$$
 SELECT topology.AsGML($1, 'gml');
$$ LANGUAGE 'sql' STABLE; -- does NOT write into visited table
-- } AsGML(TopoGeometry)

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2013-2020 Sandro Santilli <strk@kbt.io>
--
-- 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 for TopoJSON export
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -



--{
--
-- API FUNCTION
--
-- text AsTopoJSON(TopoGeometry, edgeMapTable)
--
-- Format specification here:
-- http://github.com/mbostock/topojson-specification/blob/master/README.md
--
-- }{
CREATE OR REPLACE FUNCTION topology.AsTopoJSON(tg topology.TopoGeometry, edgeMapTable regclass)
  RETURNS text AS
$$
DECLARE
  toponame text;
  json text;
  sql text;
  rec RECORD;
  rec2 RECORD;
  side int;
  arcid int;
  arcs int[];
  ringtxt TEXT[];
  comptxt TEXT[];
  edges_found BOOLEAN;
  old_search_path TEXT;
  all_faces int[];
  faces int[];
  shell_faces int[];
  visited_edges int[];
  looking_for_holes BOOLEAN;
BEGIN

  IF tg IS NULL THEN
    RETURN NULL;
  END IF;

  -- Get topology name (for subsequent queries)
  SELECT name FROM topology.topology into toponame
              WHERE id = tg.topology_id;

  -- TODO: implement scale ?

  -- Puntal TopoGeometry, simply delegate to AsGeoJSON
  IF tg.type = 1 THEN
    json := ST_AsGeoJSON(topology.Geometry(tg));
    return json;
  ELSIF tg.type = 2 THEN -- lineal

    FOR rec IN SELECT (ST_Dump(topology.Geometry(tg))).geom
    LOOP -- {

      sql := 'SELECT e.*, ST_LineLocatePoint($1'
            || ', ST_LineInterpolatePoint(e.geom, 0.2)) as pos'
            || ', ST_LineLocatePoint($1'
            || ', ST_LineInterpolatePoint(e.geom, 0.8)) as pos2 FROM '
            || quote_ident(toponame)
            || '.edge e WHERE ST_Covers($1'
            || ', e.geom) ORDER BY pos';
            -- TODO: add relation to the conditional, to reduce load ?
      FOR rec2 IN EXECUTE sql USING rec.geom
      LOOP -- {

        IF edgeMapTable IS NOT NULL THEN
          sql := 'SELECT arc_id-1 FROM ' || edgeMapTable::text || ' WHERE edge_id = $1';
          EXECUTE sql INTO arcid USING rec2.edge_id;
          IF arcid IS NULL THEN
            EXECUTE 'INSERT INTO ' || edgeMapTable::text
              || '(edge_id) VALUES ($1) RETURNING arc_id-1'
            INTO arcid USING rec2.edge_id;
          END IF;
        ELSE
          arcid := rec2.edge_id;
        END IF;

        -- edge goes in opposite direction
        IF rec2.pos2 < rec2.pos THEN
          arcid := -(arcid+1);
        END IF;

        arcs := arcs || arcid;

      END LOOP; -- }

      comptxt := comptxt || ( '[' || array_to_string(arcs, ',') || ']' );
      arcs := NULL;

    END LOOP; -- }

    json := '{ "type": "MultiLineString", "arcs": [' || array_to_string(comptxt,',') || ']}';

    return json;

  ELSIF tg.type = 3 THEN -- areal

    json := '{ "type": "MultiPolygon", "arcs": [';

    EXECUTE 'SHOW search_path' INTO old_search_path;
    EXECUTE 'SET search_path TO ' || quote_ident(toponame) || ',' || old_search_path;

    SELECT array_agg(id) as f
    FROM ( SELECT (topology.GetTopoGeomElements(tg))[1] as id ) as f
    INTO all_faces;


    visited_edges := ARRAY[]::int[];
    faces := all_faces;
    looking_for_holes := false;
    shell_faces := ARRAY[]::int[];

    CREATE TEMP TABLE _postgis_topology_astopojson_tmp_edges
    ON COMMIT DROP
    AS
    SELECT
         ROW_NUMBER() OVER (
            ORDER BY
              ST_XMin(e.geom),
              ST_YMin(e.geom),
              edge_id
         ) leftmost_index,
         e.edge_id,
         e.left_face,
         e.right_face,
         e.next_right_edge,
         e.next_left_edge
    FROM edge e
    WHERE
         ( e.left_face = ANY ( all_faces ) OR
           e.right_face = ANY ( all_faces ) )
    ;
    CREATE INDEX on _postgis_topology_astopojson_tmp_edges (edge_id);

    LOOP -- { until all edges were visited

      arcs := NULL;
      edges_found := false;


      FOR rec in -- {
WITH RECURSIVE
_edges AS (
  SELECT
     *,
     left_face = ANY ( faces ) as lf,
     right_face = ANY ( faces ) as rf
  FROM
    _postgis_topology_astopojson_tmp_edges
),
_leftmost_non_dangling_edge AS (
  SELECT e.edge_id
    FROM _edges e WHERE e.lf != e.rf
  ORDER BY
    leftmost_index
  LIMIT 1
),
_edgepath AS (
  SELECT
    CASE
      WHEN e.lf THEN lme.edge_id
      ELSE -lme.edge_id
    END as signed_edge_id,
    false as back,

    e.lf = e.rf as dangling,
    e.left_face, e.right_face,
    e.lf, e.rf,
    e.next_right_edge, e.next_left_edge

  FROM _edges e, _leftmost_non_dangling_edge lme
  WHERE e.edge_id = abs(lme.edge_id)
    UNION
  SELECT
    CASE
      WHEN p.dangling AND NOT p.back THEN -p.signed_edge_id
      WHEN p.signed_edge_id < 0 THEN p.next_right_edge
      ELSE p.next_left_edge
    END, -- signed_edge_id
    CASE
      WHEN p.dangling AND NOT p.back THEN true
      ELSE false
    END, -- back

    e.lf = e.rf, -- dangling
    e.left_face, e.right_face,
    e.lf, e.rf,
    e.next_right_edge, e.next_left_edge

  FROM _edges e, _edgepath p
  WHERE
    e.edge_id = CASE
      WHEN p.dangling AND NOT p.back THEN abs(p.signed_edge_id)
      WHEN p.signed_edge_id < 0 THEN abs(p.next_right_edge)
      ELSE abs(p.next_left_edge)
    END
)
SELECT abs(signed_edge_id) as edge_id, signed_edge_id, dangling,
        lf, rf, left_face, right_face
FROM _edgepath
      -- }

      LOOP  -- { over recursive query


        IF rec.left_face = ANY (all_faces) AND NOT rec.left_face = ANY (shell_faces) THEN
          shell_faces := shell_faces || rec.left_face;
        END IF;

        IF rec.right_face = ANY (all_faces) AND NOT rec.right_face = ANY (shell_faces) THEN
          shell_faces := shell_faces || rec.right_face;
        END IF;

        visited_edges := visited_edges || rec.edge_id;

        edges_found := true;

        -- TODO: drop ?
        IF rec.dangling THEN
          CONTINUE;
        END IF;

        IF rec.left_face = ANY (all_faces) AND rec.right_face = ANY (all_faces) THEN
          CONTINUE;
        END IF;

        IF edgeMapTable IS NOT NULL THEN
          sql := 'SELECT arc_id-1 FROM ' || edgeMapTable::text || ' WHERE edge_id = $1';
          EXECUTE sql INTO arcid USING rec.edge_id;
          IF arcid IS NULL THEN
            EXECUTE 'INSERT INTO ' || edgeMapTable::text
              || '(edge_id) VALUES ($1) RETURNING arc_id-1'
            INTO arcid USING rec.edge_id;
          END IF;
        ELSE
          arcid := rec.edge_id-1;
        END IF;

        -- Swap sign, use two's complement for negative edges
        IF rec.signed_edge_id >= 0 THEN
          arcid := - ( arcid + 1 );
        END IF;


        arcs := arcid || arcs;

      END LOOP; -- } over recursive query

      DELETE from _postgis_topology_astopojson_tmp_edges
      WHERE edge_id = ANY (visited_edges);
      visited_edges := ARRAY[]::int[];


      IF NOT edges_found THEN -- {

        IF looking_for_holes THEN
          looking_for_holes := false;
          comptxt := comptxt || ( '[' || array_to_string(ringtxt, ',') || ']' );
          ringtxt := NULL;
          faces := all_faces;
          shell_faces := ARRAY[]::int[];
        ELSE
          EXIT; -- end of loop
        END IF;

      ELSE -- } edges found {

        faces := shell_faces;
        IF arcs IS NOT NULL THEN
          ringtxt := ringtxt || ( '[' || array_to_string(arcs,',') || ']' );
        END IF;
        looking_for_holes := true;

      END IF; -- }

    END LOOP; -- }

    DROP TABLE _postgis_topology_astopojson_tmp_edges;

    json := json || array_to_string(comptxt, ',') || ']}';

    EXECUTE 'SET search_path TO ' || old_search_path;

  ELSIF tg.type = 4 THEN -- collection
    RAISE EXCEPTION 'Collection TopoGeometries are not supported by AsTopoJSON';

  END IF;

  RETURN json;

END
$$ LANGUAGE 'plpgsql' VOLATILE; -- writes into visited table
-- } AsTopoJSON(TopoGeometry, visited_table)


--=}  POSTGIS-SPECIFIC block

--  SQL/MM block
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2010-2015 Sandro Santilli <strk@kbt.io>
-- Copyright (C) 2005 Refractions Research Inc.
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- Author: Sandro Santilli <strk@kbt.io>
--
--



--={ ----------------------------------------------------------------
--  SQL/MM block
--
--  This part contains function in the SQL/MM specification
--
---------------------------------------------------------------------

--
-- Type returned by ST_GetFaceEdges
--
-- Availability: 1.1.0
--
CREATE TYPE topology.GetFaceEdges_ReturnType AS (
  sequence integer,
  edge integer
);

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.5
--
--  ST_GetFaceEdges(atopology, aface)
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.ST_GetFaceEdges(toponame varchar, face_id integer)
  RETURNS SETOF topology.GetFaceEdges_ReturnType AS
	'$libdir/postgis_topology-3', 'ST_GetFaceEdges'
  LANGUAGE 'c' STABLE;
--} ST_GetFaceEdges

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.10
--
--  ST_NewEdgeHeal(atopology, anedge, anotheredge)
--
-- Not in the specs:
-- * Refuses to heal two edges if any of the two is closed
-- * Raise an exception when trying to heal an edge with itself
-- * Raise an exception if any TopoGeometry is defined by only one
--   of the two edges
-- * Update references in the Relation table.
--
CREATE OR REPLACE FUNCTION topology.ST_NewEdgeHeal(toponame varchar, e1id integer, e2id integer)
  RETURNS int AS
  '$libdir/postgis_topology-3','ST_NewEdgeHeal'
  LANGUAGE 'c' VOLATILE;
--} ST_NewEdgeHeal

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.11
--
--  ST_ModEdgeHeal(atopology, anedge, anotheredge)
--
-- Not in the specs:
-- * Returns the id of the node being removed
-- * Refuses to heal two edges if any of the two is closed
-- * Raise an exception when trying to heal an edge with itself
-- * Raise an exception if any TopoGeometry is defined by only one
--   of the two edges
-- * Update references in the Relation table.
-- {
CREATE OR REPLACE FUNCTION topology.ST_ModEdgeHeal(toponame varchar, e1id integer, e2id integer)
  RETURNS int AS
  '$libdir/postgis_topology-3','ST_ModEdgeHeal'
  LANGUAGE 'c' VOLATILE;
--} ST_ModEdgeHeal

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.14
--
--  ST_RemEdgeNewFace(atopology, anedge)
--
-- Not in the specs:
-- * Raise an exception if any TopoGeometry is defined by only one
--   of the two faces that will dissolve.
-- * Raise an exception if any TopoGeometry is defined by
--   the edge being removed.
-- * Properly set containg_face on nodes that remains isolated by the drop
-- * Update containg_face for isolated nodes in the dissolved faces
-- * Update references in the Relation table
--
-- }{
CREATE OR REPLACE FUNCTION topology.ST_RemEdgeNewFace(toponame varchar, e1id integer)
  RETURNS int AS
	'$libdir/postgis_topology-3','ST_RemEdgeNewFace'
  LANGUAGE 'c' VOLATILE;
--} ST_RemEdgeNewFace

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.15
--
--  ST_RemEdgeModFace(atopology, anedge)
--
-- Not in the specs:
-- * Raise an exception if any TopoGeometry is defined by only one
--   of the two faces that will dissolve.
-- * Raise an exception if any TopoGeometry is defined by
--   the edge being removed.
-- * Properly set containg_face on nodes that remains isolated by the drop
-- * Update containg_face for isolated nodes in the dissolved faces
-- * Update references in the Relation table
-- * Return id of the face taking up the removed edge space
--
-- }{
CREATE OR REPLACE FUNCTION topology.ST_RemEdgeModFace(toponame varchar, e1id integer)
  RETURNS int AS
	'$libdir/postgis_topology-3','ST_RemEdgeModFace'
  LANGUAGE 'c' VOLATILE;
--} ST_RemEdgeModFace

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.16
--
--  ST_GetFaceGeometry(atopology, aface)
--
CREATE OR REPLACE FUNCTION topology.ST_GetFaceGeometry(toponame varchar, aface integer)
  RETURNS GEOMETRY AS
	'$libdir/postgis_topology-3', 'ST_GetFaceGeometry'
  LANGUAGE 'c' STABLE;
--} ST_GetFaceGeometry

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.1
--
--  ST_AddIsoNode(atopology, aface, apoint)
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.ST_AddIsoNode(atopology varchar, aface integer, apoint geometry)
  RETURNS INTEGER AS
	'$libdir/postgis_topology-3','ST_AddIsoNode'
  LANGUAGE 'c' VOLATILE;
--} ST_AddIsoNode

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.2
--
--  ST_MoveIsoNode(atopology, anode, apoint)
--
CREATE OR REPLACE FUNCTION topology.ST_MoveIsoNode(atopology character varying, anode integer, apoint geometry)
  RETURNS text AS
  '$libdir/postgis_topology-3','ST_MoveIsoNode'
  LANGUAGE 'c' VOLATILE;
--} ST_MoveIsoNode

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.3
--
--  ST_RemoveIsoNode(atopology, anode)
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.ST_RemoveIsoNode(atopology varchar, anode integer)
  RETURNS TEXT AS
	'$libdir/postgis_topology-3','ST_RemoveIsoNode'
  LANGUAGE 'c' VOLATILE;
--} ST_RemoveIsoNode

--{
-- According to http://trac.osgeo.org/postgis/ticket/798
-- ST_RemoveIsoNode was renamed to ST_RemIsoNode in the final ISO
-- document
--
CREATE OR REPLACE FUNCTION topology.ST_RemIsoNode(varchar, integer)
  RETURNS TEXT AS
	'$libdir/postgis_topology-3','ST_RemoveIsoNode'
  LANGUAGE 'c' VOLATILE;
--} ST_RemIsoNode

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.7
--
--  ST_RemoveIsoEdge(atopology, anedge)
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.ST_RemoveIsoEdge(atopology varchar, anedge integer)
  RETURNS TEXT AS
	'$libdir/postgis_topology-3','ST_RemIsoEdge'
	LANGUAGE 'c' VOLATILE;
--} ST_RemoveIsoEdge

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.8
--
--  ST_NewEdgesSplit(atopology, anedge, apoint)
--
-- Not in the specs:
-- * Update references in the Relation table.
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.ST_NewEdgesSplit(atopology varchar, anedge integer, apoint geometry)
  RETURNS INTEGER AS
  '$libdir/postgis_topology-3','ST_NewEdgesSplit'
  LANGUAGE 'c' VOLATILE;
--} ST_NewEdgesSplit

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.9
--
--  ST_ModEdgeSplit(atopology, anedge, apoint)
--
-- Not in the specs:
-- * Update references in the Relation table.
--
-- Availability: 2.0.0
-- Replaces ST_ModEdgesSplit(varchar, integer, geometry) deprecated in 2.0.0
--
CREATE OR REPLACE FUNCTION topology.ST_ModEdgeSplit(atopology varchar, anedge integer, apoint geometry)
  RETURNS INTEGER AS
	'$libdir/postgis_topology-3','ST_ModEdgeSplit'
  LANGUAGE 'c' VOLATILE;
--} ST_ModEdgesSplit

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.4
--
--  ST_AddIsoEdge(atopology, anode, anothernode, acurve)
--
-- Not in the specs:
-- * Reset containing_face for starting and ending point,
--   as they stop being isolated nodes
-- * Refuse to add a closed edge, as it would not be isolated
--   (ie: would create a ring)
--
-- Availability: 1.1.0
--
-- }{
--
CREATE OR REPLACE FUNCTION topology.ST_AddIsoEdge(atopology varchar, anode integer, anothernode integer, acurve geometry)
  RETURNS INTEGER AS
	'$libdir/postgis_topology-3','ST_AddIsoEdge'
  LANGUAGE 'c' VOLATILE;
--} ST_AddIsoEdge

-- Internal function used by ST_ChangeEdgeGeom to compare
-- adjacent edges of an edge endpoint
--
-- @param anode the node to use edge end star of
-- @param anedge the directed edge to get adjacents from
--        if positive `anode' is assumed to be its start node
--        if negative `anode' is assumed to be its end node
--
-- @todo DROP, NOT NEEDED ANYMORE (might need to go in some _drop, I guess)
-- {
CREATE OR REPLACE FUNCTION topology._ST_AdjacentEdges(atopology varchar, anode integer, anedge integer)
RETURNS integer[] AS
$$
DECLARE
  ret integer[];
BEGIN
  WITH edgestar AS (
    SELECT *, count(*) over () AS cnt
    FROM topology.GetNodeEdges(atopology, anode)
  )
  SELECT ARRAY[ (
      SELECT p.edge AS prev FROM edgestar p
      WHERE p.sequence = CASE WHEN m.sequence-1 < 1 THEN cnt
                         ELSE m.sequence-1 END
    ), (
      SELECT p.edge AS prev FROM edgestar p WHERE p.sequence = ((m.sequence)%cnt)+1
    ) ]
  FROM edgestar m
  WHERE edge = anedge
  INTO ret;

  RETURN ret;
END
$$
LANGUAGE 'plpgsql' STABLE;
--}

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.6
--
--  ST_ChangeEdgeGeom(atopology, anedge, acurve)
--
-- Not in the specs:
-- * Raise an exception if given a non-existent edge
-- * Raise an exception if movement is not topologically isomorphic
--
-- }{
CREATE OR REPLACE FUNCTION topology.ST_ChangeEdgeGeom(atopology varchar, anedge integer, acurve geometry)
  RETURNS TEXT AS
	'$libdir/postgis_topology-3','ST_ChangeEdgeGeom'
  LANGUAGE 'c' VOLATILE;
--} ST_ChangeEdgeGeom

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.12
--
--  ST_AddEdgeNewFaces(atopology, anode, anothernode, acurve)
--
-- Not in the specs:
-- * Reset containing_face for starting and ending point,
--   as they stop being isolated nodes
-- * Update references in the Relation table.
--
CREATE OR REPLACE FUNCTION topology.ST_AddEdgeNewFaces(atopology varchar, anode integer, anothernode integer, acurve geometry)
  RETURNS INTEGER AS
	'$libdir/postgis_topology-3','ST_AddEdgeNewFaces'
  LANGUAGE 'c' VOLATILE;
--} ST_AddEdgeNewFaces

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.13
--
--  ST_AddEdgeModFace(atopology, anode, anothernode, acurve)
--
-- Not in the specs:
-- * Reset containing_face for starting and ending point,
--   as they stop being isolated nodes
-- * Update references in the Relation table.
--
CREATE OR REPLACE FUNCTION topology.ST_AddEdgeModFace(atopology varchar, anode integer, anothernode integer, acurve geometry)
  RETURNS INTEGER AS
	'$libdir/postgis_topology-3','ST_AddEdgeModFace'
  LANGUAGE 'c' VOLATILE;
--} ST_AddEdgeModFace

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.17
--
--  ST_InitTopoGeo(atopology)
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.ST_InitTopoGeo(atopology varchar)
RETURNS text
AS
$$
DECLARE
  rec RECORD;
  topology_id numeric;
BEGIN
  IF atopology IS NULL THEN
    RAISE EXCEPTION 'SQL/MM Spatial exception - null argument';
  END IF;

  FOR rec IN SELECT * FROM pg_namespace WHERE text(nspname) = atopology
  LOOP
    RAISE EXCEPTION 'SQL/MM Spatial exception - schema already exists';
  END LOOP;

  FOR rec IN EXECUTE 'SELECT topology.CreateTopology('
    ||quote_literal(atopology)|| ') as id'
  LOOP
    topology_id := rec.id;
  END LOOP;

  RETURN 'Topology-Geometry ' || quote_literal(atopology)
    || ' (id:' || topology_id || ') created.';
END
$$
LANGUAGE 'plpgsql' VOLATILE;
--} ST_InitTopoGeo

--{
-- Topo-Geo and Topo-Net 3: Routine Details
-- X.3.18
--
--  ST_CreateTopoGeo(atopology, acollection)
--}{
CREATE OR REPLACE FUNCTION topology.ST_CreateTopoGeo(atopology varchar, acollection geometry)
RETURNS text
AS
$$
DECLARE
  typ char(4);
  rec RECORD;
  ret int;
  nodededges GEOMETRY;
  points GEOMETRY;
  snode_id int;
  enode_id int;
  tolerance FLOAT8;
  topoinfo RECORD;
BEGIN

  IF atopology IS NULL OR acollection IS NULL THEN
    RAISE EXCEPTION 'SQL/MM Spatial exception - null argument';
  END IF;

  -- Get topology information
  BEGIN
    SELECT * FROM topology.topology
      INTO STRICT topoinfo WHERE name = atopology;
  EXCEPTION
    WHEN NO_DATA_FOUND THEN
      RAISE EXCEPTION 'SQL/MM Spatial exception - invalid topology name';
  END;

  -- Check SRID compatibility
  IF ST_SRID(acollection) != topoinfo.SRID THEN
    RAISE EXCEPTION 'Geometry SRID (%) does not match topology SRID (%)',
      ST_SRID(acollection), topoinfo.SRID;
  END IF;

  -- Verify pre-conditions (valid, empty topology schema exists)
  BEGIN -- {

    -- Verify the topology views in the topology schema to be empty
    FOR rec in EXECUTE
      'SELECT count(*) FROM '
      || quote_ident(atopology) || '.edge_data '
      || ' UNION ' ||
      'SELECT count(*) FROM '
      || quote_ident(atopology) || '.node '
    LOOP
      IF rec.count > 0 THEN
    RAISE EXCEPTION 'SQL/MM Spatial exception - non-empty view';
      END IF;
    END LOOP;

    -- face check is separated as it will contain a single (world)
    -- face record
    FOR rec in EXECUTE
      'SELECT count(*) FROM '
      || quote_ident(atopology) || '.face '
    LOOP
      IF rec.count != 1 THEN
    RAISE EXCEPTION 'SQL/MM Spatial exception - non-empty face view';
      END IF;
    END LOOP;

  EXCEPTION
    WHEN INVALID_SCHEMA_NAME THEN
      RAISE EXCEPTION 'SQL/MM Spatial exception - invalid topology name';
    WHEN UNDEFINED_TABLE THEN
      RAISE EXCEPTION 'SQL/MM Spatial exception - non-existent view';

  END; -- }


  --
  -- Node input linework with itself
  --
  WITH components AS ( SELECT geom FROM ST_Dump(acollection) )
  SELECT ST_UnaryUnion(ST_Collect(geom)) FROM (
    SELECT geom FROM components
      WHERE ST_Dimension(geom) = 1
    UNION ALL
    SELECT ST_Boundary(geom) FROM components
      WHERE ST_Dimension(geom) = 2
  ) as linework INTO STRICT nodededges;


  --
  -- Linemerge the resulting edges, to reduce the working set
  -- NOTE: this is more of a workaround for GEOS splitting overlapping
  --       lines to each of the segments.
  --
  SELECT ST_LineMerge(nodededges) INTO STRICT nodededges;


  --
  -- Collect input points and input lines endpoints
  --
  WITH components AS ( SELECT geom FROM ST_Dump(acollection) )
  SELECT ST_Union(geom) FROM (
    SELECT geom FROM components
      WHERE ST_Dimension(geom) = 0
    UNION ALL
    SELECT ST_Boundary(geom) FROM components
      WHERE ST_Dimension(geom) = 1
  ) as nodes INTO STRICT points;


  --
  -- Further split edges by points, if needed
  --
  IF points IS NOT NULL THEN
    nodededges := ST_Split(nodededges, points);
  END IF; -- points is not null

  --
  -- Collect all nodes (from points and noded linework endpoints)
  --

  WITH edges AS ( SELECT geom FROM ST_Dump(nodededges) )
  SELECT ST_Union( -- TODO: ST_UnaryUnion ?
          COALESCE(ST_UnaryUnion(ST_Collect(geom)),
            ST_SetSRID('POINT EMPTY'::geometry, topoinfo.SRID)),
          COALESCE(points,
            ST_SetSRID('POINT EMPTY'::geometry, topoinfo.SRID))
         )
  FROM (
    SELECT ST_StartPoint(geom) as geom FROM edges
      UNION ALL
    SELECT ST_EndPoint(geom) FROM edges
  ) as endpoints INTO points;


  --
  -- Add all nodes as isolated so that
  -- later calls to AddEdgeModFace will tweak their being
  -- isolated or not...
  --
  FOR rec IN SELECT geom FROM ST_Dump(points)
  LOOP
    PERFORM topology.ST_AddIsoNode(atopology, 0, rec.geom);
  END LOOP;

  FOR rec IN SELECT geom FROM ST_Dump(nodededges)
  LOOP
    SELECT topology.GetNodeByPoint(atopology, st_startpoint(rec.geom), 0)
      INTO STRICT snode_id;
    SELECT topology.GetNodeByPoint(atopology, st_endpoint(rec.geom), 0)
      INTO STRICT enode_id;
    PERFORM topology.ST_AddEdgeModFace(atopology, snode_id, enode_id, rec.geom);
  END LOOP;

  RETURN 'Topology ' || atopology || ' populated';

END
$$
LANGUAGE 'plpgsql' VOLATILE;
--} ST_CreateTopoGeo

--=}  SQL/MM block


-- The following files needs getfaceedges_returntype, defined in sqlmm.sql
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011 2012 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- Developed by Sandro Santilli <strk@kbt.io>
-- 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)"
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--
-- Return a list of edges (sequence, id) resulting by starting from the
-- given edge and following the leftmost turn at each encountered node.
--
-- Edge ids are signed, they are negative if traversed backward.
-- Sequence numbers start with 1.
--
-- Use a negative starting_edge to follow its rigth face rather than
-- left face (to start traversing it in reverse).
--
-- Optionally pass a limit on the number of edges to traverse. This is a
-- safety measure against not-properly linked topologies, where you may
-- end up looping forever (single edge loops edge are detected but longer
-- ones are not). Default is no limit (good luck!)
--
-- GetRingEdges(atopology, anedge, [maxedges])
--
CREATE OR REPLACE FUNCTION topology.GetRingEdges(atopology varchar, anedge int, maxedges int DEFAULT null)
	RETURNS SETOF topology.GetFaceEdges_ReturnType AS
	'$libdir/postgis_topology-3', 'GetRingEdges'
LANGUAGE 'c' STABLE;
--} GetRingEdges
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2012 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--
-- Return a list of edges (sequence, id) incident to the given node.
--
-- Edge ids are signed, they are negative if the node is their endpoint.
-- Sequence numbers start with 1 ordering edges by azimuth (clockwise).
--
-- GetNodeEdges(atopology, anode)
--
CREATE OR REPLACE FUNCTION topology.GetNodeEdges(atopology varchar, anode int)
	RETURNS SETOF topology.GetFaceEdges_ReturnType
AS
$BODY$
DECLARE
  curedge int;
  nextedge int;
  rec RECORD;
  retrec topology.GetFaceEdges_ReturnType;
  n int;
  sql text;
BEGIN

  n := 0;
  sql := format(
    $$
      WITH incident_edges AS (
        SELECT
          edge_id,
          start_node,
          end_node,
          ST_RemoveRepeatedPoints(geom) as geom
        FROM %1$I.edge_data
        WHERE start_node = $1
        or end_node = $1
      )
      SELECT
        edge_id,
        ST_Azimuth(ST_StartPoint(geom), ST_PointN(geom, 2)) as az
      FROM incident_edges
      WHERE start_node = $1
        UNION ALL
      SELECT
        -edge_id,
        ST_Azimuth(ST_EndPoint(geom), ST_PointN(geom, ST_NumPoints(geom)-1))
      FROM incident_edges
      WHERE end_node = $1
      ORDER BY az
    $$,
    atopology
  );

  FOR rec IN EXECUTE sql USING anode
  LOOP -- incident edges {

    n := n + 1;
    retrec.sequence := n;
    retrec.edge := rec.edge_id;
    RETURN NEXT retrec;
  END LOOP; -- incident edges }

END
$BODY$
LANGUAGE 'plpgsql' STABLE;
--} GetNodeEdges

--general management --
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://www.postgis.net
--
-- Copyright (C) 2011 Regina Obe <lr@pcorp.us>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--  AddToSearchPath(schema_name)
--
-- Adds the specified schema to the database search path
-- if it is not already in the database search path
-- This is a helper function for upgrade/install
-- We may want to move this function as a generic helper
CREATE OR REPLACE FUNCTION topology.AddToSearchPath(a_schema_name varchar)
RETURNS text
AS
$$
DECLARE
	var_result text;
	var_cur_search_path text;
BEGIN
	SELECT reset_val INTO var_cur_search_path FROM pg_catalog.pg_settings WHERE name = 'search_path';
	IF var_cur_search_path LIKE '%' || pg_catalog.quote_ident(a_schema_name) || '%' THEN
		var_result := a_schema_name || ' already in database search_path';
	ELSE
		var_cur_search_path := var_cur_search_path || ', '
                        || pg_catalog.quote_ident(a_schema_name);
		EXECUTE 'ALTER DATABASE ' || pg_catalog.quote_ident(pg_catalog.current_database())
                              || ' SET search_path = ' || var_cur_search_path;
		var_result := a_schema_name || ' has been added to end of database search_path ';
	END IF;

	EXECUTE 'SET search_path = ' || var_cur_search_path;

  RETURN var_result;
END
$$
LANGUAGE 'plpgsql' VOLATILE STRICT;

--} AddToSearchPath
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2022 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--  AddTopoGeometryColumn(toponame, schema, table, colum, type, [child])
--
--  Add a TopoGeometry column to a table, making it a topology layer.
--  Returns created layer id.
--
-- Availability: 1.1.0
--
CREATE OR REPLACE FUNCTION topology.AddTopoGeometryColumn(toponame varchar, schema varchar, tbl varchar, col varchar, ltype varchar, child integer)
  RETURNS integer
AS
$$
DECLARE
  intltype integer;
  newlevel integer;
  topoid integer;
  rec RECORD;
  newlayer_id integer;
  query text;
BEGIN

  -- Get topology id
  SELECT id INTO topoid
    FROM topology.topology WHERE name = toponame;

  IF NOT FOUND THEN
    RAISE EXCEPTION 'Topology % does not exist', quote_literal(toponame);
  END IF;

  IF ltype ILIKE '%POINT%' OR ltype ILIKE 'PUNTAL' THEN
    intltype = 1;
  ELSIF ltype ILIKE '%LINE%' OR ltype ILIKE 'LINEAL' THEN
    intltype = 2;
  ELSIF ltype ILIKE '%POLYGON%' OR ltype ILIKE 'AREAL' THEN
    intltype = 3;
  ELSIF ltype ILIKE '%COLLECTION%' OR ltype ILIKE 'GEOMETRY' THEN
    intltype = 4;
  ELSE
    RAISE EXCEPTION 'Layer type must be one of POINT,LINE,POLYGON,COLLECTION';
  END IF;

  --
  -- Add new TopoGeometry column in schema.table
  --
  EXECUTE 'ALTER TABLE ' || quote_ident(schema)
    || '.' || quote_ident(tbl)
    || ' ADD COLUMN ' || quote_ident(col)
    || ' topology.TopoGeometry;';

  --
  -- See if child id exists and extract its level
  --
  IF child IS NOT NULL THEN
    SELECT level + 1 FROM topology.layer
      WHERE layer_id = child
        AND topology_id = topoid
      INTO newlevel;
    IF newlevel IS NULL THEN
      RAISE EXCEPTION 'Child layer % does not exist in topology "%"', child, toponame;
    END IF;
  END IF;

  --
  -- Get new layer id from sequence
  --
  EXECUTE 'SELECT nextval(' ||
    quote_literal(
      quote_ident(toponame) || '.layer_id_seq'
    ) || ')' INTO STRICT newlayer_id;

  EXECUTE 'INSERT INTO '
       'topology.layer(topology_id, '
       'layer_id, level, child_id, schema_name, '
       'table_name, feature_column, feature_type) '
       'VALUES ('
    || topoid || ','
    || newlayer_id || ',' || COALESCE(newlevel, 0) || ','
    || COALESCE(child::text, 'NULL') || ','
    || quote_literal(schema) || ','
    || quote_literal(tbl) || ','
    || quote_literal(col) || ','
    || intltype || ');';

  --
  -- Create a sequence for TopoGeometries in this new layer
  --
  EXECUTE 'CREATE SEQUENCE ' || quote_ident(toponame)
    || '.topogeo_s_' || newlayer_id;

  --
  -- Add constraints on TopoGeom column
  --
  EXECUTE 'ALTER TABLE ' || quote_ident(schema)
    || '.' || quote_ident(tbl)
    || ' ADD CONSTRAINT "check_topogeom_' || col || '" CHECK ('
       'topology_id(' || quote_ident(col) || ') = ' || topoid
    || ' AND '
       'layer_id(' || quote_ident(col) || ') = ' || newlayer_id
    || ' AND '
       'type(' || quote_ident(col) || ') = ' || intltype
    || ');';

  --
  -- Add dependency of the feature column on the topology schema
  --
  query = 'INSERT INTO pg_catalog.pg_depend SELECT '
       'fcat.oid, fobj.oid, fsub.attnum, tcat.oid, '
       'tobj.oid, 0, ''n'' '
       'FROM pg_class fcat, pg_namespace fnsp, '
       ' pg_class fobj, pg_attribute fsub, '
       ' pg_class tcat, pg_namespace tobj '
       ' WHERE fcat.relname = ''pg_class'' '
       ' AND fnsp.nspname = ' || quote_literal(schema)
    || ' AND fobj.relnamespace = fnsp.oid '
       ' AND fobj.relname = ' || quote_literal(tbl)
    || ' AND fsub.attrelid = fobj.oid '
       ' AND fsub.attname = ' || quote_literal(col)
    || ' AND tcat.relname = ''pg_namespace'' '
       ' AND tobj.nspname = ' || quote_literal(toponame);

--
-- The only reason to add this dependency is to avoid
-- simple drop of a feature column. Still, drop cascade
-- will remove both the feature column and the sequence
-- corrupting the topology anyway ...
--

  RETURN newlayer_id;
END;
$$
LANGUAGE 'plpgsql' VOLATILE;
--}{ AddTopoGeometryColumn

CREATE OR REPLACE FUNCTION topology.AddTopoGeometryColumn(varchar, varchar, varchar, varchar, varchar)
  RETURNS integer
AS
$$
  SELECT topology.AddTopoGeometryColumn($1, $2, $3, $4, $5, NULL);
$$
LANGUAGE 'sql' VOLATILE;

--
--} AddTopoGeometryColumn

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2010, 2011 Sandro Santilli <strk@kbt.io>
-- Copyright (C) 2005 Refractions Research Inc.
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- Author: Sandro Santilli <strk@kbt.io>
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--  CreateTopology(name, SRID, precision, hasZ)
--
-- Create a topology schema, add a topology info record
-- in the topology.topology relation, return it's numeric
-- id.
--
-- Availability: 2.0.0
--
CREATE OR REPLACE FUNCTION topology.CreateTopology(atopology varchar, srid integer, prec float8, hasZ boolean)
RETURNS integer
AS
$BODY$
DECLARE
  rec RECORD;
  topology_id integer;
  sql TEXT;
  zsuffix TEXT := '';
BEGIN

--  FOR rec IN SELECT * FROM pg_namespace WHERE text(nspname) = atopology
--  LOOP
--    RAISE EXCEPTION 'SQL/MM Spatial exception - schema already exists';
--  END LOOP;

  IF hasZ THEN
    zsuffix := 'z';
  END IF;

  IF srid < 0 THEN
    RAISE NOTICE 'SRID value % converted to the officially unknown SRID value 0', srid;
    srid = 0;
  END IF;

  ------ Fetch next id for the new topology
  FOR rec IN SELECT nextval('topology.topology_id_seq')
  LOOP
    topology_id = rec.nextval;
  END LOOP;

  sql := format(
    $$
      CREATE SCHEMA %1$I;

      -------------{ face CREATION
      CREATE TABLE %1$I.face(
        face_id SERIAL,
        mbr GEOMETRY(polygon, %2$L), -- 2d only mbr is good enough
        CONSTRAINT face_primary_key
        PRIMARY KEY(face_id)
      );

      -- Face standard view description
      COMMENT ON TABLE %1$I.face IS
      'Contains face topology primitives';

      -------------} END OF face CREATION

      --------------{ node CREATION

      CREATE TABLE %1$I.node(
        node_id SERIAL,
        containing_face INTEGER,
        geom GEOMETRY(point%3$s, %2$L),
        CONSTRAINT node_primary_key
          PRIMARY KEY(node_id),
        CONSTRAINT face_exists
          FOREIGN KEY(containing_face)
          REFERENCES %1$I.face(face_id)
      );

      -- Node standard view description
      COMMENT ON TABLE %1$I.node IS
      'Contains node topology primitives';

      ------- Add index on containing_face, to speed up
      ------- topology editing (adding/removing faces)
      ------- See http://trac.osgeo.org/postgis/ticket/2861
      CREATE INDEX node_containing_face_idx
        ON %1$I.node (containing_face);

      --------------} END OF node CREATION

      --------------{ edge CREATION

      -- edge_data table
      CREATE TABLE %1$I.edge_data (
        edge_id SERIAL NOT NULL PRIMARY KEY,
        start_node INTEGER NOT NULL,
        end_node INTEGER NOT NULL,
        next_left_edge INTEGER NOT NULL,
        abs_next_left_edge INTEGER NOT NULL,
        next_right_edge INTEGER NOT NULL,
        abs_next_right_edge INTEGER NOT NULL,
        left_face INTEGER NOT NULL,
        right_face INTEGER NOT NULL,
        geom GEOMETRY(linestring%3$s, %2$L) NOT NULL,

        CONSTRAINT start_node_exists
          FOREIGN KEY(start_node)
          REFERENCES %1$I.node(node_id),

        CONSTRAINT end_node_exists
          FOREIGN KEY(end_node)
          REFERENCES %1$I.node(node_id),

        CONSTRAINT left_face_exists
          FOREIGN KEY(left_face)
          REFERENCES %1$I.face(face_id),

        CONSTRAINT right_face_exists
          FOREIGN KEY(right_face)
          REFERENCES %1$I.face(face_id),

        CONSTRAINT next_left_edge_exists
          FOREIGN KEY(abs_next_left_edge)
          REFERENCES %1$I.edge_data(edge_id)
          DEFERRABLE INITIALLY DEFERRED,

        CONSTRAINT next_right_edge_exists
          FOREIGN KEY(abs_next_right_edge)
          REFERENCES %1$I.edge_data(edge_id)
          DEFERRABLE INITIALLY DEFERRED
      );

      -- edge standard view (select rule)
      CREATE VIEW %1$I.edge AS
      SELECT
        edge_id, start_node, end_node, next_left_edge,
        next_right_edge, left_face, right_face, geom
      FROM %1$I.edge_data;

      -- Edge standard view description
      COMMENT ON VIEW %1$I.edge IS
      'Contains edge topology primitives';
      COMMENT ON COLUMN %1$I.edge.edge_id IS
      'Unique identifier of the edge';
      COMMENT ON COLUMN %1$I.edge.start_node IS
      'Unique identifier of the node at the start of the edge';
      COMMENT ON COLUMN %1$I.edge.end_node IS
      'Unique identifier of the node at the end of the edge';
      COMMENT ON COLUMN %1$I.edge.next_left_edge IS
      'Unique identifier of the next edge of the face on the left (when looking in the direction from START_NODE to END_NODE), moving counterclockwise around the face boundary';
      COMMENT ON COLUMN %1$I.edge.next_right_edge IS
      'Unique identifier of the next edge of the face on the right (when looking in the direction from START_NODE to END_NODE), moving counterclockwise around the face boundary';
      COMMENT ON COLUMN %1$I.edge.left_face IS
      'Unique identifier of the face on the left side of the edge when looking in the direction from START_NODE to END_NODE';
      COMMENT ON COLUMN %1$I.edge.right_face IS
      'Unique identifier of the face on the right side of the edge when looking in the direction from START_NODE to END_NODE';
      COMMENT ON COLUMN %1$I.edge.geom IS
      'The geometry of the edge';

      -- edge standard view (insert rule)
      CREATE RULE edge_insert_rule AS
      ON INSERT TO %1$I.edge
      DO INSTEAD INSERT into %1$I.edge_data
      VALUES (
        NEW.edge_id, NEW.start_node, NEW.end_node,
        NEW.next_left_edge, abs(NEW.next_left_edge),
        NEW.next_right_edge, abs(NEW.next_right_edge),
        NEW.left_face, NEW.right_face, NEW.geom
      );

      ------- Add support indices supporting edge linking foreign keys
      ------- See https://trac.osgeo.org/postgis/ticket/2083#comment:20
      CREATE INDEX ON %1$I.edge_data(abs_next_left_edge);
      CREATE INDEX ON %1$I.edge_data(abs_next_right_edge);

      --------------} END OF edge CREATION

      --------------{ layer sequence
      CREATE SEQUENCE %1$I.layer_id_seq;
      --------------} layer sequence

      --------------{ relation CREATION
      CREATE TABLE %1$I.relation (
        topogeo_id integer NOT NULL,
        layer_id integer NOT NULL,
        element_id integer NOT NULL,
        element_type integer NOT NULL,
        UNIQUE(layer_id,topogeo_id,element_id,element_type)
      );
      ------- Add index on element_type, element_id to speed up
      ------- queries looking for TopoGeometries using specific
      ------- primitive elements of the topology
      ------- See http://trac.osgeo.org/postgis/ticket/2083
      CREATE INDEX relation_element_id_idx
        ON %1$I.relation (element_id);

      CREATE TRIGGER relation_integrity_checks
      BEFORE UPDATE OR INSERT ON %1$I.relation
      FOR EACH ROW EXECUTE PROCEDURE
      topology.RelationTrigger(%4$L, %1$L);
      --------------} END OF relation CREATION

      ------- Default (world) face
      INSERT INTO %1$I.face(face_id) VALUES(0);

      ------- GiST index on face
      CREATE INDEX face_gist ON %1$I.face
      USING gist (mbr);

      ------- GiST index on node
      CREATE INDEX node_gist ON %1$I.node
      USING gist (geom);

      ------- GiST index on edge
      CREATE INDEX edge_gist ON %1$I.edge_data
      USING gist (geom);

      ------- Indexes on left_face and right_face of edge_data
      ------- NOTE: these indexes speed up GetFaceGeometry (and thus
      -------       TopoGeometry::Geometry) by a factor of 10 !
      -------       See http://trac.osgeo.org/postgis/ticket/806
      CREATE INDEX edge_left_face_idx
        ON %1$I.edge_data (left_face);
      CREATE INDEX edge_right_face_idx
        ON %1$I.edge_data (right_face);

      ------- Indexes on start_node and end_node of edge_data
      ------- NOTE: this indexes speed up node deletion
      -------       by a factor of 1000 !
      -------       See http://trac.osgeo.org/postgis/ticket/2082
      CREATE INDEX edge_start_node_idx
        ON %1$I.edge_data (start_node);
      CREATE INDEX edge_end_node_idx
        ON %1$I.edge_data (end_node);

      -- TODO: consider also adding an index on node.containing_face

    $$,
    atopology,   -- %1
    srid,        -- %2
    zsuffix,     -- %3
    topology_id  -- %4
  );
  EXECUTE sql;

  ------- Add record to the "topology" metadata table
  INSERT INTO topology.topology (id, name, srid, precision, hasZ)
  VALUES (topology_id, atopology, srid, prec, hasZ);

  RETURN topology_id;
END
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;

--} CreateTopology

--{ CreateTopology wrappers for unspecified srid or precision or hasZ

--  CreateTopology(name, SRID, precision) -- hasZ = false
-- Availability: 1.1.0
CREATE OR REPLACE FUNCTION topology.CreateTopology(toponame varchar, srid integer, prec float8)
RETURNS integer AS
' SELECT topology.CreateTopology($1, $2, $3, false);'
LANGUAGE 'sql' VOLATILE STRICT;

--  CreateTopology(name, SRID) -- precision = 0
-- Availability: 1.1.0
CREATE OR REPLACE FUNCTION topology.CreateTopology(varchar, integer)
RETURNS integer AS
' SELECT topology.CreateTopology($1, $2, 0); '
LANGUAGE 'sql' VOLATILE STRICT;

--  CreateTopology(name) -- srid = unknown, precision = 0
-- Availability: 1.1.0
CREATE OR REPLACE FUNCTION topology.CreateTopology(varchar)
RETURNS integer AS
$$ SELECT topology.CreateTopology($1, ST_SRID('POINT EMPTY'::geometry), 0); $$
LANGUAGE 'sql' VOLATILE STRICT;

--} CreateTopology

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--  TopologySummary(name)
--
-- Print an overview about a topology
--
CREATE OR REPLACE FUNCTION topology.TopologySummary(atopology varchar)
RETURNS text
AS
$BODY$
DECLARE
  rec RECORD;
  rec2 RECORD;
  var_topology_id integer;
  n int4;
  missing int4;
  sql text;
  ret text;
  tgcount int4;
BEGIN

  ret := 'Topology ' || quote_ident(atopology) ;

  BEGIN
    SELECT * FROM topology.topology WHERE name = atopology INTO STRICT rec;
    -- TODO: catch <no_rows> to give a nice error message
    var_topology_id := rec.id;

    ret := ret || ' (id ' || rec.id || ', '
               || 'SRID ' || rec.srid || ', '
               || 'precision ' || rec.precision;
    IF rec.hasz THEN ret := ret || ', has Z'; END IF;
    ret := ret || E')\n';
  EXCEPTION
    WHEN NO_DATA_FOUND THEN
      ret := ret || E' (unregistered)\n';
  END;

  BEGIN
    EXECUTE 'SELECT count(*) FROM ' || quote_ident(atopology)
      || '.node ' INTO STRICT n;
    ret = ret || n || ' nodes, ';
  EXCEPTION
    WHEN UNDEFINED_TABLE OR INVALID_SCHEMA_NAME THEN
      IF NOT EXISTS (
          SELECT * FROM pg_catalog.pg_namespace WHERE nspname = atopology
         )
      THEN
        ret = ret || 'missing schema';
        RETURN ret;
      ELSE
        ret = ret || 'missing nodes, ';
      END IF;
  END;

  BEGIN
    EXECUTE 'SELECT count(*) FROM ' || quote_ident(atopology)
      || '.edge' INTO STRICT n;
    ret = ret || n || ' edges, ';
  EXCEPTION
    WHEN UNDEFINED_TABLE OR INVALID_SCHEMA_NAME THEN
      ret = ret || 'missing edges, ';
  END;

  BEGIN
    EXECUTE 'SELECT count(*) FROM ' || quote_ident(atopology)
      || '.face' INTO STRICT n;
    ret = ret || greatest(n-1,0) || ' faces, '; -- -1 is face=0
  EXCEPTION
    WHEN UNDEFINED_TABLE OR INVALID_SCHEMA_NAME THEN
      ret = ret || 'missing faces, ';
  END;

  BEGIN
    EXECUTE 'SELECT count(distinct layer_id) AS ln, '
      || 'count(distinct (layer_id,topogeo_id)) AS tn FROM '
      || quote_ident(atopology) || '.relation' INTO STRICT rec;
    tgcount := rec.tn;
    ret = ret || rec.tn || ' topogeoms in ' || rec.ln || E' layers\n';
  EXCEPTION
    WHEN UNDEFINED_TABLE THEN
      ret = ret || E'missing relations\n';
    WHEN UNDEFINED_COLUMN THEN
      ret = ret || E'corrupted relations\n';
  END;

  -- print information about registered layers
  FOR rec IN SELECT * FROM topology.layer l
    WHERE l.topology_id = var_topology_id
    ORDER by layer_id
  LOOP -- {
    ret = ret || 'Layer ' || rec.layer_id || ', type ';
    CASE
      WHEN rec.feature_type = 1 THEN
        ret = ret || 'Puntal';
      WHEN rec.feature_type = 2 THEN
        ret = ret || 'Lineal';
      WHEN rec.feature_type = 3 THEN
        ret = ret || 'Polygonal';
      WHEN rec.feature_type = 4 THEN
        ret = ret || 'Mixed';
      ELSE
        ret = ret || '???';
    END CASE;

    ret = ret || ' (' || rec.feature_type || '), ';

    BEGIN

      EXECUTE 'SELECT count(*) FROM ( SELECT DISTINCT topogeo_id FROM '
        || quote_ident(atopology)
        || '.relation r WHERE r.layer_id = ' || rec.layer_id
        || ' ) foo ' INTO STRICT n;

      ret = ret || n || ' topogeoms' || E'\n';

    EXCEPTION WHEN UNDEFINED_TABLE OR UNDEFINED_COLUMN THEN
      n := NULL;
      ret = ret || 'X topogeoms' || E'\n';
    END;

      IF rec.level > 0 THEN
        ret = ret || ' Hierarchy level ' || rec.level
                  || ', child layer ' || rec.child_id || E'\n';
      END IF;

      ret = ret || ' Deploy: ';
      IF rec.feature_column != '' THEN
        ret = ret || quote_ident(rec.schema_name) || '.'
                  || quote_ident(rec.table_name) || '.'
                  || quote_ident(rec.feature_column);

        IF n > 0 THEN
          sql := format(
            $$
SELECT count(*) FROM (
  SELECT topogeo_id
  FROM %1$I.relation r
  WHERE r.layer_id = %2$L
    EXCEPT
  SELECT DISTINCT id(%3$I)
  FROM %4$I.%5$I
  WHERE layer_id(%3$I) = %2$L
    AND topology_id(%3$I) = %6$L
) as foo
            $$,
            atopology,
            rec.layer_id,
            rec.feature_column,
            rec.schema_name,
            rec.table_name,
            var_topology_id
          );
          BEGIN
            RAISE DEBUG 'Executing %', sql;
            EXECUTE sql INTO STRICT missing;
            IF missing > 0 THEN
              ret = ret || ' (' || missing || ' missing topogeoms)';
            END IF;
          EXCEPTION
            WHEN UNDEFINED_TABLE THEN
              ret = ret || ' ( unexistent table )';
            WHEN UNDEFINED_COLUMN THEN
              ret = ret || ' ( unexistent column )';
          END;
        END IF;
        ret = ret || E'\n';

      ELSE
        ret = ret || E'NONE (detached)\n';
      END IF;

  END LOOP; -- }

  -- print information about unregistered layers containing topogeoms
  IF tgcount > 0 THEN -- {

    sql := 'SELECT layer_id FROM '
        || quote_ident(atopology) || '.relation EXCEPT SELECT layer_id'
        || ' FROM topology.layer WHERE topology_id = $1 ORDER BY layer_id';
    --RAISE DEBUG '%', sql;
    FOR rec IN  EXECUTE sql USING var_topology_id
    LOOP -- {
      ret = ret || 'Layer ' || rec.layer_id::text || ', UNREGISTERED, ';

      EXECUTE 'SELECT count(*) FROM ( SELECT DISTINCT topogeo_id FROM '
        || quote_ident(atopology)
        || '.relation r WHERE r.layer_id = ' || rec.layer_id
        || ' ) foo ' INTO STRICT n;

      ret = ret || n || ' topogeoms' || E'\n';

    END LOOP; -- }

  END IF; -- }

  RETURN ret;
END
$BODY$
LANGUAGE 'plpgsql' STABLE STRICT;

--} TopologySummary
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2011 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
-- CopyTopology(name_source, name_target)
--
-- Makes a copy of a topology (primitives + topogeometry collections) .
-- Returns the new topology id.
--
CREATE OR REPLACE FUNCTION topology.CopyTopology(atopology varchar, newtopo varchar)
RETURNS int
AS
$BODY$
DECLARE
  rec RECORD;
  rec2 RECORD;
  oldtopo_id integer;
  newtopo_id integer;
  n int4;
  ret text;
  sql text;
BEGIN

  SELECT * FROM topology.topology where name = atopology
  INTO strict rec;
  oldtopo_id = rec.id;
  -- TODO: more gracefully handle unexistent topology

  SELECT topology.CreateTopology(newtopo, rec.SRID, rec.precision, rec.hasZ)
  INTO strict newtopo_id;

  sql := format(
    $$
      -- Copy faces
      INSERT INTO %1$I.face
      SELECT * FROM %2$I.face
      WHERE face_id != 0;

      -- Update face sequence
      SELECT setval(
        '%1$I.face_face_id_seq',
        (SELECT last_value FROM %2$I.face_face_id_seq)
      );

      -- Copy nodes
      INSERT INTO %1$I.node
      SELECT * FROM %2$I.node;

      -- Update node sequence
      SELECT setval(
        '%1$I.node_node_id_seq',
        (SELECT last_value FROM %2$I.node_node_id_seq)
      );

      -- Copy edges
      INSERT INTO %1$I.edge_data
      SELECT * FROM %2$I.edge_data;

      -- Update edge sequence
      SELECT setval(
        '%1$I.edge_data_edge_id_seq',
        (SELECT last_value FROM %2$I.edge_data_edge_id_seq)
      );
    $$,
    newtopo,
    atopology
  );
  EXECUTE sql;

  -- Copy layers and their TopoGeometry sequences
  FOR rec IN SELECT * FROM topology.layer WHERE topology_id = oldtopo_id
  LOOP
    INSERT INTO topology.layer (topology_id, layer_id, feature_type,
      level, child_id, schema_name, table_name, feature_column)
      VALUES (newtopo_id, rec.layer_id, rec.feature_type,
              rec.level, rec.child_id, newtopo,
              'LAYER' ||  rec.layer_id, '');

    -- Create layer's TopoGeometry sequences
    EXECUTE format(
      $$
        CREATE SEQUENCE %1$I.topogeo_s_%2$s;
        SELECT setval(
          '%1$I.topogeo_s_%2$s',
          (SELECT last_value FROM %3$I.topogeo_s_%2$s)
        );
      $$,
      newtopo,
      rec.layer_id,
      atopology
    );
  END LOOP;

  -- Copy TopoGeometry definitions
  EXECUTE format(
    $$
      INSERT INTO %1$I.relation
      SELECT * FROM %2$I.relation
    $$,
    newtopo,
    atopology
  );

  RETURN newtopo_id;
END
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;

--} TopologySummary
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2021 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--  FindTopology(TopoGeometry)
--
-- Return a topology.topology record from a TopoGeometry
--
CREATE OR REPLACE FUNCTION topology.FindTopology(topology.TopoGeometry)
RETURNS topology.topology
AS $$
    SELECT * FROM topology.topology
    WHERE id = topology_id($1);
$$ LANGUAGE 'sql';
--}

--{
--  FindTopology(layerRegclass, layerColumn)
--
-- Return a topology.topology record from a layer
--
CREATE OR REPLACE FUNCTION topology.FindTopology(regclass, name)
RETURNS topology.topology
AS $$
    SELECT t.*
    FROM topology.topology t
    JOIN topology.layer l ON (t.id = l.topology_id)
    WHERE format('%I.%I', l.schema_name, l.table_name)::regclass = $1
    AND l.feature_column = $2;
$$ LANGUAGE 'sql';
--}

--{
--  FindTopology(layerSchema, layerTable, layerColumn)
--
-- Return a topology.topology record from a layer
--
CREATE OR REPLACE FUNCTION topology.FindTopology(name, name, name)
RETURNS topology.topology
AS $$
    SELECT t.*
    FROM topology.topology t
    JOIN topology.layer l ON (t.id = l.topology_id)
    WHERE l.schema_name = $1
    AND l.table_name = $2
    AND l.feature_column = $3;
$$ LANGUAGE 'sql';
--}

--{
--  FindTopology(topoName)
--
-- Return a topology.topology record from a topology name
--
CREATE OR REPLACE FUNCTION topology.FindTopology(text)
RETURNS topology.topology
AS $$
    SELECT *
    FROM topology.topology
    WHERE name = $1
$$ LANGUAGE 'sql';
--}

--{
--  FindTopology(topoId)
--
-- Return a topology.topology record from a topology id
--
CREATE OR REPLACE FUNCTION topology.FindTopology(integer)
RETURNS topology.topology
AS $$
    SELECT *
    FROM topology.topology
    WHERE id = $1
$$ LANGUAGE 'sql';
--}
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2021 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--{
--  FindLayer(TopoGeometry)
--
-- Return a topology.layer record from a TopoGeometry
--
CREATE OR REPLACE FUNCTION topology.FindLayer(tg topology.TopoGeometry)
RETURNS topology.layer
AS $$
    SELECT * FROM topology.layer
    WHERE topology_id = topology_id($1)
    AND layer_id = layer_id($1)
$$ LANGUAGE 'sql';
--}

--{
--  FindLayer(layerRegclass, layerColumn)
--
-- Return a topology.layer record from a layer table/column
--
CREATE OR REPLACE FUNCTION topology.FindLayer(layer_table regclass, feature_column name)
RETURNS topology.layer
AS $$
    SELECT l.*
    FROM topology.layer l, pg_class c, pg_namespace n
    WHERE l.schema_name = n.nspname
    AND l.table_name = c.relname
    AND c.oid = $1
    AND c.relnamespace = n.oid
    AND l.feature_column = $2
$$ LANGUAGE 'sql';
--}

--{
--  FindLayer(layerSchema, layerTable, layerColumn)
--
-- Return a topology.layer record from a layer schema/table/column
--
CREATE OR REPLACE FUNCTION topology.FindLayer(schema_name name, table_name name, feature_column name)
RETURNS topology.layer
AS $$
    SELECT * FROM topology.layer
    WHERE schema_name = $1
    AND table_name = $2
    AND feature_column = $3;
$$ LANGUAGE 'sql';
--}

--{
--  FindLayer(topoName)
--
-- Return a topology.layer record from a topology id and layer id
--
CREATE OR REPLACE FUNCTION topology.FindLayer(topology_id integer, layer_id integer)
RETURNS topology.layer
AS $$
    SELECT * FROM topology.layer
    WHERE topology_id = $1
      AND layer_id = $2
$$ LANGUAGE 'sql';
--}
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2010, 2021 Sandro Santilli <strk@kbt.io>
-- Copyright (C) 2005 Refractions Research Inc.
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- Author: Sandro Santilli <strk@kbt.io>
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

--
-- Type returned by ValidateTopology
--
CREATE TYPE topology.ValidateTopology_ReturnType AS (
  error varchar,
  id1 integer,
  id2 integer
);

--{
-- Return the exterior ring of a topology face
--
--
CREATE OR REPLACE FUNCTION topology._ValidateTopologyGetFaceShellMaximalEdgeRing(atopology varchar, aface int)
RETURNS GEOMETRY AS
$BODY$
DECLARE
  sql TEXT;
  outsidePoint GEOMETRY;
  shell GEOMETRY;
BEGIN

  sql := format(
    $$
      SELECT
        ST_Translate(
          ST_StartPoint( ST_BoundingDiagonal(mbr) ),
          -1,
          -1
        )
      FROM %1$I.face
      WHERE face_id = $1
    $$,
    atopology
  );
  EXECUTE sql USING aface INTO outsidePoint;

  sql := format(
    $$
      WITH
      outside_point AS (
        SELECT ST_Translate(
          ST_StartPoint( ST_BoundingDiagonal(mbr) ),
          -1,
          -1
        )
        FROM %1$I.face
        WHERE face_id = $1
      ),
      leftmost_edge AS (
        SELECT
          CASE WHEN left_face = $1
          THEN
            edge_id
          ELSE
            -edge_id
          END ring_id
        FROM %1$I.edge
        WHERE left_face = $1 or right_face = $1
        ORDER BY
          geom <-> $2
        LIMIT 1
      ),
      edgering AS (
        SELECT *
        FROM
          GetRingEdges(
            %1$L,
            (SELECT ring_id FROM leftmost_edge)
          )
      )
      SELECT
        ST_MakeLine(
          CASE WHEN r.edge > 0 THEN
            e.geom
          ELSE
            ST_Reverse(e.geom)
          END
          ORDER BY r.sequence
        ) outerRing
      FROM edgering r, %1$I.edge e
      WHERE e.edge_id = abs(r.edge)
    $$,
    atopology
  );

  --RAISE DEBUG 'SQL: %', sql;

  EXECUTE sql USING aface, outsidePoint
  INTO shell;

  -- TODO: check if the ring is not closed

  shell := ST_MakePolygon(shell);

  RETURN shell;
END;
$BODY$ LANGUAGE 'plpgsql' STABLE; --}

-- Assumes search_path has topology schema first
--{
CREATE OR REPLACE FUNCTION topology._ValidateTopologyGetRingEdges(starting_edge int)
RETURNS int[]
AS
$BODY$
DECLARE
  ret int[];
BEGIN
  WITH RECURSIVE edgering AS (
    SELECT
      starting_edge as signed_edge_id,
      edge_id,
      next_left_edge,
      next_right_edge
    FROM edge_data
    WHERE edge_id = abs(starting_edge)
      UNION
    SELECT
      CASE WHEN p.signed_edge_id < 0 THEN
        p.next_right_edge
      ELSE
        p.next_left_edge
      END,
      e.edge_id,
      e.next_left_edge,
      e.next_right_edge
    FROM edge_data e, edgering p
    WHERE e.edge_id =
      CASE WHEN p.signed_edge_id < 0 THEN
        abs(p.next_right_edge)
      ELSE
        abs(p.next_left_edge)
      END
  )
  SELECT array_agg(signed_edge_id)
  FROM edgering
  INTO ret;

  RETURN ret;
END;
$BODY$ LANGUAGE 'plpgsql';
--}

CREATE OR REPLACE FUNCTION topology._CheckEdgeLinking(curedge_edge_id INT, prevedge_edge_id INT, prevedge_next_left_edge INT, prevedge_next_right_edge INT)
RETURNS topology.ValidateTopology_ReturnType
AS
$BODY$
DECLARE
  retrec topology.ValidateTopology_ReturnType;
BEGIN
  IF prevedge_edge_id > 0
  THEN -- previous was outgoing, this one should be next-right
    IF prevedge_next_right_edge != curedge_edge_id THEN
      retrec.error = 'invalid next_right_edge';
      retrec.id1 = abs(prevedge_edge_id);
      retrec.id2 = curedge_edge_id; -- we put the expected one here, for convenience
      RETURN retrec;
    END IF;
  ELSE -- previous was incoming, this one should be next-left
    IF prevedge_next_left_edge != curedge_edge_id THEN
      retrec.error = 'invalid next_left_edge';
      retrec.id1 = abs(prevedge_edge_id);
      retrec.id2 = curedge_edge_id; -- we put the expected one here, for convenience
      RETURN retrec;
    END IF;
  END IF;

  RETURN retrec;
END;
$BODY$
LANGUAGE 'plpgsql' IMMUTABLE STRICT;

--
-- Check that the edges incident to topology nodes
-- (as advertised by their start_node/end_node)
-- correctly link to the next incident node on each
-- side (CW and CCW)
--
-- NOTE: if start_node/end_node values are incorrect the behavior
--       of this function is undefined
--
-- NOTE: assumes search_path was set before calling this function
--
CREATE OR REPLACE FUNCTION topology._ValidateTopologyEdgeLinking(bbox geometry DEFAULT NULL)
RETURNS SETOF topology.ValidateTopology_ReturnType
AS --{
$BODY$
DECLARE
  retrec topology.ValidateTopology_ReturnType;
  rec RECORD;
  last_node_id int;
  last_node_first_edge RECORD;
  last_node_prev_edge RECORD;
BEGIN
  RAISE NOTICE 'Checking edge linking';
  -- NOTE: this check relies on correct start_node and end_node
  --       for edges, if those are not correct the results
  --       of this check do not make much sense.
  FOR rec IN --{
      WITH
      nodes AS (
        SELECT node_id
        FROM node
        WHERE containing_face IS NULL
        AND (
          bbox IS NULL
          OR geom && bbox
        )
      ),
      incident_edges AS (
        SELECT
          n.node_id,
          e.edge_id,
          e.start_node,
          e.end_node,
          e.next_left_edge,
          e.next_right_edge,
          ST_RemoveRepeatedPoints(e.geom) as edge_geom
        FROM edge_data e, nodes n
        WHERE e.start_node = n.node_id
        or e.end_node = n.node_id
      ),
      edge_star AS (
        SELECT
          node_id,
          edge_id,
          next_left_edge,
          next_right_edge,
          ST_Azimuth(ST_StartPoint(edge_geom), ST_PointN(edge_geom, 2)) as az
        FROM incident_edges
        WHERE start_node = node_id
          UNION ALL
        SELECT
          node_id,
          -edge_id,
          next_left_edge,
          next_right_edge,
          ST_Azimuth(ST_EndPoint(edge_geom), ST_PointN(edge_geom, ST_NumPoints(edge_geom)-1))
        FROM incident_edges
        WHERE end_node = node_id
      ),
      sequenced_edge_star AS (
        SELECT
          row_number() over (partition by node_id order by az, edge_id) seq,
          *
        FROM edge_star
      )
      SELECT * FROM sequenced_edge_star
      ORDER BY node_id, seq
  LOOP --}{
    IF last_node_id IS NULL OR last_node_id != rec.node_id
    THEN --{
      IF last_node_id IS NOT NULL
      THEN
        -- Check that last edge (CW from prev one) is correctly linked
        retrec := topology._CheckEdgeLinking(
          last_node_first_edge.edge_id,
          last_node_prev_edge.edge_id,
          last_node_prev_edge.next_left_edge,
          last_node_prev_edge.next_right_edge
        );
        IF retrec IS NOT NULL
        THEN
          RETURN NEXT retrec;
        END IF;
      END IF;
      last_node_id = rec.node_id;
      last_node_first_edge = rec;
    ELSE --}{
      -- Check that this edge (CW from last one) is correctly linked
      retrec := topology._CheckEdgeLinking(
        rec.edge_id,
        last_node_prev_edge.edge_id,
        last_node_prev_edge.next_left_edge,
        last_node_prev_edge.next_right_edge
      );
      IF retrec IS NOT NULL
      THEN
        RETURN NEXT retrec;
      END IF;
    END IF; --}
    last_node_prev_edge = rec;
  END LOOP; --}
  IF last_node_id IS NOT NULL THEN --{
    -- Check that last edge (CW from prev one) is correctly linked
    retrec := topology._CheckEdgeLinking(
      last_node_first_edge.edge_id,
      last_node_prev_edge.edge_id,
      last_node_prev_edge.next_left_edge,
      last_node_prev_edge.next_right_edge
      );
    IF retrec IS NOT NULL
    THEN
      RETURN NEXT retrec;
    END IF;
  END IF; --}


END;
$BODY$ --}
LANGUAGE 'plpgsql' VOLATILE;

--
-- Check that the edges forming all rings
-- (as advertised by their next_right_edge/next_left_edge)
-- consistently advertise the same face on the walking side
-- (CW or CCW)
--
-- NOTE: if next_right_edge/next_left_edge values are incorrect
--       the behavior of this function is undefined, use
--       _ValidateTopologyEdgeLinking to verify that
--
-- NOTE: assumes search_path was set before calling this function
--
-- Creates a pg_temp.hole_check table
-- Creates a pg_temp.shell_check table
--
CREATE OR REPLACE FUNCTION topology._ValidateTopologyRings(bbox geometry DEFAULT NULL)
RETURNS SETOF topology.ValidateTopology_ReturnType
AS --{
$BODY$
DECLARE
  retrec topology.ValidateTopology_ReturnType;
  rec RECORD;
  ring_poly GEOMETRY;
  is_shell BOOLEAN;
  found_rings INT := 0;
  found_shells INT := 0;
  found_holes INT := 0;
BEGIN

  CREATE TEMP TABLE shell_check (
    face_id int PRIMARY KEY,
    ring_geom geometry
  );

  CREATE TEMP TABLE hole_check (
    ring_id int,
    hole_mbr geometry, -- point
    hole_point geometry, -- point
    in_shell int
  );

  RAISE NOTICE 'Building edge rings';

  -- Find all rings that can be formed on both sides
  -- of selected edges
  FOR rec IN
    WITH --{
    considered_edges AS (
      SELECT e.* FROM edge_data e, node n
      WHERE
        ( e.start_node = n.node_id OR e.end_node = n.node_id )
        AND
        ( bbox IS NULL OR n.geom && bbox )
    ),
    forward_rings AS (
      SELECT topology._ValidateTopologyGetRingEdges(e.edge_id) edges
      FROM considered_edges e
    ),
    forward_rings_with_id AS (
      SELECT
        (select min(e) FROM unnest(edges) e) ring_id,
        *
      FROM forward_rings
    ),
    distinct_forward_rings AS (
      SELECT
        DISTINCT ON (ring_id)
        *
      FROM forward_rings_with_id
    ),
    backward_rings AS (
      SELECT topology._ValidateTopologyGetRingEdges(-e.edge_id) edges
      FROM considered_edges e
      WHERE -edge_id NOT IN (
        SELECT x FROM (
          SELECT unnest(edges) x
          FROM distinct_forward_rings
        ) foo
      )
    ),
    backward_rings_with_id AS (
      SELECT
        (select min(e) FROM unnest(edges) e) ring_id,
        *
      FROM backward_rings
    ),
    distinct_backward_rings AS (
      SELECT
        DISTINCT ON (ring_id)
        *
      FROM backward_rings_with_id
    ),
    all_rings AS (
      SELECT * FROM distinct_forward_rings
      UNION
      SELECT * FROM distinct_backward_rings
    ),
    all_rings_with_ring_ordinal_edge AS (
      SELECT
        r.ring_id,
        e.seq,
        e.edge signed_edge_id
      FROM all_rings r
      LEFT JOIN LATERAL unnest(r.edges) WITH ORDINALITY AS e(edge, seq)
      ON TRUE
    ),
    all_rings_with_ring_geom AS (
      SELECT
        r.ring_id,
        ST_MakeLine(
          CASE WHEN signed_edge_id > 0 THEN
            e.geom
          ELSE
            ST_Reverse(e.geom)
          END
           -- TODO: how to make sure rows are ordered ?
          ORDER BY seq
        ) geom,
        array_agg(
          DISTINCT
          CASE WHEN signed_edge_id > 0 THEN
            e.left_face
          ELSE
            e.right_face
          END
        ) side_faces,
        count(signed_edge_id) num_edges,
        count(distinct abs(signed_edge_id)) distinct_edges
      FROM
        all_rings_with_ring_ordinal_edge r,
        edge_data e
      WHERE e.edge_id = abs(r.signed_edge_id)
      GROUP BY ring_id
    ) --}{
    SELECT ring_id, geom as ring_geom, side_faces, distinct_edges, num_edges
    FROM all_rings_with_ring_geom
  LOOP --}{

    found_rings := found_rings + 1;


    -- Check that there's a single face advertised
    IF array_upper(rec.side_faces,1) != 1
    THEN --{

      retrec.error = 'mixed face labeling in ring';
      retrec.id1 = rec.ring_id;
      retrec.id2 = NULL;
      RETURN NEXT retrec;
      CONTINUE;

    END IF; --}

    --RAISE DEBUG 'Ring geom: %', ST_AsTexT(rec.ring_geom);
    --RAISE DEBUG 'Distinct edges: %', rec.distinct_edges;
    --RAISE DEBUG 'Num edges: %', rec.num_edges;

    IF NOT ST_Equals(
      ST_StartPoint(rec.ring_geom),
      ST_EndPoint(rec.ring_geom)
    )
    THEN --{
      -- This should have been reported before,
      -- on the edge linking check
      retrec.error = 'non-closed ring';
      retrec.id1 = rec.ring_id;
      retrec.id2 = NULL;
      RETURN NEXT retrec;
      CONTINUE;
    END IF; --}

    -- Ring is valid, save it.
    is_shell := false;
    IF ST_NPoints(rec.ring_geom) > 3 AND
       rec.num_edges != rec.distinct_edges * 2
    THEN
      ring_poly := ST_MakePolygon(rec.ring_geom);
      IF ST_IsPolygonCCW(ring_poly) THEN
        is_shell := true;
      END IF;
    END IF;

    IF is_shell THEN --{ It's a shell (CCW)
      -- Check that a single face is ever used
      --       for each distinct CCW ring (shell)
      BEGIN
        INSERT INTO shell_check VALUES (
          rec.side_faces[1],
          ring_poly
        );
        found_shells := found_shells + 1;
      EXCEPTION WHEN unique_violation THEN
        retrec.error = 'face has multiple shells';
        retrec.id1 = rec.side_faces[1];
        retrec.id2 = rec.ring_id;
        RETURN NEXT retrec;
      END;
    ELSE -- }{ It's an hole (CW)
    -- NOTE: multiple CW rings (holes) can exist for a given face
      INSERT INTO hole_check VALUES (
        rec.ring_id,
        ST_Envelope(rec.ring_geom),
        ST_PointN(rec.ring_geom, 1),
        -- NOTE: we don't incurr in the risk
        --       of a ring touching the shell
        --       because in those cases the
        --       intruding "hole" will not really
        --       be considered an hole as its ring
        --       will not be CW
        rec.side_faces[1]
      );
      found_holes := found_holes + 1;
    END IF; --} hole

  END LOOP; --}

  RAISE NOTICE 'Found % rings, % valid shells, % valid holes',
    found_rings, found_shells, found_holes
  ;


END;
$BODY$ --}
LANGUAGE 'plpgsql' VOLATILE;

--{
--  ValidateTopology(toponame, [bbox])
--
--  Return a Set of ValidateTopology_ReturnType containing
--  informations on all topology inconsistencies
--
-- Availability: 1.0.0
-- Changed: 3.2.0 - add bbox optional parameter
--
CREATE OR REPLACE FUNCTION topology.ValidateTopology(toponame varchar, bbox geometry DEFAULT NULL)
  RETURNS setof topology.ValidateTopology_ReturnType
AS
$$
DECLARE
  retrec topology.ValidateTopology_ReturnType;
  rec RECORD;
  rec2 RECORD;
  affected_rows integer;
  invalid_edges integer[];
  invalid_faces integer[];
  has_invalid_edge_linking BOOLEAN := false;
  has_invalid_rings BOOLEAN := false;
  search_path_backup text;
  containing_face integer;
BEGIN

  IF NOT EXISTS (
    SELECT oid
    FROM pg_catalog.pg_namespace
    WHERE nspname = toponame
  )
  THEN
    RAISE EXCEPTION 'Topology schema % does not exist', toponame;
  END IF;

  IF NOT EXISTS (
    SELECT id
    FROM topology.topology
    WHERE name = toponame
  )
  THEN
    RAISE WARNING 'Topology % is not registered in topology.topology', toponame;
  END IF;

  EXECUTE 'SHOW search_path' INTO search_path_backup;
  EXECUTE 'SET search_PATH TO ' || quote_ident(toponame) || ','
                                || search_path_backup;

  IF bbox IS NOT NULL THEN
    RAISE NOTICE 'Limiting topology checking to bbox %', ST_AsEWKT(ST_Envelope(bbox));
  END IF;


  -- Check for coincident nodes
  RAISE NOTICE 'Checking for coincident nodes';
  FOR rec IN
    SELECT a.node_id as id1, b.node_id as id2
    FROM
      node a,
      node b
    WHERE a.node_id < b.node_id
    AND ST_DWithin(a.geom, b.geom, 0) -- NOTE: see #1625 and #1789
    AND (
      bbox IS NULL
      OR (
        a.geom && bbox
        AND
        b.geom && bbox
      )
    )
  LOOP
    retrec.error = 'coincident nodes';
    retrec.id1 = rec.id1;
    retrec.id2 = rec.id2;
    RETURN NEXT retrec;
  END LOOP;

  -- Check for edge crossed nodes
  -- TODO: do this in the single edge loop
  RAISE NOTICE 'Checking for edges crossing nodes';
  FOR rec IN
    SELECT n.node_id as nid, e.edge_id as eid
    FROM
      node n,
      edge e
    WHERE e.start_node != n.node_id
    AND e.end_node != n.node_id
    AND ST_Within(n.geom, e.geom)
    AND (
      bbox IS NULL
      OR (
        n.geom && bbox
        AND
        e.geom && bbox
      )
    )
  LOOP
    retrec.error = 'edge crosses node';
    retrec.id1 = rec.eid; -- edge_id
    retrec.id2 = rec.nid; -- node_id
    RETURN NEXT retrec;
  END LOOP;

  -- Scan all edges
  RAISE NOTICE 'Checking for invalid or not-simple edges';
  FOR rec IN
    SELECT e.geom, e.edge_id as id1, e.left_face, e.right_face
    FROM edge e
    WHERE (
      bbox IS NULL
      OR e.geom && bbox
    )
    ORDER BY edge_id
  LOOP --{

    -- Any invalid edge becomes a cancer for higher level complexes
    IF NOT ST_IsValid(rec.geom) THEN

      retrec.error = 'invalid edge';
      retrec.id1 = rec.id1;
      retrec.id2 = NULL;
      RETURN NEXT retrec;
      invalid_edges := array_append(invalid_edges, rec.id1);

      IF invalid_faces IS NULL OR NOT rec.left_face = ANY ( invalid_faces )
      THEN
        invalid_faces := array_append(invalid_faces, rec.left_face);
      END IF;

      IF rec.right_face != rec.left_face AND ( invalid_faces IS NULL OR
            NOT rec.right_face = ANY ( invalid_faces ) )
      THEN
        invalid_faces := array_append(invalid_faces, rec.right_face);
      END IF;

      CONTINUE;

    END IF;

    -- Check edge being simple (ie: not self-intersecting)
    IF NOT ST_IsSimple(rec.geom) THEN
      retrec.error = 'edge not simple';
      retrec.id1 = rec.id1;
      retrec.id2 = NULL;
      RETURN NEXT retrec;
    END IF;

  END LOOP; --}

  -- Check for edge crossing
  RAISE NOTICE 'Checking for crossing edges';
  FOR rec IN
    SELECT
      e1.edge_id as id1,
      e2.edge_id as id2,
      e1.geom as g1,
      e2.geom as g2,
      ST_Relate(e1.geom, e2.geom) as im
    FROM
      edge e1,
      edge e2
    WHERE
      e1.edge_id < e2.edge_id
      AND e1.geom && e2.geom
      AND (
        invalid_edges IS NULL OR (
          NOT e1.edge_id = ANY (invalid_edges)
          AND
          NOT e2.edge_id = ANY (invalid_edges)
        )
      )
      AND (
        bbox IS NULL
        OR (
          e1.geom && bbox
          AND
          e2.geom && bbox
        )
      )
  LOOP --{

    IF ST_RelateMatch(rec.im, 'FF1F**1*2') THEN
      CONTINUE; -- no interior intersection

    --
    -- Closed lines have no boundary, so endpoint
    -- intersection would be considered interior
    -- See http://trac.osgeo.org/postgis/ticket/770
    -- See also full explanation in topology.AddEdge
    --

    ELSIF ST_RelateMatch(rec.im, 'FF10F01F2') THEN
      -- first line (g1) is open, second (g2) is closed
      -- first boundary has puntual intersection with second interior
      --
      -- compute intersection, check it equals second endpoint
      IF ST_Equals(ST_Intersection(rec.g2, rec.g1),
                   ST_StartPoint(rec.g2))
      THEN
        CONTINUE;
      END IF;

    ELSIF ST_RelateMatch(rec.im, 'F01FFF102') THEN
      -- second line (g2) is open, first (g1) is closed
      -- second boundary has puntual intersection with first interior
      --
      -- compute intersection, check it equals first endpoint
      IF ST_Equals(ST_Intersection(rec.g2, rec.g1),
                   ST_StartPoint(rec.g1))
      THEN
        CONTINUE;
      END IF;

    ELSIF 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
      -- and equals first StartPoint _and_ second StartPoint
      IF ST_Equals(ST_Intersection(rec.g1, rec.g2),
                   ST_StartPoint(rec.g1)) AND
         ST_Equals(ST_StartPoint(rec.g1), ST_StartPoint(rec.g2))
      THEN
        CONTINUE;
      END IF;

    END IF;

    retrec.error = 'edge crosses edge';
    retrec.id1 = rec.id1;
    retrec.id2 = rec.id2;
    RETURN NEXT retrec;
  END LOOP; --}

  -- Check for edge start_node geometry mis-match
  -- TODO: move this in the first edge table scan
  RAISE NOTICE 'Checking for edges start_node mismatch';
  FOR rec IN
    SELECT e.edge_id as id1, n.node_id as id2
    FROM
      edge e,
      node n
    WHERE e.start_node = n.node_id
    AND NOT ST_Equals(ST_StartPoint(e.geom), n.geom)
    AND (
      bbox IS NULL
      OR e.geom && bbox
    )
  LOOP --{
    retrec.error = 'edge start node geometry mis-match';
    retrec.id1 = rec.id1;
    retrec.id2 = rec.id2;
    RETURN NEXT retrec;
  END LOOP; --}

  -- Check for edge end_node geometry mis-match
  -- TODO: move this in the first edge table scan
  RAISE NOTICE 'Checking for edges end_node mismatch';
  FOR rec IN
    SELECT e.edge_id as id1, n.node_id as id2
    FROM
      edge e,
      node n
    WHERE e.end_node = n.node_id
    AND NOT ST_Equals(ST_EndPoint(e.geom), n.geom)
    AND (
      bbox IS NULL
      OR e.geom && bbox
    )
  LOOP --{
    retrec.error = 'edge end node geometry mis-match';
    retrec.id1 = rec.id1;
    retrec.id2 = rec.id2;
    RETURN NEXT retrec;
  END LOOP; --}

  -- Check for faces w/out edges
  RAISE NOTICE 'Checking for faces without edges';
  FOR rec IN
    SELECT face_id as id1
    FROM face
    WHERE face_id > 0
    AND (
      bbox IS NULL
      OR mbr && bbox
    )
    EXCEPT (
      SELECT left_face FROM edge
      UNION
      SELECT right_face FROM edge
    )
  LOOP --{
    retrec.error = 'face without edges';
    retrec.id1 = rec.id1;
    retrec.id2 = NULL;
    RETURN NEXT retrec;
  END LOOP; --}

  -- Validate edge linking
  -- NOTE: relies on correct start_node/end_node on edges
  FOR rec IN SELECT * FROM topology._ValidateTopologyEdgeLinking(bbox)
  LOOP
    RETURN next rec;
    has_invalid_edge_linking := true;
  END LOOP;

  IF has_invalid_edge_linking THEN
    DROP TABLE IF EXISTS pg_temp.hole_check;
    DROP TABLE IF EXISTS pg_temp.shell_check;
    RETURN; -- does not make sense to continue
  END IF;

  --- Validate edge rings
  FOR rec IN SELECT * FROM topology._ValidateTopologyRings(bbox)
  LOOP
    RETURN next rec;
    has_invalid_rings := true;
  END LOOP;

  IF has_invalid_rings THEN
    DROP TABLE IF EXISTS pg_temp.hole_check;
    DROP TABLE IF EXISTS pg_temp.shell_check;
    RETURN; -- does not make sense to continue
  END IF;

  -- Now create a temporary table to construct all face geometries
  -- for checking their consistency

  RAISE NOTICE 'Constructing geometry of all faces';
  -- TODO: only construct exterior ring

  CREATE TEMP TABLE face_check ON COMMIT DROP AS
  SELECT
    sc.face_id,
    sc.ring_geom AS shell,
    f.mbr
  FROM
    pg_temp.shell_check sc, face f
  WHERE
    f.face_id = sc.face_id
  ;

  DROP TABLE pg_temp.shell_check;


  IF bbox IS NOT NULL
  THEN --{
    INSERT INTO pg_temp.face_check
    SELECT face_id,
      topology._ValidateTopologyGetFaceShellMaximalEdgeRing(toponame, face_id),
      mbr
    FROM face
    WHERE mbr && bbox
    AND (
      CASE WHEN invalid_faces IS NOT NULL THEN
        NOT face_id = ANY(invalid_faces)
      ELSE
        TRUE
      END
    )
    AND face_id NOT IN (
      SELECT face_id FROM pg_temp.face_check
    )
    ;
  END IF; --}

  -- Build a gist index on geom
  CREATE INDEX ON face_check USING gist (shell);

  -- Build a btree index on id
  CREATE INDEX ON face_check (face_id);

  -- Scan the table looking for NULL geometries
  -- or geometries with wrong MBR consistency
  RAISE NOTICE 'Checking faces';
  affected_rows := 0;
  FOR rec IN
    SELECT * FROM face_check
  LOOP --{

    affected_rows := affected_rows + 1;

    IF rec.shell IS NULL OR ST_IsEmpty(rec.shell)
    THEN
      -- Face missing !
      retrec.error := 'face has no rings';
      retrec.id1 := rec.face_id;
      retrec.id2 := NULL;
      RETURN NEXT retrec;
    END IF;

    IF NOT ST_Equals(rec.mbr, ST_Envelope(rec.shell))
    THEN
      -- Inconsistent MBR!
      retrec.error := 'face has wrong mbr';
      retrec.id1 := rec.face_id;
      retrec.id2 := NULL;
      RETURN NEXT retrec;
    END IF;

  END LOOP; --}

  RAISE NOTICE 'Checked % faces', affected_rows;

  -- Check edges are covered by their left-right faces (#4830)
  RAISE NOTICE 'Checking for holes coverage';
  affected_rows := 0;
  FOR rec IN
    SELECT * FROM hole_check
  LOOP --{
    SELECT f.face_id
    FROM face_check f
    WHERE rec.hole_mbr @ f.shell
    AND _ST_Contains(f.shell, rec.hole_point)
    ORDER BY ST_Area(f.shell) ASC
    LIMIT 1
    INTO rec2;

    IF ( NOT FOUND AND rec.in_shell != 0 )
       OR ( rec2.face_id != rec.in_shell )
    THEN
        retrec.error := 'hole not in advertised face';
        retrec.id1 := rec.ring_id;
        retrec.id2 := NULL;
        RETURN NEXT retrec;
    END IF;
    affected_rows := affected_rows + 1;

  END LOOP; --}

  RAISE NOTICE 'Finished checking for coverage of % holes', affected_rows;

  -- Check nodes have correct containing_face (#3233)
  -- NOTE: relies on correct edge linking
  RAISE NOTICE 'Checking for node containing_face correctness';
  FOR rec IN
    SELECT
      n.node_id,
      n.geom geom,
      n.containing_face,
      e.edge_id
    FROM node n
    LEFT JOIN edge e ON (
      e.start_node = n.node_id OR
      e.end_node = n.node_id
    )
    WHERE
     ( bbox IS NULL OR n.geom && bbox )
  LOOP --{

    IF rec.edge_id IS NOT NULL
    THEN --{
      -- Node is not isolated, make sure it
      -- advertises itself as such
      IF rec.containing_face IS NOT NULL
      THEN --{
        -- node is not really isolated
        retrec.error := 'not-isolated node has not-null containing_face';
        retrec.id1 := rec.node_id;
        retrec.id2 := NULL;
        RETURN NEXT retrec;
      END IF; --}
    ELSE -- }{
      -- Node is isolated, make sure it
      -- advertises itself as such
      IF rec.containing_face IS NULL
      THEN --{
        -- isolated node advertises itself as non-isolated
        retrec.error := 'isolated node has null containing_face';
        retrec.id1 := rec.node_id;
        retrec.id2 := NULL;
        RETURN NEXT retrec;
      ELSE -- }{
        -- node is isolated and advertising a containing_face
        -- now let's check it's really in contained by it
        BEGIN
          containing_face := topology.GetFaceContainingPoint(toponame, rec.geom);
        EXCEPTION WHEN OTHERS THEN
          RAISE NOTICE 'Got % (%)', SQLSTATE, SQLERRM;
          retrec.error := format('got exception trying to find face containing node: %s', SQLERRM);
          retrec.id1 := rec.node_id;
          retrec.id2 := NULL;
          RETURN NEXT retrec;
        END;
        IF containing_face != rec.containing_face THEN
          retrec.error := 'isolated node has wrong containing_face';
          retrec.id1 := rec.node_id;
          retrec.id2 := NULL; -- TODO: write expected containing_face here ?
          RETURN NEXT retrec;
        END IF;
      END IF; --}
    END IF; --}

  END LOOP; --}


  DROP TABLE pg_temp.hole_check;
  DROP TABLE pg_temp.face_check;

  EXECUTE 'SET search_PATH TO ' || search_path_backup;

  RETURN;
END
$$
LANGUAGE 'plpgsql' VOLATILE; -- NOTE: we need VOLATILE to use SHOW
--} ValidateTopology(toponame, bbox)

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.net
--
-- Copyright (C) 2021 Sandro Santilli <strk@kbt.io>
--
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
--
-- Author: Sandro Santilli <strk@kbt.io>
--
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-- Invalidities are:
--   - Referenced element does not exist
--   - Same element listed twice
--
CREATE OR REPLACE FUNCTION topology.ValidateTopologyRelation(toponame varchar)
RETURNS TABLE(error text, layer_id int, topogeo_id int, element_id int)
AS
$BODY$
DECLARE
  layerrec RECORD;
  rel RECORD;
  search_path_backup text;
BEGIN
  IF NOT EXISTS (
    SELECT oid
    FROM pg_catalog.pg_namespace
    WHERE nspname = toponame
  )
  THEN
    RAISE EXCEPTION 'Topology schema % does not exist', toponame;
  END IF;

  IF NOT EXISTS (
    SELECT id
    FROM topology.topology
    WHERE name = toponame
  )
  THEN
    RAISE WARNING 'Topology % is not registered in topology.topology', toponame;
  END IF;

  EXECUTE 'SHOW search_path' INTO search_path_backup;
  EXECUTE 'SET search_PATH TO ' || quote_ident(toponame) || ','
                                || search_path_backup;

  FOR layerrec IN SELECT * FROM topology.layer
  LOOP --{
    IF layerrec.child_id IS NULL
    THEN --{ Layer is simple

      -- Check that all referenced nodes exist
      FOR rel IN
        SELECT r.layer_id, r.topogeo_id, r.element_id
        FROM relation r
        WHERE r.layer_id = layerrec.layer_id
        AND r.element_type = 1
        AND r.element_id NOT IN (
          SELECT node_id FROM node
        )
      LOOP
        error := 'TopoGeometry references unexistent node';
        layer_id := rel.layer_id;
        topogeo_id := rel.topogeo_id;
        element_id := rel.element_id;
        RETURN NEXT;
      END LOOP;

      -- Check that all referenced edges exist
      FOR rel IN
        SELECT r.layer_id, r.topogeo_id, r.element_id
        FROM relation r
        WHERE r.layer_id = layerrec.layer_id
        AND r.element_type = 2
        AND abs(r.element_id) NOT IN (
          SELECT edge_id FROM edge_data
        )
      LOOP
        error := 'TopoGeometry references unexistent edge';
        layer_id := rel.layer_id;
        topogeo_id := rel.topogeo_id;
        element_id := rel.element_id;
        RETURN NEXT;
      END LOOP;

      -- Check that all referenced faces exist
      FOR rel IN
        SELECT r.layer_id, r.topogeo_id, r.element_id
        FROM relation r
        WHERE r.layer_id = layerrec.layer_id
        AND r.element_type = 3
        AND r.element_id NOT IN (
          SELECT face_id FROM face
        )
      LOOP
        error := 'TopoGeometry references unexistent face';
        layer_id := rel.layer_id;
        topogeo_id := rel.topogeo_id;
        element_id := rel.element_id;
        RETURN NEXT;
      END LOOP;

    ELSE -- }{ Layer is hierarchical

      --RAISE DEBUG 'Checking hierarchical layer %', layerrec.layer_id;

      FOR rel IN
        SELECT r.layer_id, r.topogeo_id, r.element_id
        FROM relation r
        WHERE r.layer_id = layerrec.layer_id
        AND abs(r.element_id) NOT IN (
          SELECT r2.topogeo_id
          FROM relation r2
          WHERE r2.layer_id = layerrec.child_id
        )
      LOOP
        error := 'TopoGeometry references unexistent child';
        layer_id := rel.layer_id;
        topogeo_id := rel.topogeo_id;
        element_id := rel.element_id;
        RETURN NEXT;
      END LOOP;

    END IF; --} Layer is hierarchical
  END LOOP; --}

  EXECUTE 'SET search_PATH TO ' || search_path_backup;
END;
$BODY$
LANGUAGE 'plpgsql' VOLATILE STRICT;

-- Cleanup functions
-- This function finds primitives intersecting a bounding box
-- that are not used in the definition of TopoGeometry objects
-- and remove them all.
--
-- Returns number of deleted primitives
--
CREATE OR REPLACE FUNCTION topology.RemoveUnusedPrimitives(
  atopology text,
  bbox GEOMETRY DEFAULT NULL
)
RETURNS INT AS
$BODY$
DECLARE
  topo topology.topology;
  deletedNodes INT := 0;
  deletedEdges INT := 0;
  deletedNodesDeg2 INT := 0;
  sql TEXT;
  rec RECORD;
  edgeMap JSONB := '{}';
  edge1 INT;
  edge2 INT;
  removedNode INT;
  ok BOOLEAN;
  fixedLinks INT := 0;
  mergedFaces INT[];
  moreMergedFaces INT[];
BEGIN

  topo := findTopology(atopology);
  IF topo.id IS NULL THEN
    RAISE EXCEPTION 'Could not find topology "%"', atopology;
  END IF;

  RAISE NOTICE 'Removing unused edges';

  RAISE DEBUG 'Determining edges not referenced by linear TopoGeoms';
  -- Delete edges not used in non-hierarchical TopoGeometry
  -- from linear typed layers
  sql := format(
    $$
      CREATE TEMPORARY TABLE deleted_edges AS
      SELECT
        edge_id,
        next_right_edge,
        next_left_edge,
        left_face,
        right_face,
        start_node,
        end_node
      FROM %1$I.edge_data e
      WHERE ( $1 IS NULL OR ST_Intersects(geom, $1) )
      AND NOT EXISTS (
        SELECT 1
        FROM %1$I.relation r, topology.layer l
        WHERE r.layer_id = l.layer_id
        AND l.topology_id = $2
        AND l.child_id IS NULL
        AND l.feature_type IN (2, 4)
        AND r.element_id in ( e.edge_id, -e.edge_id )
      )
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  EXECUTE sql USING bbox, topo.id;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Found % edges not referenced by linear TopoGeoms', fixedLinks;

  -- remove from deleted_edges the edges binding
  -- faces that individually (not both) take part
  -- of the definition of an areal TopoGeometry
  RAISE DEBUG 'Determining edges not binding areal TopoGeoms';
  sql := format(
    $$
      WITH breaking_merges AS (
        SELECT
          DISTINCT
          de.edge_id
          --, ARRAY[r.layer_id, r.topogeo_id] topogeo
          --, array_agg(r.element_id) faces
        FROM
          topology.layer l,
          %1$I.relation r,
          pg_temp.deleted_edges de
        WHERE l.topology_id = %2$L
        AND l.child_id IS NULL -- non-hierarchical layer
        AND l.feature_type IN (3, 4) -- areal or mixed layer
        AND r.layer_id = l.layer_id
        AND r.element_type = 3 -- face primitive
        AND de.left_face != de.right_face -- non-dangling edges
        AND ( r.element_id = de.left_face OR r.element_id = de.right_face )
        GROUP BY de.edge_id, r.layer_id, r.topogeo_id
        HAVING count(DISTINCT r.element_id) != 2
      )
      --SELECT * FROM breaking_merges
      DELETE FROM pg_temp.deleted_edges de
      WHERE edge_id IN (
        SELECT edge_id FROM breaking_merges
      )
    $$,
    topo.name,
    topo.id
  );
  --RAISE DEBUG 'SQL: %', sql;
  EXECUTE sql;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Retained % edges binding areal TopoGeoms', fixedLinks;
  --FOR rec IN EXECUTE sql LOOP
    --RAISE NOTICE 'Should retain edges % binding areal TopoGeom % in layer %', rec.edge_id, rec.topogeo_id, rec.layer_id;
  --END LOOP;


  RAISE DEBUG 'Deleting unused edges';
  sql := format(
    $$
      DELETE FROM %1$I.edge_data e
      WHERE e.edge_id IN (
        SELECT edge_id FROM pg_temp.deleted_edges
      )
    $$,
    topo.name
  );
  EXECUTE sql;
  GET DIAGNOSTICS deletedEdges = ROW_COUNT;
  RAISE DEBUG 'Deleted % unused edges', deletedEdges;


  RAISE DEBUG 'Fixing broken next_right_edge links';
  sql := format(
    $$
      UPDATE %1$I.edge_data e
      SET
        next_right_edge =
          CASE
          WHEN e.next_right_edge = ne.edge_id THEN
            ne.next_right_edge
          ELSE
            ne.next_left_edge
          END,
        abs_next_right_edge =
          CASE
          WHEN e.next_right_edge = ne.edge_id THEN
            abs(ne.next_right_edge)
          ELSE
            abs(ne.next_left_edge)
          END
      FROM pg_temp.deleted_edges ne
      WHERE e.abs_next_right_edge = ne.edge_id
        AND e.next_right_edge !=
        CASE
        WHEN e.next_right_edge = ne.edge_id THEN
          ne.next_right_edge
        ELSE
          ne.next_left_edge
        END
      RETURNING e.*
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  LOOP
    EXECUTE sql;
    GET DIAGNOSTICS fixedLinks = ROW_COUNT;
--    fixedLinks := 0;
--    FOR rec IN EXECUTE sql LOOP
--      fixedLinks := fixedLinks + 1;
--      RAISE DEBUG 'Updated next_right_edge link for edge %, now having next_right_edge=% and abs_next_right_edge=%', rec.edge_id, rec.next_right_edge, rec.abs_next_right_edge;
--    END LOOP;
    IF fixedLinks = 0 THEN
      RAISE DEBUG 'No (more) broken next_right_edge links';
      EXIT;
    END IF;
    RAISE DEBUG 'Updated % broken next_right_edge links', fixedLinks;
  END LOOP;

  RAISE DEBUG 'Fixing broken next_left_edge links';
  sql := format(
    $$
      UPDATE %1$I.edge_data e
      SET
        next_left_edge =
          CASE
          WHEN e.next_left_edge = ne.edge_id THEN
            ne.next_right_edge
          ELSE
            ne.next_left_edge
          END,
        abs_next_left_edge =
          CASE
          WHEN e.next_left_edge = ne.edge_id THEN
            abs(ne.next_right_edge)
          ELSE
            abs(ne.next_left_edge)
          END
      FROM pg_temp.deleted_edges ne
      WHERE e.abs_next_left_edge = ne.edge_id
        -- Avoid updating records which do not need
        -- to be updated (alternatively we could DELETE
        -- those records from deleted_edges before next iteration)
        AND e.next_left_edge !=
        CASE
        WHEN e.next_left_edge = ne.edge_id THEN
          ne.next_right_edge
        ELSE
          ne.next_left_edge
        END
      RETURNING e.*
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  LOOP
    EXECUTE sql;
    GET DIAGNOSTICS fixedLinks = ROW_COUNT;
--    fixedLinks := 0;
--    FOR rec IN EXECUTE sql LOOP
--      fixedLinks := fixedLinks + 1;
--      RAISE DEBUG 'Updated next_left_edge link for edge %, now having next_left_edge=% and abs_next_left_edge=%', rec.edge_id, rec.next_left_edge, rec.abs_next_left_edge;
--    END LOOP;
    IF fixedLinks = 0 THEN
      RAISE DEBUG 'No (more) broken next_left_edge links found';
      EXIT;
    END IF;
    RAISE DEBUG 'Updated % broken next_left_edge links', fixedLinks;
  END LOOP;

  --
  -- Build arrays of faces to be merged
  --

  RAISE DEBUG 'Building face merge sets';

  CREATE TEMPORARY TABLE mergeable_faces AS
  WITH merges AS (
    SELECT
      DISTINCT
      ARRAY[
        LEAST(left_face, right_face),
        GREATEST(left_face, right_face)
      ] faceset
    FROM deleted_edges
    WHERE left_face != right_face
  )
  SELECT faceset
  FROM merges;

  CREATE TEMPORARY TABLE merged_faces (keep INT, merge INT[]);

  LOOP -- {

    -- Fetch next merge
    DELETE FROM mergeable_faces
    WHERE ctid = (SELECT ctid FROM mergeable_faces LIMIT 1)
    RETURNING faceset
    INTO mergedFaces;
    IF mergedFaces IS NULL THEN
      EXIT;
    END IF;

    RAISE DEBUG 'Next merged faces start with: %', mergedFaces;

    LOOP --{
      WITH deleted AS (
        DELETE FROM mergeable_faces
        WHERE faceset && mergedFaces
        RETURNING faceset
      ), flood_faces AS (
        SELECT DISTINCT unnest(faceset) merged
        FROM deleted
      )
      SELECT array_agg(merged)
      FROM flood_faces
      INTO moreMergedFaces;

      IF moreMergedFaces IS NULL THEN
        EXIT;
      END IF;

      RAISE DEBUG 'There is more merged faces: %', moreMergedFaces;
      SELECT array_agg(x) FROM (
        SELECT unnest(mergedFaces) x
          UNION
        SELECT unnest(moreMergedFaces)
      ) foo
      INTO mergedFaces;
      RAISE DEBUG 'Merged faces grows to: %', mergedFaces;

    END LOOP; --}

    mergedFaces := array_agg(distinct x ORDER BY x) FROM unnest(mergedFaces) x;
    RAISE DEBUG 'Storing merged faceset: %', mergedFaces;

    INSERT INTO pg_temp.merged_faces VALUES (
      mergedFaces[1],
      array_remove(mergedFaces, mergedFaces[1])
    );

  END LOOP; --}

  DROP TABLE pg_temp.mergeable_faces;

  --
  -- Fix face labels
  --

  RAISE DEBUG 'Fixing broken left_face labels';
  sql := format(
    $$
      UPDATE %1$I.edge_data e
      SET left_face = mf.keep
      FROM pg_temp.merged_faces mf
      WHERE e.left_face != mf.keep
        AND e.left_face = ANY(mf.merge)
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  EXECUTE sql;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Updated % broken left_face links', fixedLinks;

  RAISE DEBUG 'Fixing broken right_face labels';
  sql := format(
    $$
      UPDATE %1$I.edge_data e
      SET right_face = mf.keep
      FROM pg_temp.merged_faces mf
      WHERE e.right_face != mf.keep
        AND e.right_face = ANY(mf.merge)
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  EXECUTE sql;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Updated % broken right_face links', fixedLinks;


  RAISE DEBUG 'Updating containing_face labels for merged faces';
  sql := format(
    $$
      UPDATE %1$I.node n
      SET containing_face = mf.keep
      FROM pg_temp.merged_faces mf
      WHERE n.containing_face != mf.keep
      AND n.containing_face = ANY(mf.merge)
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  EXECUTE sql;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Updated % containing_face labels for nodes', fixedLinks;

  --
  -- Fix face table (delete/update mbr)
  --

  RAISE DEBUG 'Updating merged faces MBR';
  sql := format(
    $$
      WITH merged_mbr AS (
        SELECT
          mf.keep,
          ST_Envelope(
            ST_Collect(
              f.mbr
            )
          ) mbr
        FROM pg_temp.merged_faces mf
        JOIN %1$I.face f ON (
          f.face_id = mf.keep OR
          f.face_id = ANY( mf.merge )
        )
        WHERE mf.keep != 0
        GROUP by mf.keep
      )
      UPDATE %1$I.face f
      SET mbr = m.mbr
      FROM merged_mbr m
      WHERE f.face_id = m.keep
    $$,
    topo.name
  );
  EXECUTE sql;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Updated % merged faces MBR', fixedLinks;

  RAISE DEBUG 'Deleting removed faces';
  sql := format(
    $$
      DELETE FROM %1$I.face
      USING pg_temp.merged_faces mf
      WHERE face_id = ANY (mf.merge)
    $$,
    topo.name
  );
  EXECUTE sql;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Deleted % merged faces', fixedLinks;

  --
  -- Fix TopoGeometry
  --

  RAISE DEBUG 'Updating areal TopoGeometry definitions';
  -- We remove the merged faces from the definition
  -- of areal TopoGeometry objects
  sql := format(
    $$
      WITH deleted AS (
        DELETE FROM %1$I.relation r
        USING topology.layer l, pg_temp.merged_faces mf
        WHERE l.topology_id = %2$L
        AND l.feature_type IN (3, 4)
        AND l.child_id IS NULL
        AND r.layer_id = l.layer_id
        AND r.element_id = ANY (mf.merge)
        RETURNING
          r.topogeo_id,
          r.layer_id,
          l.schema_name,
          l.table_name,
          l.feature_column,
          mf.merge,
          mf.keep,
          r.element_id
      )
      SELECT
        topogeo_id,
        layer_id,
        schema_name,
        table_name,
        feature_column,
        merge,
        keep,
        array_agg(element_id) lost_faces
      FROM deleted
      GROUP BY 1,2,3,4,5,6,7
    $$,
    topo.name,
    topo.id
  );
  --RAISE NOTICE 'SQL: %', sql;
  FOR rec IN EXECUTE sql
  LOOP
    RAISE DEBUG 'Areal TopoGeometry % in layer %.%.% '
      'lost faces % (kept %) in its composition',
      rec.topogeo_id, rec.schema_name,
      rec.table_name, rec.feature_column,
      rec.lost_faces, rec.keep
    ;
  END LOOP;

  --
  -- Mark newly isolated nodes as such
  --

  RAISE DEBUG 'Determining newly isolated nodes';
  sql := format(
    $$
      WITH unlinked_nodes AS (
        SELECT start_node node_id FROM pg_temp.deleted_edges
          UNION
        SELECT end_node FROM pg_temp.deleted_edges
      ), isolated AS (
        SELECT node_id FROM unlinked_nodes
        EXCEPT SELECT start_node FROM %1$I.edge_data
        EXCEPT  SELECT end_node FROM %1$I.edge_data
      ), incident_faces AS (
        SELECT
          node_id,
          array_agg(DISTINCT face_id) incident_faces
        FROM (
          SELECT DISTINCT node_id, unnest(face_id) face_id
          FROM (
            SELECT
              i.node_id,
              ARRAY[e.left_face, e.right_face] face_id
            FROM isolated i, deleted_edges e
            WHERE e.start_node = i.node_id
              UNION
            SELECT
              i.node_id,
              ARRAY[e.left_face, e.right_face] face_id
            FROM isolated i, deleted_edges e
            WHERE e.end_node = i.node_id
          ) foo
        ) bar
        GROUP BY node_id
      ), containing_faces AS (
        SELECT
          inc.node_id,
          COALESCE(mf.keep, incident_faces[1]) face_id
        FROM incident_faces inc
        LEFT JOIN pg_temp.merged_faces mf
          ON ( inc.incident_faces && mf.merge )
      )
      UPDATE %1$I.node n
      SET containing_face = cf.face_id
      FROM containing_faces cf
      WHERE n.node_id = cf.node_id
      AND n.containing_face IS DISTINCT FROM cf.face_id
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  EXECUTE sql;
  GET DIAGNOSTICS fixedLinks = ROW_COUNT;
  RAISE DEBUG 'Isolated % nodes', fixedLinks;

  RAISE NOTICE 'Removed % unused edges', deletedEdges;


  --
  -- Clean isolated nodes
  --

  -- Cleanup isolated nodes
  -- (non-isolated ones would have become isolated by now)
  sql := format(
    $$
      SELECT
        n.node_id
      FROM
        %1$I.node n
      WHERE ( $1 IS NULL OR ST_Intersects(n.geom, $1) )
      AND n.containing_face IS NOT NULL
      AND NOT EXISTS (
        SELECT 1
        FROM %1$I.relation r, topology.layer l
        WHERE r.layer_id = l.layer_id
        AND l.topology_id = $2
        AND l.child_id IS NULL
        AND l.feature_type = 1
        AND r.element_id = n.node_id
      )
    $$,
    topo.name
  );
  RAISE NOTICE 'Removing isolated nodes';
  FOR rec in EXECUTE sql USING bbox, topo.id
  LOOP --{
    BEGIN
      PERFORM topology.ST_RemIsoNode(topo.name, rec.node_id);
      RAISE DEBUG 'Removed isolated node %', rec.node_id;
      deletedNodes := deletedNodes + 1;
    EXCEPTION WHEN OTHERS
    THEN
      RAISE WARNING 'Isolated node % could not be removed: %', rec.node_id, SQLERRM;
    END;
  END LOOP; --}
  RAISE NOTICE 'Removed % isolated nodes', deletedNodes;

  -- Remove nodes connecting only 2 edges if
  -- no lineal TopoGeometry exists that is defined
  -- by only one of them
  sql := format(
    $$
      WITH
      unused_connected_nodes_in_bbox AS (
        SELECT
          n.node_id
        FROM %1$I.node n
        WHERE
          ( $1 IS NULL OR ST_Intersects(n.geom, $1) )
          AND n.containing_face IS NULL

          EXCEPT

        SELECT r.element_id
        FROM %1$I.relation r
        JOIN topology.layer l ON ( r.layer_id = l.layer_id )
        WHERE l.child_id IS NULL
        AND l.topology_id = $2
        AND l.feature_type IN ( 1, 4 ) -- puntual or mixed layer
        AND r.element_type = 1 -- node primitive
      ),
      removable_nodes_of_degree_2_in_bbox AS (
        SELECT
          n.node_id,
          array_agg(e.edge_id) edges
        FROM
          unused_connected_nodes_in_bbox n,
          %1$I.edge e
        WHERE (
            n.node_id = e.start_node
            OR n.node_id = e.end_node
          )
        GROUP BY n.node_id
        HAVING count(e.edge_id) = 2
      ),
      breaking_heals AS (
        SELECT
          DISTINCT
          n.node_id
          -- , ARRAY[r.layer_id, r.topogeo_id]
          -- , array_agg(r.element_id) edges
        FROM
          removable_nodes_of_degree_2_in_bbox n,
          %1$I.relation r,
          topology.layer l
        WHERE l.topology_id = $2
          AND l.child_id IS NULL
          AND l.feature_type IN ( 2, 4 ) -- lineal or mixed layer
          AND r.layer_id = l.layer_id
          AND r.element_type = 2 -- edge primitive
          AND r.element_id IN (
            n.edges[1], -n.edges[1],
            n.edges[2], -n.edges[2]
          )
          GROUP BY n.node_id, r.layer_id, r.topogeo_id
          HAVING count(DISTINCT abs(r.element_id)) != 2
      )
      SELECT
        node_id,
        edges[1] edge1,
        edges[2] edge2
      FROM removable_nodes_of_degree_2_in_bbox
      WHERE node_id NOT IN (
        SELECT node_id FROM breaking_heals
      )
    $$,
    topo.name
  );
  --RAISE DEBUG 'SQL: %', sql;
  RAISE NOTICE 'Removing unneeded nodes of degree 2';
  EXECUTE sql USING bbox, topo.id;
  FOR rec in EXECUTE sql USING bbox, topo.id
  LOOP --{
    RAISE DEBUG 'edgeMap: %', edgeMap;
    -- Edges may have changed name
    edge1 := COALESCE( (edgeMap -> rec.edge1::text)::int, rec.edge1);
    edge2 := COALESCE( (edgeMap -> rec.edge2::text)::int, rec.edge2);

    RAISE DEBUG 'Should heal edges % (now %) and % (now %) bound by node %',
      rec.edge1, edge1, rec.edge2, edge2, rec.node_id;

    IF edge1 = edge2 THEN
      -- Nothing to merge here, continue
      CONTINUE;
    END IF;
    ok := false;

    BEGIN
      -- TODO: replace ST_ModEdgeHeal loop with a faster direct deletion and healing
      removedNode := topology.ST_ModEdgeHeal(topo.name, edge1, edge2);
      IF rec.node_id != removedNode THEN
        RAISE EXCEPTION 'Healing of edges % and % was reported '
                        'to remove node % while we expected % instead',
                        edge1, edge2, removedNode, rec.node_id;
      END IF;
      RAISE DEBUG 'Edge % merged into %, dropping node %', edge2, edge1, removedNode;
      ok := 1;
    EXCEPTION WHEN OTHERS
    THEN
      RAISE WARNING 'Edges % and % joined by node % could not be healed: %', edge1, edge2, rec.node_id, SQLERRM;
    END;
    IF ok THEN
      -- edge2 was now renamed to edge1, update map
      edgeMap := jsonb_set(edgeMap, ARRAY[edge2::text], to_jsonb(edge1));
      deletedNodesDeg2 := deletedNodesDeg2 + 1;
    END IF;
  END LOOP; --}
  RAISE NOTICE 'Removed % unneeded nodes of degree 2', deletedNodesDeg2;

  DROP TABLE pg_temp.deleted_edges;
  DROP TABLE pg_temp.merged_faces;

  RETURN deletedEdges + deletedNodes + deletedNodesDeg2;
END;
$BODY$ LANGUAGE 'plpgsql' VOLATILE;



CREATE OR REPLACE FUNCTION topology.postgis_topology_scripts_installed() RETURNS text
	AS $$ SELECT trim('3.3.2'::text || $rev$ 3.3.2 $rev$) AS version $$
	LANGUAGE 'sql' IMMUTABLE;

-- Make sure topology is in database search path --
SELECT topology.AddToSearchPath('topology');

-- Make metadata tables publically visible
-- See https://trac.osgeo.org/postgis/ticket/4575
GRANT SELECT ON topology.topology TO PUBLIC;
GRANT SELECT ON topology.layer TO PUBLIC;
-- See https://trac.osgeo.org/postgis/ticket/5124
GRANT USAGE ON SCHEMA topology TO PUBLIC;

COMMIT;
