[postgis-users] Topology: cannot delete slivers (or gaps)

[Guillaume Drolet]

Hi!

I'm am turning to you all hoping you'll have some tricks to help with
this problem:

I created a topology for 2290 polygons and populated a topogeom column
using these commands:

SELECT topology.CreateTopology('de_20k_topo', find_srid('syshiera',
'de_20k', 'geom_32198'), 1e-4);
SELECT topology.AddTopoGeometryColumn('de_20k_topo', 'syshiera',
'de_20k', 'topogeom', 'MULTIPOLYGON');

DO $$DECLARE r record;
BEGIN
FOR r IN SELECT * FROM syshiera.de_20k LOOP
BEGIN
UPDATE syshiera.de_20k SET topogeom = toTopoGeom(geom_32198,
'de_20k_topo', 1, 1e-4)
WHERE gid = r.gid;
EXCEPTION
WHEN OTHERS THEN
RAISE WARNING 'Loading of record % failed: %', r.gid, SQLERRM;
END;
END LOOP;
END$$;

-- For testing:
SELECT topology.CopyTopology('de_20k_topo', 'de_20k_topo_tests');


I need to fix 1) extra nodes along edges that are not at an
intersection and 2) sliver polygons.

For cases in 1, I was successful using 'ST_NewEdgeHeal' and a plPgsql script.

For cases in 2, I can remove some slivers with 'ST_RemEdgeModFace' but
some of them I just can't delete:

I get either:

ERROR: TopoGeom 1123 in layer 1 (de_20k_topo_tests.LAYER1.) cannot be
represented healing faces 1514 and 1502

Or, when trying with remove a node with 'topology.ST_NewEdgeHeal':

ERROR: SQL/MM Spatial exception - other edges connected (3912)

In this particular example, two edges are connected two the same two
nodes, similar to what was discussed here:

http://postgis.17.x6.nabble.com/PostGIS-1955-Exception-when-2-edges-passed-to-ST-ModEdgeHeal-that-are-attached-to-the-same-2-nodes-td4999383.html

But according to that discussion thread, this was fixed in an earlier
version (I'm using topology version 2.1.3).

Maybe my problem is different but I'll need your help to get it
sorted. Here's an example showing the small face (1502) I want remove:

https://dl.dropboxusercontent.com/u/5196336/example1.bmp

Thanks for your help!
_______________________________________________
postgis-users mailing list
postgi...@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[Sandro Santilli]

On Fri, Oct 31, 2014 at 10:17:48AM -0400, Guillaume Drolet wrote:

> I created a topology for 2290 polygons and populated a topogeom column

[...]

> I need to fix 1) extra nodes along edges that are not at an
> intersection and 2) sliver polygons.

[...]

> For cases in 2, I can remove some slivers with 'ST_RemEdgeModFace' but
> some of them I just can't delete:
>
> I get either:
>
> ERROR: TopoGeom 1123 in layer 1 (de_20k_topo_tests.LAYER1.) cannot be
> represented healing faces 1514 and 1502

This error means that a TopoGeometry object exists (with id=1123)
that's defined as containing only one of the faces separated by the
edge you're trying to remove.

The fix here would be to find faces that take part in the composition
of multiple TopoGeometry objects and those NOT taking part in the
composition of any TopoGeometry object. I've seen GRASS GUIs referring
to those 2 classes as "polygon0" and "polygon2", enough that I've been
thinking it could be useful to have them shown as layers in the QGIS
viewer for PostGIS Topology. You should be able to extract those faces
with a query to the relation and the layers tables.

So basically you first have to fix the TopoGeometry objects and then
can safely drop the no-more-needed edges.

> Or, when trying with remove a node with 'topology.ST_NewEdgeHeal':
>
> ERROR: SQL/MM Spatial exception - other edges connected (3912)
>
> In this particular example, two edges are connected two the same two
> nodes, similar to what was discussed here:
>
> http://postgis.17.x6.nabble.com/PostGIS-1955-Exception-when-2-edges-passed-to-ST-ModEdgeHeal-that-are-attached-to-the-same-2-nodes-td4999383.html
>
> But according to that discussion thread, this was fixed in an earlier
> version (I'm using topology version 2.1.3).
>
> Maybe my problem is different but I'll need your help to get it
> sorted. Here's an example showing the small face (1502) I want remove:
>
> https://dl.dropboxusercontent.com/u/5196336/example1.bmp

Which edge ids did you pass to ST_NewEdgeHeal ?
Is the topology in a good state according to topology.ValidateTopology ?

--strk;

() Free GIS & Flash consultant/developer
/\ http://strk.keybit.net/services.html

[Guillaume Drolet]

Hi Sandro,

Thanks for your help!

topology.ValidateTopology gives this:

error id1 id2
character varying integer integer
------------------ ----------- ----------------
face has no rings 1831

If I calculate the area of face 1831
(st_area(st_getfacegeometry('de_20k_topo', face_id))), I get a NULL value.

I also find I have some 650 of these slivers with an area below 75 m2 while
most of my correct polygons have an area above 300 000 000 m2.

Re the edges passed to ST_NewEdgeHeal, in the example I sent it was edges
3908 and 3909.

If you think about something else that could be problematic, let me know.
Meanwhile I'll dig into the approach you suggested in your reply (i.e.
fixing TopoGeometry objects).

Thanks!

PS. Forgot to mention that in the first topology I had created, I used a
tolerance of 2.0 meters instead of 1e-4 meters. While that yielded much less
slivers, I was getting way more of those nodes not located at an
intersection. What is best?





--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007257.html
Sent from the PostGIS - User mailing list archive at Nabble.com.

[Sandro Santilli]

On Fri, Oct 31, 2014 at 12:19:31PM -0700, Guillaume Drolet wrote:
> Hi Sandro,
>
> Thanks for your help!
>
> topology.ValidateTopology gives this:
>
> error id1 id2
> character varying integer integer
> ------------------ ----------- ----------------
> face has no rings 1831
>
> If I calculate the area of face 1831
> (st_area(st_getfacegeometry('de_20k_topo', face_id))), I get a NULL value.

I suggest you take a look at edges having face 1831 on either
the right or left edge. Do any edge exist ?

> I also find I have some 650 of these slivers with an area below 75 m2 while
> most of my correct polygons have an area above 300 000 000 m2.

The correct procedure to drop those slivers would be to assign
those small areas to the (closest|biggest|leftmost|youpickit)
TopoGeometry object, to one and only one, and then, if you really
need it, drop the edge internal to the TopoGeometry.

