Avoid intersection of new polygons.. IS BROKEN

[pinux]

In the last few month I have used a lot the "Avoid intersection of new polygons.." function i qgis. Yesterday I needed to do merge all the polygons in a postgis table, I tried it with the merge tool in qgis and it was not working. So I tried it in postgis And the result was one polygon with a lot of holes. I understood that the "Avoid intersection of new polygons.." function was not working well. After I corrected all the errors on the polygons in may table I was able to merge all the polygons even with qgis. My next step was to create a test table with a set of test views to seen all the error that the "Avoid intersection of new polygons.." function is creating in qgis.

the test table:

CREATE TABLE av_intersection (

gid int UNIQUE,
t_name varchar(50),
PRIMARY KEY (gid)

) ;
SELECT AddGeometryColumn ('public', 'av_intersection', 'the_geom', '900913', 'MULTIPOLYGON', 2);
CREATE INDEX av_intersection_gidx ON av_intersection USING GIST (the_geom) ;

the union view:

CREATE OR REPLACE VIEW av_intersection_union AS
SELECT
min(p.gid) as gid,
t_name,
st_union(p.the_geom) AS the_geom
FROM av_intersection p
GROUP BY p.t_name;

view to see if there are multiple rings:

CREATE OR REPLACE VIEW av_intersection_multiple_ring AS
SELECT
gid,
the_geom,
(ST_NRings(the_geom) - ST_NumInteriorRings(the_geom)) as nrings
FROM av_intersection
WHERE ST_NRings(the_geom) > (1 + ST_NumInteriorRings(the_geom));

and a view with all the polygons that are overlapping:

CREATE OR REPLACE VIEW av_intersection_overlaps AS
SELECT DISTINCT
(a.gid*2+b.gid*10) as gid,
a.gid a_gid,
b.gid b_gid,
a.the_geom
FROM
av_intersection a,
av_intersection b
WHERE
a.gid != b.gid AND ST_Overlaps(b.the_geom, a.the_geom);

The result with 11 polygons:

• union with 2 holes
  
• 9 polygons are overlapping
  
• 2 polygon with multiple ring
  

The result of the test you can see in the attachment with all the code to test yourself and the insert command for the polygons present in my test.

I have done the test on 1.5 but the problem is even in trunk.

[mhugent]

I cannot entirely follow your description. Do you mean that the merge tool and avoid intersection tools don't work with invalid geometries? And after you fixed the geometry errors on the polygons, it did work? If yes, that's the behaviour I would expect. Because merge and avoid intersection use geos and I'm pretty sure it does not do union or difference on invalid geometries. Maybe postgis tests the geometries first and tries to fix them automatically?

[pinux]

Replying to mhugent:

I cannot entirely follow your description. Do you mean that the merge tool and avoid intersection tools don't work with invalid geometries? And after you fixed the geometry errors on the polygons, it did work? If yes, that's the behaviour I would expect. Because merge and avoid intersection use geos and I'm pretty sure it does not do union or difference on invalid geometries. Maybe postgis tests the geometries first and tries to fix them automatically?

My problem is not with the merge, when I was trying to merge polygons that I had digitized with the option "Avoid intersection of new polygons" I had some errors. So i have made some view to understand what was happening with my data. And I have found that the option "Avoid intersection of new polygons" is not working correctly:

1. If I do an union in postgis I have 2 holes,
   
2. Polygons are overlapping;
   
3. Polygons with multiple ring, even if I don't have digitized polygons with multiple rings.
   

If you load the sql that I have uploaded you will see what I'm talking about.
The result with 12 polygons:

• union with 2 holes
  
• 9 polygons are overlapping
  
• 2 polygon with multiple ring
  

[mhugent]

I'm having problems importing the attached .sql into postgres: ERROR: AddGeometryColumns() - invalid SRID (Seems my db does not know the google crs)

Could you attach shape or spatialite? The individual (overlaping) polygons should be enough. All I need is a set of polygons to replay the problematic intersection steps.

[pinux]

Replying to mhugent:

I'm having problems importing the attached .sql into postgres: ERROR: AddGeometryColumns() - invalid SRID (Seems my db does not know the google crs)

Could you attach shape or spatialite? The individual (overlaping) polygons should be enough. All I need is a set of polygons to replay the problematic intersection steps.

I have uploaded a new version of the file with the SRID -1.

[mhugent]

Ok, I can reproduce the bug now.

[mhugent]

This was a rounding problem, should be fixed in r15087. Please test and reopen the bug if the problem still persist.

[pinux]

The problem is still present

[pcav]

Is this really so bad to be blocking for 1.7?

[pavolsk]

Yes, it really is. I spent few hours trying to understand, where can be the problem. Apparently this function itself is not broken, removePolygonIntersections function is broken (something in the internal logic of this - something inside is antiGEOS). More or less, it is not really necessary to have function like this. Complex digitising tasks can be solved using polygonize function (polygonize plugin, polygonize in PostGIS), It would be great to have own polygonize function in QGIS core. Bugs like this are the worst promotion for digitising eager newcomers to qgis. I think, this function and removePolygonIntersections function should be droped in the new release.

[mayeulk]

Hi,

I can confirm most of the above, with some comments which might help debugging.

You can try this query:

