[postgis-users] ST_CreateTopoGeo

[paul.malm at lfv.se]

Hi again,
I can see in the edge data of topo2 that the main polygon of Egypt intersects with an island next to the main polygon. If I delete the island in the original layer before running the SQL sequence the main polygon of Egypt is included in the result layer. Is there any solution of this problem?

Kind regards,
Paul

-----Ursprungligt meddelande-----
Från: postgis-users [mailto:postgis-users-bounces at lists.osgeo.org] För Malm, Paul (Operations AIM)
Skickat: den 28 februari 2020 11:36
Till: postgis-users at lists.osgeo.org; strk at kbt.io
Ämne: Re: [postgis-users] ST_CreateTopoGeo

Thanks, Sandro!
It worked with your LOOP! Great!
The problem I can see on the result (I'm simplifying country polygons) is that the large polygon of Finland, Germany, Italy and Egypt is missing.
I've tried to change the different tolerances (Simplify tolerance and your loop snap-tolerance). If I use simplifying tolerance 2500 m the missing polygons are there, but not when using 5000 m. The countries included in the  original layer is from Sweden, Norway and Finland down to Jordan, Egypt, Libya and Tunisia.
Ref sys =WGS 84_UTM-33N


If you have the time and want to look at the SQL sequence, here it is;
Thank you,
Paul

CREATE TABLE "simple_countries2nd" as (
       SELECT "NAME_ZH", "OGR_STYLE", fid, (st_dump(the_geom)).*
       FROM "countries2nd");

CREATE INDEX "simple_countries2nd_geom_gist" ON "simple_countries2nd" USING gist(geom);

ALTER TABLE "simple_countries2nd" ADD COLUMN simple_geom geometry(POLYGON, 32633);

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

DO $$
DECLARE
  rec RECORD;
BEGIN
  FOR rec in SELECT "NAME_ZH", "OGR_STYLE", "fid", (ST_Dump(the_geom)).geom FROM "countries2nd"
  LOOP
    BEGIN
      IF GeometryType(rec.geom) = 'POLYGON' THEN
         PERFORM topology.TopoGeo_AddPolygon('topo1', rec.geom, 400);
      ELSIF GeometryType(rec.geom) = 'LINESTRING' THEN
        PERFORM topology.TopoGeo_AddLinestring('topo1', rec.geom, 400);
      ELSIF GeometryType(rec.geom) = 'POINT' THEN
        PERFORM topology.TopoGeo_AddPoint('topo1', rec.geom, 400);
      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_countries2nd" 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;
	

-----Ursprungligt meddelande-----
Från: postgis-users [mailto:postgis-users-bounces at lists.osgeo.org] För Sandro Santilli
Skickat: den 25 februari 2020 15:38
Till: PostGIS Users Discussion
Ämne: Re: [postgis-users] ST_CreateTopoGeo

On Tue, Feb 25, 2020 at 01:22:39PM +0000, paul.malm at lfv.se wrote:
> Hi,
> Where do you mean I can play with the tolerance?

In TopoGeo_addLinestring
https://postgis.net/docs/manual-3.1/TopoGeo_AddLineString.html

> This is what I have done before the ST_createTopoGeo
> SELECT topology.CreateTopology('topo1', 4326)";
> SELECT topology.ST_CreateTopoGeo('topo1', ST_Collect(geom)) from "countries_first";
> Btw, I'm intending to simplify later on in my SQL command list.

You could do something like this:

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

You can tweak the above to do something different rather than raising
a WARNING (for example store ID of offending geoms in a table).
Then you can look at the offending geometries in isolation, possibly
tweaking the "tol" variable.

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