I made a QGIS plugin which may help you do some of these things
manually, which is useful in some cases, mainly to understand the
issues.

> Re the edges passed to ST_NewEdgeHeal, in the example I sent it was edges
> 3908 and 3909.

So the error was right about there being another edge attached to the node:

ERROR: SQL/MM Spatial exception - other edges connected (3912)

I'm guessing 3912 and 3908 are the two parallel edges in the picture
you attached: https://dl.dropboxusercontent.com/u/5196336/example1.bmp

Again, QGIS may help you getting more info out of the visual, enable
the DBManager core plugin, select your topology schema and select
Schema->TopoViewer from the menu.

> PS. Forgot to mention that in the first topology I had created, I used a
> tolerance of 2.0 meters instead of 1e-4 meters. While that yielded much less
> slivers, I was getting way more of those nodes not located at an
> intersection. What is best?

It's probably easier to clean up nodes for an areal topology than to
clean up faces. With my QGIS plugin (PostGIS Topology Editor) you can
select all nodes and click on a button which will drop all the ones
that are safe to drop.

--strk;

() Free GIS & Flash consultant/developer
/\ http://strk.keybit.net/services.html

[Guillaume Drolet]

Sandro Santilli wrote

>> I suggest you take a look at edges having face 1831 on either
>> the right or left edge. Do any edge exist ?
>

> edge 5187 has face 1831 as its right face

>
>> The correct procedure to drop those slivers would be to assign
>> those small areas to the (closest|biggest|leftmost|youpickit)
>> TopoGeometry object, to one and only one, and then, if you really
>> need it, drop the edge internal to the TopoGeometry.
>
>> I made a QGIS plugin which may help you do some of these things
>> manually, which is useful in some cases, mainly to understand the
>> issues.
>

> The plugin you made is a great tool (thanks). I will first try to get the
> number of slivers down to
> a manageable number, which I had with setting the tolerance at 2.0 m, and
> then apply the
> TopoGeometry approach above.

>
>> So the error was right about there being another edge attached to the
>> node:
>
>> ERROR: SQL/MM Spatial exception - other edges connected (3912)
>
>> I'm guessing 3912 and 3908 are the two parallel edges in the picture
>> you attached: https://dl.dropboxusercontent.com/u/5196336/example1.bmp
>

> You're right, they are the two parallel edges. And it's the same for most
> of
> the some 650 slivers I have.

>
>> Again, QGIS may help you getting more info out of the visual, enable
>> the DBManager core plugin, select your topology schema and select
>> Schema->TopoViewer from the menu.
>

> I didn't know about the TopoViewer tool. I was importing the topology
> tables (node, edge, face)
> like any other PG geometry tables. Will explore with this new tool.

>
>
>> It's probably easier to clean up nodes for an areal topology than to
>> clean up faces. With my QGIS plugin (PostGIS Topology Editor) you can
>> select all nodes and click on a button which will drop all the ones
>> that are safe to drop.
>

> You're right! I'll create my topology again but with the 2.0 m tolerance,
> then run
> my function for getting rid of nodes that are not at intersections and
> finally,
> I will use the approach you suggested of working with the TopoGeometry
> objects and the
> relation table to eliminate the remaining slivers.
>
> Will keep you posted on how it goes. Again, I really appreciate your help.
> Thanks.
>
> Guillaume.





--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007266.html

Sent from the PostGIS - User mailing list archive at Nabble.com.

[Rémi Cura]

Just for completeness sake,

if your input geom is really ambiguous with the precision you choose,

it may be faster in both human and machine time to use grass 7 to import and clean the topology,
then export from grass 7 to postgis topology with the built-in grass function.

I ressort sometime to this trick because grass cleaning functionality are very advanced.

Cheers,
Rémi-C

[John Abraham]

Just FYI in case you’re interested I found it much easier to use the triangle algorithm here https://github.com/tudelft3d/pprepair on an exported shape file 

It’s called Planar Partition Repair.  It got rid of all the slivers and overlaps with ease.

Here is the paper http://www.gdmc.nl/ken/files/12_pfg.pdf and thanks to Martijn Meijers and team.

Then I rebuilt the topology layer in PostGIS from the fixed up Simple Feature Class multi polygon layer.

--

John Abraham

j...@hbaspecto.com

403-232-1060

[Guillaume Drolet]

Hi Sandro and others,

First, thanks Jon for the hint on pprepair: I ran pprepair on my shapefile
and I'm now in the process of
importing it into PostGIs for topology creation. Then I'll be able to tell
if pprepair really fixed everything.

But because I'm stubborn and never give up, I still want to succeed on
fixing TopoGeometry objects in PostGIS (bear with me, what follows is a bit
long):

Sandro, you said:

> The fix here would be to find faces that take part in the composition
> of multiple TopoGeometry objects and those NOT taking part in the
> composition of any TopoGeometry object. I've seen GRASS GUIs referring
> to those 2 classes as "polygon0" and "polygon2", enough that I've been
> thinking it could be useful to have them shown as layers in the QGIS
> viewer for PostGIS Topology. You should be able to extract those faces
> with a query to the relation and the layers tables.

> So basically you first have to fix the TopoGeometry objects and then
> can safely drop the no-more-needed edges.

First, I looked at faces NOT taking part in the composition of any
TopoGeometry object:

/SELECT face_id FROM de_20k_topo.face WHERE face_id NOT IN (SELECT
element_id

FROM de_20k_topo.relation);/
face_id
0
2418
2650
2663
2766
3077
3118
3146
3267
3270

None of these 10 faces are involved with my three remaining slivers... And
when I look at them, they ARE related to edges, as a next_face or a
right_face. Is there any useful information in that for my purpose?

Second, I tried to work along the lines suggested by Sandro, using a new
topology (tolerance 2.0 m). As in one of my previous attempts, I
successfully fixed 20 of the 23 slivers using ST_RemEdgeModFace. For the
remaining three slivers, however, I tried modyfying TopoGeometry objects but
with no success. I must admit I'm a bit lost and not too sure if I really
understand what TopoGeometry objects are:

