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

[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 | |
 | +---+ |
 +-------+


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