select *, st_area(geom), st_isvalid(geom)  from (select a.gid + b.gid * 100 AS gid, a.gid AS a_gid, b.gid AS b_gid,
st_setsrid((st_intersection(a.the_geom, b.the_geom)), 4326) as geom
FROM av_intersection a, av_intersection b
  WHERE a.gid < b.gid AND st_overlaps(b.the_geom, a.the_geom)) as the_intersects;

You'll see that 2 thirds of the intersections have a null area. For the other ones, their area is always (on my machine) 0.00048828125 and I believe it is a Postgis/geos error, see below:

WITH av_intersect_911 as (
select 91104 AS gid,
st_setsrid(ST_GeometryN((st_intersection(a.the_geom, b.the_geom)),4), 4326) as geom
FROM av_intersection a, av_intersection b
  WHERE a.gid=9 and b.gid=11
limit 1)

SELECT gid, geom, GeometryType (av_intersect_911.geom), st_area(av_intersect_911.geom), ST_AsEWKT (av_intersect_911.geom)

FROM av_intersect_911;

The above query shows a polygon which is (part of) the intersection of polygons 9 and 11. It's area is 0.00048828125. This polygon has this EWKT definition:

"SRID=4326;POLYGON((
1517034.97851789 5634807.73109361,
1517034.97851789 5634807.73109361,
1517034.11662417 5634808.29056848,
1517034.97851789 5634807.73109361))"

Since 3 of its four points are identical, it is trivial to say that this polygon is flat and its area is zero, hence the area computation is false.

I am not 100% sure that there are "real" intersections, as several of the intersections are not polygons but lines at the boundary of 2 polygons. Still, I agree that the topology is incorrect: the topology is broken.

Still, this does not corrupt existing data nor create invalid geometry: they can be displayed in qgis and the following is true in postgis:

select st_isvalid(the_geom) from av_intersection;

Given this postgis bug, I tried to reproduce this with a shapefile (to make sure that we are hunting a QGIS bug, not a postgis bug).

Below an ASCCI art explanation of what I found, please see also attached screenshot and attached shapefile.

You have A: 
____
    |
A   |
-----

You add B:

____________
    |B     /
A   |     /
----|----/

Everything is OK until you add C:

   ____
  / C  \
_/______\___
    |B     /
A   |     /
----|----/

You expect C to be like this:
   ____
  / C  \
 /______\

The topological error that I can confirm is the following: C is like this:
   ____
  / C  \
 /______\
    |  
    |     

My conclusion:
When you add a polygon C where two polygons A and B are
touching (or are supposed to), then C will have a leg going
in between A and B. It is not clear whether A and B are 
strictly contiguous or if there is an (unwanted) hole between 
them.

[mayeulk]

Screenshot showing "legs" (wrong topology) in red (one selected polygon)

[mayeulk]

Shapefile illustrating wrong typology (see ​http://trac.osgeo.org/qgis/ticket/2921#comment:10)

[mhugent]

The removePolygonIntersction method from QGIS just calls difference() method from geos repeatedly.

Taking the attached shape and removing the last polygon, it seems to me that the polygons in that dataset are not really adjacent to each other (there are small legs between). You could test that with the union function in QGIS directly.

Is it possible those 'legs' have been introduced by adding polygons with the 'topological editing' mode in QGIS? That's where r15087 enhanced the precision. Now I don't manage to introduce legs like that anymore. If you do, please send me a shapefile, a screenshot how to add the new polygon and the settings you used for it (e.g. snapping tolerance) used.

[mhugent]

Wait, I managed to reproduce the mentioned issue with the test shape. So probably using a more sensitive rounding tolerance should be enough.

[mayeulk]

I can reproduce the bug on the shapefile I attached with those parameters:

For the shapefile I attached:

tolerance: 0.000001 (in map units)

Mode:to vertex and segment

It is still unclear to me whether there is originally a hole between A and B, or if they are OK and only C is wrong. Is there any reliable way to check? Visible inspection is useless even at scale>1:1 and postgis does not seem reliable here.

[mhugent]

You could check if the polygons are really have a common boundary with the merge edit tool in QGIS.

Could you try if it works with r15741 (creation of the whole polygon collection, not only the last step)?

Thanks, Marco

[mayeulk]

Hi Marco, I will look at it. In the meantime, I found this:

same bug with tolerance: 0.0001(in map units)

same bug with tolerance: 0.1(in map units)

However, with tolerance: 1 (in map units, Mode:to vertex and segment ), then there is snapping during digitization of each point, and with snapping I could not reproduce the bug anymore with r15718. (In fact, snapping makes the set operation unnecessary and probably the bug cannot appear).

[mhugent]

enabling 'topological editing' also plays an important role. With it enabled, qgis inserts new vertices and the bug is less likely to appear.

[mayeulk]

Hi,

Unfortunately the bug is still there with r15741 (with 0.000001 as snapping option in map units).

And I got as often as earlier the following message:

"The feature could not be added because removing the polygon intersections would change the geometry type."

(especially when trying to create a polygon over areas with many polygons and where many polygons touch). Which I guess might be due to some holes in the polygons: the difference create multipart polygons, hence the message (postgis is more flexible than the shapefile format here and allows the buggy geometries to be created).