For example, let's look at problematic face_id: 2227
(https://dl.dropboxusercontent.com/u/5196336/example_2227.bmp). In this
case, the light grey line is where a single edge should be. Blue edges 6381
and 6341 should not be there. The light grey line should start from the node
at the bottom right and connect at its other end somewhere in the middle of
edge 6340:

I thought, let's first remove edge 6346, which should not exist in the first
place, and then split edge 6340 with a node and connect it to node 4753 on
the lower right corner of the image:

SELECT ST_RemEdgeModFace('de_20k_topo', 6346);

As in my previous post, this yields:

ERROR: TopoGeom 1546 in layer 1 (syshiera.de_20k.topogeom) cannot be
represented healing faces 2225 and 2227
CONTEXT: instruction SQL « SELECT topology._ST_RemEdgeCheck(toponame,
topoid, e1id, e1rec.left_face, e1rec.right_face) »
function PL/pgsql st_remedgemodface(character varying,integer), line 43 at
PERFORM

So I looked at TopoGeom 1546:

/SELECT * FROM de_20k_topo.relation WHERE topogeo_id = 1546;/

topogeo_id;layer_id;element_id;element_type
1546;1;2225;3
1546;1;2233;3

Edge 2227 is missing from TopoGeom 1546: Where is it?

/SELECT * FROM de_20k_topo.relation WHERE element_id = 2227;/

topogeo_id;layer_id;element_id;element_type
1551;1;2227;3

There it is, in TopoGeom 1551. What else is there:

/SELECT * FROM de_20k_topo.relation WHERE topogeo_id = 1551;/

topogeo_id;layer_id;element_id;element_type
1551;1;2225;3
1551;1;2234;3
1551;1;2227;3

So edges 2225 and 2227 are associated with TopoGeom 1551. Will change that:

/
UPDATE de_20k_topo.relation SET topogeo_id = 1546
WHERE topogeo_id = 1551 AND element_id = 2227;/

Query returned successfully: one row affected, 481 ms execution time.

And then again:

/SELECT ST_RemEdgeModFace('de_20k_topo', 6346); /

ERROR: TopoGeom 1551 in layer 1 (syshiera.de_20k.topogeom) cannot be
represented healing faces 2225 and 2227
CONTEXT: instruction SQL « SELECT topology._ST_RemEdgeCheck(toponame,
topoid, e1id, e1rec.left_face, e1rec.right_face) »
function PL/pgsql st_remedgemodface(character varying,integer), line 43 at
PERFORM

Obviously, as mentionned above, I don't quite understand how TopoGeometry
objects work and the link between TopoGeoms, geometries and primitives. I
need a bit of education here.

Thanks a lot for your patience and your help.

PS. Remi, thanks for suggesting GRASS but I've never used it and the
learning curve seems pretty steep, no?

--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007273.html

Sent from the PostGIS - User mailing list archive at Nabble.com.

[Sandro Santilli]

On Tue, Nov 04, 2014 at 03:08:15PM -0800, Guillaume Drolet wrote:

> But because I'm stubborn and never give up, I still want to succeed on
> fixing TopoGeometry objects in PostGIS (bear with me, what follows is a bit
> long):

I like you ! :)

> First, I looked at faces NOT taking part in the composition of any
> TopoGeometry object:
>
> /SELECT face_id FROM de_20k_topo.face WHERE face_id NOT IN (SELECT
> element_id
>
> FROM de_20k_topo.relation);/
> face_id
> 0
> 2418
> 2650
> 2663
> 2766
> 3077
> 3118
> 3146
> 3267
> 3270
>
> None of these 10 faces are involved with my three remaining slivers... And
> when I look at them, they ARE related to edges, as a next_face or a
> right_face. Is there any useful information in that for my purpose?

Except face 0, which is the "universe" face, all the others
should be "gaps" in your topology. You can zoom to the face's
bounding box ("mbr" field) to take a closer look (do you have qgis
set up by now?).

What you could do programmatically is for each face to fetch the list
of edges bounding it and pick the face on the other side as a candidate;
then for each candidate face you find TopoGeometry objects containing
it in their definition and pick one to assign that orphaned face
(you could use smaller TopoGeometry, or larger, or whatever).

It might be better first to deal with the faces that are associated
with multiple TopoGeometry, just to avoid assigning an orphaned face
to a TopoGeometry based on it's containing a face which will later
be removed from the definition because duplicated.

> Second, I tried to work along the lines suggested by Sandro, using a new
> topology (tolerance 2.0 m). As in one of my previous attempts, I
> successfully fixed 20 of the 23 slivers using ST_RemEdgeModFace. For the
> remaining three slivers, however, I tried modyfying TopoGeometry objects but
> with no success. I must admit I'm a bit lost and not too sure if I really
> understand what TopoGeometry objects are:

Eh, that part is a still a bit weak in documantation, maybe.
A TopoGeometry is a Geometry defined by its components.

In your case (areal TopoGeometry objects) it is a (multi)Polygon
defined by a set of faces. Note that QGIS is also able to
show you the TopoGeometry objects, and lets you edit them "spatially"
(less robust than editing them by adding/removing components but with
some care easier to use).

