[postgis-users] postgis topology

[Ludovic Granjon]

Hi all

I try to import a polygon layer to postgis and build topology with
tolerance.
I try something like that

SELECT CreateTopology('ec_topo2', 27572);
SELECT AddTopoGeometryColumn('ec_topo2', 'public', 'ec', 'topogeom2',
'MULTIPOLYGON');
UPDATE ec SET topogeom2 = toTopoGeom(geom, 'ec_topo', 1, 1.0);

But when I do that, I have

********** Erreur **********

ERREUR: Spatial exception - geometry intersects edge 262
État SQL :P0001
Contexte : fonction PL/pgsql « topogeo_addlinestring », ligne 124 à
affectation
SQL statement "SELECT array_cat(edges, array_agg(x)) FROM ( select
topology.TopoGeo_addLinestring(atopology, rec.geom, tol) as x ) as foo"
fonction PL/pgsql « topogeo_addpolygon », ligne 24 à affectation
fonction PL/pgsql « totopogeom », ligne 112 à FOR sur des lignes de SELECT
fonction PL/pgsql « totopogeom », ligne 94 à affectation

I try to modify the tolerance parameter but it still doesn't work

Have you a solution for that ?

Thanks a lot

Regards

Ludovic

[Christophe Vergon]

Hello, (bonjour)

When you udpdate a topology by using a SET statement if a polygon
intersect an other you will have this error.
If you want to create a topology from polygons try to use the
topogeo_addpolygon function, the intersection between two polygons will
be a new face.

If, as I mean, you use a french survey "cadastre" data set, it's the
best way to do that.



Ludovic Granjon a écrit :

> _______________________________________________
> postgis-users mailing list
> postgi...@lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users


_______________________________________________
postgis-users mailing list
postgi...@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[Ludovic Granjon]

Hello (Bonjour Christophe)

I don't use french cadastre, only specific archaelogical data.

Thanks for this, I try but, I don't realy understand how worked topogeo_addpolygon.
If I do :

CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
SET search_path = topology,public;

I import my layer with QGIS

SELECT topology.CreateTopology('ec_topo', 27572);
SELECT topology.AddTopoGeometryColumn('ec_topo','public','ec','topo_geom','POLYGON');

DO $$DECLARE r record;
BEGIN
  FOR r IN SELECT * FROM ec LOOP
    BEGIN
    SELECT topology.TopoGeo_AddPolygon('ec_topo',st_geometryn(r.geom,1),1);
    EXCEPTION
      WHEN OTHERS THEN
        RAISE WARNING 'Loading of record % failed: %', r.id, SQLERRM;
    END;
  END LOOP;
END$$;

I have :
ATTENTION:  Loading of record 1 failed: la requête n'a pas de destination pour les données résultantes

I don't know where I'm wrong

If you have an idea

Thanks

Regards

Ludovic

