[postgis-users] ST_CreateTopoGeo

[paul.malm at lfv.se]

In both topo1- and topo2.face shows a box where Egypt is located, as for all other countries (type of bounding box?).
In both topo1- and topo2.edge_data shows all countries as well as Egypt. Egypt is however intersecting an island in topo2. The island is formed like a banana in topo1 (not intersecting the Egypt polygon) and like a triangle in topo2, due to the simplification (5000 m). Therefore a headland in the Egypt-polygon is intersecting the island in topo2.edge_data.
I think that's why I can't create a polygon geometry for Egypt. It is the same for Finland (also has an island that becomes intersecting the main polygon of Finland in topo2).

Here is the whole SQL sequence:

DROP TABLE  IF EXISTS "simple_countries";
DELETE FROM topology.topology WHERE name = 'topo1';
DELETE FROM topology.topology WHERE name = 'topo2';
DROP SCHEMA IF EXISTS topo2 CASCADE;
DROP SCHEMA IF EXISTS topo1 CASCADE;
CREATE TABLE "simple_countries" AS (
       SELECT "featurecla", "scalerank",   "fid", (st_dump(the_geom)).*
       FROM "countries");
CREATE INDEX "simple_countries_geom_gist" on "simple_countries" using gist(geom);
SELECT topology.CreateTopology('topo1', 32633);

DO $$
DECLARE
  rec RECORD;
  tol FLOAT8;
BEGIN
  tol := 10;
  FOR rec in SELECT featurecla, scalerank,  fid, (ST_Dump(the_geom)).geom FROM "countries"
  LOOP
    BEGIN
      IF GeometryType(rec.geom) = 'POLYGON' THEN
         PERFORM topology.TopoGeo_AddPolygon('topo1', rec.geom, 10);
      ELSIF GeometryType(rec.geom) = 'LINESTRING' THEN
        PERFORM topology.TopoGeo_AddLinestring('topo1', rec.geom, 10);
      ELSIF GeometryType(rec.geom) = 'POINT' THEN
        PERFORM topology.TopoGeo_AddPoint('topo1', rec.geom, 10);
      END IF;
    EXCEPTION WHEN OTHERS THEN
      RAISE WARNING 'For geometry % we got exception % (%)', rec.id, SQLERRM, SQLSTATE;
    END;
  END LOOP;
END;
$$ LANGUAGE 'plpgsql';

SELECT topology.CreateTopology('topo2', 32633);

SELECT topology.ST_CreateTopoGeo('topo2', the_geom)
FROM (
       SELECT ST_Collect(st_simplifyPreserveTopology(geom, 5000)) AS the_geom
       FROM topo1.edge_data
) AS foo;

WITH simple_face AS (
       SELECT topology.st_getFaceGeometry('topo2', face_id) AS the_geom
       FROM topo2.face
       WHERE face_id > 0
) 
UPDATE "simple_countries" d set geom = sf.the_geom
FROM simple_face sf
WHERE st_intersects(d.geom, sf.the_geom)
AND st_area(st_intersection(sf.the_geom, d.geom))/st_area(sf.the_geom) > 0.5;

ALTER TABLE "simple_countries" RENAME COLUMN geom TO the_geom;

/Paul
-----Ursprungligt meddelande-----
Från: postgis-users [mailto:postgis-users-bounces at lists.osgeo.org] För Sandro Santilli
Skickat: den 3 mars 2020 18:23
Till: PostGIS Users Discussion
Ämne: Re: [postgis-users] ST_CreateTopoGeo

On Tue, Mar 03, 2020 at 02:30:11PM +0000, paul.malm at lfv.se wrote:
> Hi,
> I'm not sure of what you mean, Strk.
> I create a topology (topo1) with all geometries from a polygon layer.
> I can see the island and the Egypt polygon in the edge_data, they are not intersecting.
> 
> Then I create a new topology (topo2) with:
> select topology.ST_CreateTopoGeo('topo2', the_geom) from ( 
>        select ST_Collect(st_simplifyPreserveTopology(geom, 5000)) as the_geom
>        from topo1.edge_data) as foo;

Is Egypt still there when you look at topology faces ?

> Here in the topo2.edge.data I can see that the island and the polygon of Egypt is intersecting due to the simplifying.
> And when I create a new polygon layer from this edge_data the Egypt polygon is left out.
> 
> This is done with:
> with simple_face as ( select topology.st_getFaceGeometry('topo2', face_id) as the_geom 
>        from topo2.face 
>        where face_id > 0 ) update newlayer d set geom = sf.the_geom 
> from simple_face sf 
> where st_intersects(d.geom, sf.the_geom)
> st_area(st_intersection(sf.the_geom, d.geom))/st_area(sf.the_geom) > 0.5";

I suspect you're just filtering out too many faces, due
to small intersection area..

--strk;
_______________________________________________
postgis-users mailing list
postgis-users at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users