> For example, let's look at problematic face_id: 2227
> (https://dl.dropboxusercontent.com/u/5196336/example_2227.bmp). In this
> case, the light grey line is where a single edge should be. Blue edges 6381
> and 6341 should not be there. The light grey line should start from the node
> at the bottom right and connect at its other end somewhere in the middle of
> edge 6340:

Uhm, edge identifiers are missing from that dump.
I take it 6381 is the lower blue and 6341 the upper blue.

> I thought, let's first remove edge 6346, which should not exist in the first
> place, and then split edge 6340 with a node and connect it to node 4753 on
> the lower right corner of the image:
>
> SELECT ST_RemEdgeModFace('de_20k_topo', 6346);
>
> As in my previous post, this yields:
>
> ERROR: TopoGeom 1546 in layer 1 (syshiera.de_20k.topogeom) cannot be
> represented healing faces 2225 and 2227

I take it that edge 6346 is the one on the right of the yellow face.
The message is telling you that a TopoGeometry object exists (id=1546)
that is composed only by one of the two faces 2225 and 2227.

I suggest you load the TopoGeometry layer with qgis and see for yourself.
Unfortunately (that'd be useful to add in qgis) it won't be easy to
find TopoGeometry 1546 as extracting the id of a TopoGeometry requires
running an expression on the database: id(topogeom).

You could create a view selecting just that topoGeometry, or a subset
of them, or you could just load them all, hoping it won't be too expensive
(another spot where qgis support for postgis topology could be improved).

Chances are it would be enough to assign face 2227 to the TopoGeometry
with id 1546 (likely the one not-assigned). But as face 2227 was not
in your initial set (unassigned faces) I guess there's another TopoGeometry
which does have 2227 in his definition, but not the one on the right.
In that case you have to pick which of the two contending TopoGeometry
objects should get face 2227 assigned, before you drop it. Basically, which
TopoGeometry will get assigned the space occupied by that face.

Given QGIS editing support for TopoGeometry objects you should be able
to see all of this by loading the layer and setting 50% transparency.
Enabling topological editing and editing those two TopoGeometry objects
should help you with manual resolution of this issue.

> So I looked at TopoGeom 1546:
>
> /SELECT * FROM de_20k_topo.relation WHERE topogeo_id = 1546;/
>
> topogeo_id;layer_id;element_id;element_type
> 1546;1;2225;3
> 1546;1;2233;3
>
> Edge 2227 is missing from TopoGeom 1546: Where is it?

You meant _face_ 2227.
Right, it's not assigned to TopoGeometry 1546,
but to another one.

> /SELECT * FROM de_20k_topo.relation WHERE element_id = 2227;/
>
> topogeo_id;layer_id;element_id;element_type
> 1551;1;2227;3
>
> There it is, in TopoGeom 1551. What else is there:
>
> /SELECT * FROM de_20k_topo.relation WHERE topogeo_id = 1551;/
>
> topogeo_id;layer_id;element_id;element_type
> 1551;1;2225;3
> 1551;1;2234;3
> 1551;1;2227;3
>
> So edges 2225 and 2227 are associated with TopoGeom 1551. Will change that:

Again, you mean faces (last number is element_type, 3==face).
So TopoGeometry objects 1546 and 1551 overlap, and their intersection
is the union of faces 2225 and 2227.

> UPDATE de_20k_topo.relation SET topogeo_id = 1546
> WHERE topogeo_id = 1551 AND element_id = 2227;/

There's already a record associating face 2227 with TopoGeometry 1546
so you don't need another record, just drop the one associating it
with TopoGeometry 1551:

DELETE FROM de_20k_topo.relation
WHERE topogeo_id = 1551 AND element_id = 2227;

> And then again:
>
> /SELECT ST_RemEdgeModFace('de_20k_topo', 6346); /
>
> ERROR: TopoGeom 1551 in layer 1 (syshiera.de_20k.topogeom) cannot be
> represented healing faces 2225 and 2227

Right, you only resolved the overlap of face 2227, but there's still
the overlap of face 2225 :)

DELETE FROM de_20k_topo.relation
WHERE topogeo_id = 1551 AND element_id = 2225;

> Obviously, as mentionned above, I don't quite understand how TopoGeometry
> objects work and the link between TopoGeoms, geometries and primitives. I
> need a bit of education here.

It looks like you're on the right track so far.

> Thanks a lot for your patience and your help.

Thank you for your stubbornnes in using PostGIS Topology !

--strk;

() Free GIS & Flash consultant/developer
/\ http://strk.keybit.net/services.html

[Guillaume Drolet]

Sandro Santilli wrote

>> Except face 0, which is the "universe" face, all the others
>> should be "gaps" in your topology. You can zoom to the face's
>> bounding box ("mbr" field) to take a closer look (do you have qgis
>> set up by now?).
>

> Ok. I had already looked at some of them but as face-derived geometries
> (using ST_GetFaceGeometry) and not as real faces. As geometries, to me
> they look like polygons and I didn't realize they were in fact holes. See
> this example in QGIS for face_id 2418 turned into a polygon:
> https://dl.dropboxusercontent.com/u/5196336/geomFromFace_id_2418.bmp. Now,
> look at its mbr representation:
> https://dl.dropboxusercontent.com/u/5196336/mbrFace_id_2418.bmp. How do
> you tell it's a hole and not just a regular face?

>
>> What you could do programmatically is for each face to fetch the list
>> of edges bounding it and pick the face on the other side as a candidate;
>> then for each candidate face you find TopoGeometry objects containing
>> it in their definition and pick one to assign that orphaned face
>> (you could use smaller TopoGeometry, or larger, or whatever).
>

> Let's try this or face_id 2418 (which is, btw, a meaningful unit and not
> an error...):
/
> SELECT f.face_id, e.edge_id, e.left_face, e.right_face
> FROM de_20k_topo.face f JOIN
> de_20k_topo.edge_data e ON (f.face_id = e.left_face OR
> f.face_id = e.right_face)
> WHERE face_id = 2418;
/
>
*
> face_id, edge_id, left_face, right_face
*
> 2418;2663;2418;950
> 2418;6986;2418;340
> 2418;2672;2418;951
> 2418;2665;2418;344
> 2418;6989;2418;2417
>
> I pick the first edge (2663) and take it's right_face, 950, as a
> candidate:
/
> SELECT topogeo_id
> FROM de_20k_topo.relation
> WHERE element_id = 950;
/
>
*
> topogeo_id
*
> 814
>
> Only TopoGeom 814 has face_id 950 as one of its elements. Will assign
> face_id 2418
> to this TopoGeom:
/
> INSERT INTO de_20k_topo.relation
> VALUES (814, 1, 2418, 3);
/
>
> Then what? Is that it? Does that mean that that big hole 2418 is now
> unioned with face_id 950 (in which case it isn't what I want as it's a
> meaningful ecological unit)?

>
>
>> It might be better first to deal with the faces that are associated
>> with multiple TopoGeometry, just to avoid assigning an orphaned face
>> to a TopoGeometry based on it's containing a face which will later
>> be removed from the definition because duplicated.
>

> Good point.

>
>> Eh, that part is a still a bit weak in documantation, maybe.
>> A TopoGeometry is a Geometry defined by its components.
>
>> In your case (areal TopoGeometry objects) it is a (multi)Polygon
>> defined by a set of faces. Note that QGIS is also able to
>> show you the TopoGeometry objects, and lets you edit them "spatially"
>> (less robust than editing them by adding/removing components but with
>> some care easier to use).
>

>>> In this case, the light grey line is where a single edge should be. Blue
>>> edges 6381
>>> and 6341 should not be there. The light grey line should start from the
>>> node
>>> at the bottom right and connect at its other end somewhere in the middle
>>> of
>>> edge 6340:
>
>> Uhm, edge identifiers are missing from that dump.
>> I take it 6381 is the lower blue and 6341 the upper blue.
>

> You're correct here.

>
>> I suggest you load the TopoGeometry layer with qgis and see for yourself.
>

> When you say load the TopoGeometry layer, I assume you mean the topogeom
> column I populated with:
/