PS : (si j'avais sur que c'était toi qui répondrait, j'aurai écrit sur nos forums francophone ;-) )

[Christophe Vergon]

Hello,

Is st_geometryN(r.geomp,1)  a valid geometry ?



Ludovic Granjon a écrit :

[Ludovic Granjon]

Hello

I think, I can do
SELECT st_geometryN(geom,1) FROM ec;
without error

I get

"POLYGON((504935.641400002 2310827.805,504952.193899997 2310825.6349,504964.4124 2310824.033,504964.191799998 2310823.56,504964.037199996 2310822.9832,504963.987400003 2310822.5117,504948.949900001 2310822.7941,504940.407399997 2310823.3589,504932.090300001 2310824.5769,504931.797399998 2310824.7718,504931.360100001 2310824.9757,504930.894100003 2310825.1006,504930.735399999 2310825.1239,504929.7368 2310825.2555,504929.7368 2310825.8138,504929.847599998 2310828.5646,504935.641400002 2310827.805))"

Thanks

Ludovic

[Rémi Cura]

You can use ST_IsValid() to know if it is "valid" for postgres.

Even best, you can try to see what it looks like in qgis.

Cheers,

Rémi-C

2014/1/15 Ludovic Granjon <ludovic...@u-bourgogne.fr>

[Ludovic Granjon]

Hi
thanks
With
SELECT ST_IsValid(st_geometryN(geom,1)) FROM ec;
all is true
no problem in qgis

regards

Ludovic

[Sandro Santilli]

On Tue, Jan 14, 2014 at 04:10:45PM +0100, Ludovic Granjon wrote:

> UPDATE ec SET topogeom2 = toTopoGeom(geom, 'ec_topo', 1, 1.0);

...

> ERREUR: Spatial exception - geometry intersects edge 262

...

> Have you a solution for that ?

Do the update in chunks. That way you'll import what can be and isolate
the offending record (if still offending) and then fight with that one
individually.

These kind of error sometime happen when the tolerance you specify
is bigger than the smallest distance between any two vertices in
your input dataset. An ST_SnapToGrid with given tolerance might also
give you better results, but the chunked loading is always the best
way to do topology building.

--strk;

[Rémi Cura]

Last chance :

this function is a a sloghtly modified topogeom constructor that won't stop when habving an error.

This way you can import everything that works, and then deal with the unworking ...

2014/1/15 Sandro Santilli <st...@keybit.net>

[Rémi Cura]

-- Function: public.rc_totopogeomwithouterrors(geometry, character varying, integer, double precision)

-- DROP FUNCTION public.rc_totopogeomwithouterrors(geometry, character varying, integer, double precision);

CREATE OR REPLACE FUNCTION public.rc_totopogeomwithouterrors(ageom geometry, atopology character varying, alayer integer, atolerance double precision DEFAULT 0)

  RETURNS topogeometry AS

$BODY$

DECLARE

  layer_info RECORD;

  topology_info RECORD;

  rec RECORD;

  tg topology.TopoGeometry;

  elems topology.TopoElementArray = '{{0,0}}';

  sql TEXT;

  typ TEXT;

  tolerance FLOAT8;

BEGIN

raise notice 'ageom : %',ST_AsText(ageom);

  -- 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 tolerance, if 0 was given

  tolerance := COALESCE( NULLIF(atolerance, 0), topology._st_mintolerance(atopology, 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 = 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;

  -- Now that we have a 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

--modified to catch errors:

BEGIN

  FOR rec IN SELECT DISTINCT id(tg), alayer as lyr,-- geom,

    CASE WHEN ST_Dimension(geom) = 0 THEN 1

         WHEN ST_Dimension(geom) = 1 THEN 2

         WHEN ST_Dimension(geom) = 2 THEN 3

    END as type,

    CASE WHEN ST_Dimension(geom) = 0 THEN

           topology.topogeo_addPoint(atopology, geom, tolerance)

         WHEN ST_Dimension(geom) = 1 THEN

           topology.topogeo_addLineString(atopology, geom, tolerance)

         WHEN ST_Dimension(geom) = 2 THEN

           topology.topogeo_addPolygon(atopology, geom, tolerance)

    END as primitive

    FROM (SELECT (ST_Dump(ageom)).geom) as f

    WHERE NOT ST_IsEmpty(geom)

  LOOP

  --raise notice 'coucou, tg : %, geom %',rec.id, ST_AsText(rec.geom);

    -- 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.type

        || ',' || rec.primitive || ')';

    EXECUTE sql;

  END LOOP;

  RETURN tg;

  

EXCEPTION

    WHEN SQLSTATE 'P0001' THEN

        RAISE NOTICE 'this geometry can''t be converted to topogeom : %

doing nothing',ageom;

RETURN NULL;

END;

END

$BODY$

  LANGUAGE plpgsql VOLATILE STRICT

  COST 100;

ALTER FUNCTION public.rc_totopogeomwithouterrors(geometry, character varying, integer, double precision)

  OWNER TO postgres;

2014/1/15 Rémi Cura <remi...@gmail.com>

[Ludovic Granjon]

Thanks all, so

if I try like Sandro said :

UPDATE ec SET geom = ST_SnapToGrid(geom, 1);
UPDATE ec SET topo_geom = toTopoGeom(geom, 'ec_topo', 1, 1.0);

I have this error :
SQL/MM Spatial exception - curve not simple

If i increase tolerance parameter, I have :
ERREUR: SQL/MM Spatial exception - point not on edge

With Rémi's solution :

UPDATE ec SET topo_geom = rc_totopogeomwithouterrors(geom, 'ec_topo', 1, 1.0);

It works with a lot (76) of NOTICE:  this geometry can't be converted to topogeom :
but other geometry worked

so it's better

Now, I have to know why some geometry doesn't work

Maybe I have to explain more what I'm trying to do

My data overlap and there are not clean. Sometimes polygon should have jointly vertex but it's not the case. (I join a image)

I want to snap when the distance is smaller than tolerance and intersect polygon when they overlap.

Do you think that use topology is the best way to do that ?

Thanks a lot

regards

Ludovic

[Ludovic Granjon]

The image