> UPDATE syshiera.de_20k SET topogeom = toTopoGeom(geom_32198,

> 'de_20k_topo', 1, 2.0)
/
>
> Right?
>
> When I load it in QGIS using the PostGIS table loader, I can visualize it
> and see there are holes. E.g.
> https://dl.dropboxusercontent.com/u/5196336/holeInTopogeometry.bmp

>
>> Unfortunately (that'd be useful to add in qgis) it won't be easy to
>> find TopoGeometry 1546 as extracting the id of a TopoGeometry requires
>> running an expression on the database: id(topogeom).
>

> I can either
/
> SELECT *
> FROM syshiera.de_20k
> WHERE id(topogeom) = 1546;
/
>
> or, as I mentioned before, I can superimpose my 'faces' layer, which I
> derived from de_20k_topo.face
> using ST_GetFaceGeometry, to 'topogeom' to get the topogeo_id. No?

>
>> You could create a view selecting just that topoGeometry, or a subset
>> of them, or you could just load them all, hoping it won't be too
>> expensive
>> (another spot where qgis support for postgis topology could be improved).
>

> I agree, just loading the TopoViewer freezes QGIS for a while and I don't
> lack ressources on my machine.

>
>> Chances are it would be enough to assign face 2227 to the TopoGeometry
>> with id 1546 (likely the one not-assigned). But as face 2227 was not
>> in your initial set (unassigned faces) I guess there's another
>> TopoGeometry
>> which does have 2227 in his definition, but not the one on the right.
>> In that case you have to pick which of the two contending TopoGeometry
>> objects should get face 2227 assigned, before you drop it. Basically,
>> which
>> TopoGeometry will get assigned the space occupied by that face.
>

> TopoGeometry 1551 has 2227. More explanation about this in my previous
> post and below

>
>> Given QGIS editing support for TopoGeometry objects you should be able
>> to see all of this by loading the layer and setting 50% transparency.
>> Enabling topological editing and editing those two TopoGeometry objects
>> should help you with manual resolution of this issue.
>

> Will try this, thanks

>
>>> Edge 2227 is missing from TopoGeom 1546: Where is it?
>
>> You meant _face_ 2227.
>> Right, it's not assigned to TopoGeometry 1546,
>> but to another one.
>

> That's what I meant yes.

>
>>> So edges 2225 and 2227 are associated with TopoGeom 1551. Will change
>>> that:
>
>> Again, you mean faces (last number is element_type, 3==face).
>> So TopoGeometry objects 1546 and 1551 overlap, and their intersection
>> is the union of faces 2225 and 2227.
>
>>> UPDATE de_20k_topo.relation SET topogeo_id = 1546
>>> WHERE topogeo_id = 1551 AND element_id = 2227;
>

>> There's already a record associating face 2227 with TopoGeometry 1546
>> so you don't need another record, just drop the one associating it
>> with TopoGeometry 1551:
>

> Ok, now I see where you're going..

>
>> DELETE FROM de_20k_topo.relation
>> WHERE topogeo_id = 1551 AND element_id = 2227;
>
>>> And then again:
>>>
>>> /SELECT ST_RemEdgeModFace('de_20k_topo', 6346); /
>>>
>>> ERROR: TopoGeom 1551 in layer 1 (syshiera.de_20k.topogeom) cannot be
>>> represented healing faces 2225 and 2227
>
>> Right, you only resolved the overlap of face 2227, but there's still
>> the overlap of face 2225 :)
>
>> DELETE FROM de_20k_topo.relation
>> WHERE topogeo_id = 1551 AND element_id = 2225;
>

> Now it's clearer in my head what I need to do, thanks to you!

>
>> It looks like you're on the right track so far.
>

> Will keep trying and I will win. Again, thanks so much for your help.

--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007278.html

Sent from the PostGIS - User mailing list archive at Nabble.com.

[Guillaume Drolet]

Hi John,

I tried pprepair and imported the fixed shapefile into PostGIS. Although
pprepair said all polygons were fixed, there are still 22 slivers in my
topology. Maybe due to the value I used for tolerance? Finding the right
tolerance for building a topology isn't an easy task...



--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007279.html

[Sandro Santilli]

On Wed, Nov 05, 2014 at 08:32:01AM -0800, Guillaume Drolet wrote:
> Sandro Santilli wrote:
>>
>> Except face 0, which is the "universe" face, all the others
>> should be "gaps" in your topology. You can zoom to the face's
>> bounding box ("mbr" field) to take a closer look (do you have qgis
>> set up by now?).
>
> Ok. I had already looked at some of them but as face-derived geometries
> (using ST_GetFaceGeometry) and not as real faces.

I don't understand the distinction you're making.
A face-derived geometry IS the face's geometry ...

> As geometries, to me
> they look like polygons and I didn't realize they were in fact holes.
>
> See
> this example in QGIS for face_id 2418 turned into a polygon:
> https://dl.dropboxusercontent.com/u/5196336/geomFromFace_id_2418.bmp. Now,
> look at its mbr representation:
> https://dl.dropboxusercontent.com/u/5196336/mbrFace_id_2418.bmp. How do
> you tell it's a hole and not just a regular face?

The role of a face can be different for different TopoGeometry objects.
It could be an hole (Vatican state for the Rome polygon) or a polygon
(Vatican state for the Vatican state).

A face may be bounded by a single ring (having no holes) or by multiple
rings (a shell + a number of holes).

+-------+ F1 is bounded by 2 rings (shell and 1 hole)
| F1 | F2 is bounded by 1 ring (only shell)
| +---+ |
| |F2 | |
| +---+ |
+-------+

It's not "unioned with face 950", they are still distinct faces,
but they both partecipate in the definition of TopoGeometry 814.

Schematic example:

+---------------+
| TG 814 |
+---------------+

composed by
|
v

+-------+-------+
| F950 | F2418 |
+-------+-------+

>> I suggest you load the TopoGeometry layer with qgis and see for yourself.
>
> When you say load the TopoGeometry layer, I assume you mean the topogeom
> column I populated with:
>

> UPDATE syshiera.de_20k SET topogeom = toTopoGeom(geom_32198,
> 'de_20k_topo', 1, 2.0)

Right

> When I load it in QGIS using the PostGIS table loader, I can visualize it
> and see there are holes. E.g.
> https://dl.dropboxusercontent.com/u/5196336/holeInTopogeometry.bmp

Nice. Was there an hole in the original layer (geom_32198) ?
Or could it be you DELETEd one face too much from the relation table ?
In any case you can always re-create the TopoGeometry with another
UPDATE. Find the gid of the missing geometry and do something like:

UPDATE syshiera.de_20k SET topogeom = toTopoGeom(geom_32198,

'de_20k_topo', 1, 2.0) WHERE gid = :the_gid;

>> Unfortunately (that'd be useful to add in qgis) it won't be easy to
>> find TopoGeometry 1546 as extracting the id of a TopoGeometry requires
>> running an expression on the database: id(topogeom).
>
> I can either
/
> SELECT *
> FROM syshiera.de_20k
> WHERE id(topogeom) = 1546;
/
>
> or, as I mentioned before, I can superimpose my 'faces' layer, which I
> derived from de_20k_topo.face
> using ST_GetFaceGeometry, to 'topogeom' to get the topogeo_id. No?

Yes

>> You could create a view selecting just that topoGeometry, or a subset
>> of them, or you could just load them all, hoping it won't be too
>> expensive
>> (another spot where qgis support for postgis topology could be improved).
>
> I agree, just loading the TopoViewer freezes QGIS for a while and I don't
> lack ressources on my machine.

I have some ideas to speed that up, but no time/fundings to realize it.
It would be about allowing for layer definition to contain a parametrized
SQL query so that depending on the extent/scale a different query can
be run..

> Now it's clearer in my head what I need to do, thanks to you!
>
>> It looks like you're on the right track so far.
>
> Will keep trying and I will win. Again, thanks so much for your help.

Please consider helping back by improving the available tools to deal
with PostGIS Topology.

The QGIS PostGIS Topology Editor is very young and could be added some
tools for resolving TopoGeometry overlaps and gaps. Just takes some
python code and stubborness :)

--strk;

() Free GIS & Flash consultant/developer
/\ http://strk.keybit.net/services.html

[Guillaume Drolet]

Sandro,

Using your help, I filled most holes but three I couldn't fill using the
"UPDATE ... SET topogeom ..." approach. Will explore other approaches later.

Re the presence of holes, I may have deleted faces when running my function
for removing extra nodes. This needs checking.

Advancing in my work, I also got to drop edge 6341 and I end up with this:

https://dl.dropboxusercontent.com/u/5196336/edge6340.bmp

So the next step is to split edge 6340 about its middle, with the goal of
then connecting the newly created node with node 4753, thus creating an edge
where there should be one (i.e. on the light grey line).

However, I get the error: *point not on edge* when running ST_ModEdgeSplit.
First, I tried to find the new node location using the coordinate capture
tool in QGIS but because I got that error I thought maybe I really need to
find a point that's on the edge so I used this approach:

/WITH edge AS (SELECT geom line
FROM de_20k_topo.edge_data
WHERE edge_id = 6340
), point AS (SELECT e.line, ST_ClosestPoint(e.line,
ST_GeomFromText('POINT(-208737.346 811681.163)', 32198)) point
FROM edge e
)
SELECT ST_ModEdgeSplit('de_20k_topo', 6340, point)
FROM point;

/

I still get the same error: point not on node.

How is one supposed to find a point that IS on an edge?



--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007281.html

Sent from the PostGIS - User mailing list archive at Nabble.com.

[Rémi Cura]

Just to answer about Grass.
I'm not a Grass user, and I do use PostGIS topology every day (heart of my work)
i know only 3 commands :
v.in.ogr dsn="PG:host=localhost port=5433 dbname=... user=... password=..." layer=schema.table out=layer_name
v.clean ...
v.out.postgis -l -2 input=layer_name dsn="PG:host=localhost dbname=... port=5433 user=... password=..." olayer=my_postgis_topo_layer -t

From a practical standpoint,
I like to choose the best tool for each job !

the focus of PostGIS topology is to store a topology model very close the standart and allow easy conversion to geometry.

It is not a tool designed to clean an input topology, nor merge, batch change, etc. Similarly, it is very slow to compute the topology from geometry. 

On the opposite Grass is the only GIS entirely topological since the beginning. This means that they have powerfull and honed topological capabilities, because topology is not an option in this soft.

However it is not a tool designed to get easy access to data like postgis topology.

I prefer to chose the best tool for each job !

advanced topo computing/cleaning/batch operations : Grass

topo storing/relational usage/conversion to geom : PostGis topology.

Clean correct topological operations with a given precision is something incredibly difficult, so don't ask too much to PostGIS topology.

Cheers,

Rémi-C

[Sandro Santilli]

On Wed, Nov 05, 2014 at 03:37:32PM -0800, Guillaume Drolet wrote:

> So the next step is to split edge 6340 about its middle, with the goal of
> then connecting the newly created node with node 4753, thus creating an edge
> where there should be one (i.e. on the light grey line).
>
> However, I get the error: *point not on edge* when running ST_ModEdgeSplit.

[...]

> How is one supposed to find a point that IS on an edge?

This is indeed tricky. The easiest way could be adding an explicit node
with TopoGeo_addPoint and a tolerance value large enough. Something like:

SELECT
TopoGeo_addPoint('de_20k_topo',
ST_ClosestPoint(
e.line,
ST_GeomFromText('POINT(-208737.346 811681.163)', 32198)
),
1.0 -- tolerance
)
FROM edge e
WHERE edge_id = 6340;

The query will return the id of the nelwy added node, which would
split edge 6340 if it's closer than 1.0 units from it.

You can verify if the split happened by checking if either
start_node or end_node of edge 6340 are equal to the returned id.

--strk;

() Free GIS & Flash consultant/developer
/\ http://strk.keybit.net/services.html

[Guillaume Drolet]

Sandro: TopoGeo_addPoint sounds like the appropriate choice. I hadn't seen
this alternative topology function. Will try that, thanks.

Remi: Thanks for your advice on when to use what tools. I should give GRASS
a shot at topology cleaning. Since topology-related work isn't something I
do everyday and PostGIS is a tool I use frequently, I naturally went in that
direction.



--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007284.html

Sent from the PostGIS - User mailing list archive at Nabble.com.

[Guillaume Drolet]

I finally finished cleaning my topology, using a mix of SQL commands,
functions, and also using Sandro's QGIS topology editor plugin.

In parallel, I also tried going the GRASS route, importing my original
shapefile in GRASS and applying the cleaning function (v.clean ...) as
suggested by Remi. I injected the result into PostGIS and created a
topology. Most of the 23 slivers were there in the topology so one of two
things: 1) GRASS didn't fix them (maybe I didn't choose an appropriate
tolerance) or 2) they were created by PostGIS when building the topology.

Anyway, I got what I finally got what I need, thanks a lot to all of you who
helped me.

I'm gonna ask you one last thing if I may: I want to replace the original
geom column (i.e. that used to build the topogeom) with a topologically
correct one, and keep all the associated ecological attributes. Is this the
right way to do it:

UPDATE syshiera.de_20k SET geom = topogeom::geometry; ?

Thanks.



--
View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007286.html

[Sandro Santilli]

On Fri, Nov 07, 2014 at 01:14:31PM -0800, Guillaume Drolet wrote:
> I finally finished cleaning my topology, using a mix of SQL commands,
> functions, and also using Sandro's QGIS topology editor plugin.

Congratulations!

> In parallel, I also tried going the GRASS route, importing my original
> shapefile in GRASS and applying the cleaning function (v.clean ...) as
> suggested by Remi. I injected the result into PostGIS and created a
> topology. Most of the 23 slivers were there in the topology so one of two
> things: 1) GRASS didn't fix them (maybe I didn't choose an appropriate
> tolerance) or 2) they were created by PostGIS when building the topology.

I can't help here as I've never used GRASS for that yet.
Would be interested in knowing more though. Maybe Remi can help there.

> I'm gonna ask you one last thing if I may: I want to replace the original
> geom column (i.e. that used to build the topogeom) with a topologically
> correct one, and keep all the associated ecological attributes. Is this the
> right way to do it:
>
> UPDATE syshiera.de_20k SET geom = topogeom::geometry; ?

Yes, that gives you a topologically correct "snapshot" of current
topology state (assuming you have no more unassigned or multiply-assigned
faces).

Note that there might still be cases in which the universal face enters
between two faces creating small gaps. I guess a way to check that would
be by computing angle between edges around each node and ensuring none
is below a given threshold (just ideas for more development to improve
coverage of useful functions ;)

--strk;

() Free GIS & Flash consultant/developer
/\ http://strk.keybit.net/services.html

[Humberto Cereser Ibanez]

I was with the same problem, getting gaps in some adjacent polygons
while trying to simplify polygons with Postgis.
Thank Guillaume / Rémi / Sandro by the contributions.
Illuminated by this thread, I created a script in which I avoided gaps
from the bad polygons. Really the most appropriate tool, GRASS, became
this work easier.

On Fri, 2014-11-07 at 13:14 -0800, Guillaume Drolet wrote:
> I finally finished cleaning my topology, using a mix of SQL commands,
> functions, and also using Sandro's QGIS topology editor plugin.
>
> In parallel, I also tried going the GRASS route, importing my original
> shapefile in GRASS and applying the cleaning function (v.clean ...) as
> suggested by Remi. I injected the result into PostGIS and created a
> topology. Most of the 23 slivers were there in the topology so one of two
> things: 1) GRASS didn't fix them (maybe I didn't choose an appropriate
> tolerance) or 2) they were created by PostGIS when building the topology.

Simplifying polygons with QGIS and GRASS
Note: SRID = 4326
1) Open shapefile in QGIS
2) Create New Mapset
Import Extent from opened Shapefile
3) Check the smaller area of the polygons involved
select ST_Area (geo_polygon) from table order by
ST_Area(geo_polygon);
Smaller area = 0.005471
4) Import shapefile into GRASS
v.in.ogr dsn = /home/humberto/shapefiles/pr/41MUE250GC_SIR.shp
output = prGrass snap = 0.0001 min_area = 0.004 -o
Settings of parameters: Snapping threshold for boundaries = 0.0001
and Minimum size of area to be imported (squere units) = 0.004
Note: The 0.004 value was chosen to allow importing of the smaller
area = 0.005 and all areas greater than it, avoiding any small ring that
would be arising from this importing.
5) Simplify the Shapefile polygon imported
v.generalize input=prGrass@elimSlivers type=area layer=1 -c
type=boundary method=douglas threshold=0.001 look_ahead=7 reduction=50
slide=0.5 angle_thresh=3 degree_thresh=0 closeness_thresh=0
betweeness_thresh=0 alpha=1.0 beta=1.0 iterations=1
output=prGrassGeneralized
Settings of parameters: threshold = 0.001
Note: the threshold value 0.001 was chosen by trial
6) Export the simplified polygon to Postgis
v.out.ogr.pg.py input=prGrassGeneralized@elimSlivers type=area
layer=1 olayer=prgeneralized database=simp22 host=127.0.0.1 port=5432
user=name password=nonono

>
> Anyway, I got what I finally got what I need, thanks a lot to all of you who
> helped me.
>
> I'm gonna ask you one last thing if I may: I want to replace the original
> geom column (i.e. that used to build the topogeom) with a topologically
> correct one, and keep all the associated ecological attributes. Is this the
> right way to do it:
>
> UPDATE syshiera.de_20k SET geom = topogeom::geometry; ?
>
> Thanks.
>
>
>
> --
> View this message in context: http://postgis.17.x6.nabble.com/Topology-cannot-delete-slivers-or-gaps-tp5007250p5007286.html
> Sent from the PostGIS - User mailing list archive at Nabble.com.
> _______________________________________________
> postgis-users mailing list
> postgi...@lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Thanks,

Humberto

[Rémi Cura]

Note you could have directly exported into postgis topology, which is slow but very convenient (stil faster than postgis topology conversion).

I'm not a grass user,

from what I understood you have to use the "-l" switch in the command "v.out.postgis"

Glad it helped =)

Cheers,
Rémi-C

[Humberto Cereser Ibanez]

Hi Rémi,

On Wed, 2014-11-12 at 10:51 +0100, Rémi Cura wrote:
> Note you could have directly exported into postgis topology, which is
> slow but very convenient (stil faster than postgis topology
> conversion).
>
> I'm not a grass user,
>
> from what I understood you have to use the "-l" switch in the command
> "v.out.postgis"

Thanks for your tip. It will be useful for make good use of Postgis
topology.

>
>
> Glad it helped =)
>
> Cheers,
> Rémi-C

Cheers,

Humberto

[Sandro Santilli]

On Wed, Nov 12, 2014 at 10:51:34AM +0100, Rémi Cura wrote:
> Note you could have directly exported into postgis topology, which is slow
> but very convenient (stil faster than postgis topology conversion).
> I'm not a grass user,
> from what I understood you have to use the "-l" switch in the command "

> *v.out.postgis*"

I'm trying this out today, and indeed is very slow, much slower than
building the topology, which is surprising to me.

A look at database activity shows that the command (v.out.postgis -l)
is writing _simple_ features in the output layer. Do you have an idea
about how to skip that part ? (TopoGeometry would be enough for PostGIS
to reconstruct the SFS on demand).

--strk;

[Rémi Cura]

This is suprising indeed.

Im my test using the v.out.postgis was still way faster than constructing with postgis topology (about x3 I think).

Are you using v.out.postgis -l, with grass 7?

Grass is a topology that is more node centered, opposite to postgis topology which is more edge centered.

It is possible to output a node centered topology with grass as text (v.info -t),

but I found it unpractical to try to convert this into a propoer postgis topology with pure SQL.

Now I think I should have used a python way from within grass (grass offers full binding, so it should be doable to loop trough grass internal topology and convert it to postgis topology)

I tried to reach grass dev several time to try to improve the postgis export, I had no answer.

My intuition is that the slowness comes from the fact that grass tries to write  edge by edge , and each time the edge_data trigger are triggered.

The correct way would be to compute full topology structure, then deactivate trigger, then fill table , then reactivate trigger.

Cheers,
Rémi-C

[Sandro Santilli]

On Tue, Nov 18, 2014 at 08:26:29PM +0100, Rémi Cura wrote:
> This is suprising indeed.
> Im my test using the v.out.postgis was still way faster than constructing
> with postgis topology (about x3 I think).
>
> Are you using v.out.postgis -l, with grass 7?

Yes.

The topology was pretty big, maybe that's the difference.
Here's the dataset (not sure if still accessible):
http://lists.osgeo.org/pipermail/postgis-devel/2014-January/024085.html

GRASS 7.0.0beta3 (world):~ > v.info -t map=million_poly_topo1
nodes=3247661
points=0
lines=0
boundaries=5102420
centroids=1161248
areas=1855936
islands=1177
primitives=6263668
map3d=0

> Now I think I should have used a python way from within grass (grass offers
> full binding, so it should be doable to loop trough grass internal topology
> and convert it to postgis topology)
>
> I tried to reach grass dev several time to try to improve the postgis
> export, I had no answer.

That's unfortunate.

> My intuition is that the slowness comes from the fact that grass tries to
> write edge by edge , and each time the edge_data trigger are triggered.
> The correct way would be to compute full topology structure, then
> deactivate trigger, then fill table , then reactivate trigger.

Or maybe it has to build the next_{left,right}_face info too.
Enabling PostgreSQL logging should tell (but better not done with
a ~2 million polygons topology).

Oh how better it would be to speed up topology building in PostGIS :)

--strk;

() Free GIS & Flash consultant/developer
/\ http://strk.keybit.net/services.html

[Rémi Cura]

Hey Sandro,
in my opinion we don't need a better speed for current construction, we need a batch constructing mode !

Optimally it could be done by using the geos, because it constructs all the topology very very fast.

Adding one feature is actually quite fast, even on already big topology.

Its when you want to add a lot's that it becomes increasingly slow (maybe because indexes are not updated,or because we are in one transaction?)
The slowing seems to be very non linear, probably following n^2, where n is the number of feature already constructed in the transaction.

Cheers,
Rémi

[Sandro Santilli]

On Wed, Nov 19, 2014 at 12:50:09PM +0100, Rémi Cura wrote:
> Hey Sandro,
> in my opinion we don't need a better speed for current construction, we
> need a batch constructing mode !

This is for sure, as the current one just does the intial noding
but finally invokes the incremental functions.

> Optimally it could be done by using the geos, because it constructs all the
> topology very very fast.

Possibly, yes.
Then I guess attributes association would be left to a second step...

> Adding one feature is actually quite fast, even on already big topology.
>
> Its when you want to add a lot's that it becomes increasingly slow (maybe
> because indexes are not updated,or because we are in one transaction?)
> The slowing seems to be very non linear, probably following n^2, where n is
> the number of feature already constructed in the transaction.

An issue with index use was recently fixed.
There might be another one hiding somewhere.

--strk;

[Sandro Santilli]

On Wed, Nov 19, 2014 at 04:47:48PM +0100, Sandro Santilli wrote:
> On Wed, Nov 19, 2014 at 12:50:09PM +0100, Rémi Cura wrote:
> >
> > Adding one feature is actually quite fast, even on already big topology.
> >
> > Its when you want to add a lot's that it becomes increasingly slow (maybe
> > because indexes are not updated,or because we are in one transaction?)
> >
> > The slowing seems to be very non linear, probably following n^2, where n is
> > the number of feature already constructed in the transaction.
>
> An issue with index use was recently fixed.
> There might be another one hiding somewhere.

On a closer look, I'm thinking the single-transaction is what commonly
hits during topology building (UPDATE .. SET tg = toTopoGeom ..)

Starting from an empty topology and running a single statement
invoking toTopoGeom for each of many inputs result in no stats ever
being visible by the planner within the transaction. In turn this
is likely to opt for sequencial scans (an empty table is quicker to
scan sequencially).

This would explain why populating in chunks works better, using
a transaction for each chunk
(UPDATE .. SET .. WHERE gid >= N AND gid < N+chunksize)

It could be interesting to try a wrapper function taking care of
running ANALYZE on the primitive tables every N calls to toTopoGeom
(or N primitives being created, regardless of number of simple inputs).

--strk;

() ASCII ribbon campaign -- Keep it simple !
/\ http://strk.keybit.net/rants/ascii_mails.txt

[Rémi Cura]

Sadly I think it can't be done with pure plpgsql,
because every function is wrapped in a transaction no matter what.

You can only do it using the trick to connect to the same database from within with the extension "dblink"

But I find difficult to understand how transactions and sub transactions affects performance.

Morevover , the transaction think is not the only problem. It is more a design problem.
Even with CGAL, building topology one by one or with the batch mode changes radically the time of building ( n to n^2 at least).

That's why I truly think perf is going to come from a batch mode. Tweaking the current process is just damage control in my opinion.

This is not so hard do to if we rely a bit on GEOS.

1. cut the input geom into a space partition (for line, ST_Node, for poly ST_Polygonize)

2. populate node table, and create a temp table with list of line for each node

3. Populate edge_data

4. fill next / left for edge_data

5. compute area (Polygonize, Geos? )

6. Map the input geom to generated topology (to be able to use attributes)

I already tested 1,2,3,6.

It can be fast (not to that building full topo in geos and converting it to postgis_topology I'm afraid), and it will scale very well.

Cheers,

Rémi-C