Bug report #14504
The QGIS intersection tool yields incomplete results
|Affected QGIS version:||2.12.3||Regression?:||No|
|Operating System:||Easy fix?:||No|
|Pull Request or Patch supplied:||No||Resolution:||fixed/implemented|
|Crashes QGIS or corrupts data:||No||Copied to github as #:||22477|
In a land cover study I want to derive country-wise land cover units in order to determine the area of each land cover type in each country. The source data is a land cover raster that I have vectorized and then clipped by a layer containing the national boundaries of the three countries in the study.
I then use Vector/Geoprocessing Tools/Intersect with the vectorized and clipped land cover layer and the national boundary layer as inputs. As I have made sure that both layers cover exactly the same area the result of the intersection should be land cover units divided by country whose total area should be the same as the total area of all land cover units prior to the intersection.
It turns out that this is only the case for a few land cover types. For most of them the total area after the intersection is significantly lower than the original area, indicating that the intersection operation does not yield complete results. As an example the total shrubland area is reduced with about 35 % in the intersection operation.
In order to provide a possibility to reproduce the problem I attach an archive containing the national boundary layer and an extract from the land cover layer containing shrublands only.
I have tested the analysis using both QGIS 2.8.7 and 2.12.3 and the problem is the same in both versions. I am running QGIS using Windows 8.
When checking already submitted bug reports it appears that #13246 may be related.
#1 Updated by Maximilian Krambach over 4 years ago
Related to #14494.
The processing tool calls "geom.intersection(tmpGeom)", which only works as expected if the intersecting lines share a vertex on intersection.
The problems trace back to geos::intersection, which only "Returns a Geometry representing the points shared by this Geometry and other." If the intersection is not at a vertex, nothing is returned.
I think there needs to be a API function "compare two line segments and if these intersect, return both as lines with a vertex on the intersection"
#2 Updated by Maximilian Krambach over 4 years ago
Geos has such a function, geos::algorithm::LineIntersector:computeIntersection.
should check if there is a geos:geometry:intersection, and if not, try geos:algorithm:LineIntersector:computeIntersection. Other functions, like makeDifference (in the same file) may profit from this check, too.
This is a bit of a design issue: Do two overlapping lines intersect? I can think of situations without intersection (bridges etc.), but for every day usage, most lines DO intersect. Maybe some switch ("only intersect if there is a vertex on intersection") could be added?
#4 Updated by Maximilian Krambach over 4 years ago
Sorry, I was wrong. Upon further investigation, Intersection seems not to be the culprit.
However, only polygons1 that touch the outer line of the intersection polygon2 are affected, (but not all of them). I zoomed all the way in, and they still touch.
#5 Updated by Maximilian Krambach over 4 years ago
I looked further on a failing geometry1 (FID 2883 in shrublands). This geometry1 has no vertices in common with the geometry2 it intersects, and returns NoGeometry on intersection, even if it is the only feature in a shapefile.
If I snap two vertices1 to their respective nearest vertex2, the intersection returns a multipolygon.
#6 Updated by Mattias Lindman over 4 years ago
Hello Maximilian and Jürgen
Thank's for your input into this matter.
I think the problem originates from the fact that I first made a clip of the land cover layer with the same layer I later used for the intersection (NationalBoundaries). If I skip the clipping and just do the intersection of the original land cover layer the result is as it should - all land cover units that intersect the layer NationalBoundaries are obtained.
Back to the results of using the clipped land cover layer in the intersection. I find it strange that some but not all of the land cover polygons where a part of their boundary is identical with polygon boundaries in the intersection layer are lost in the intersection. It seems that this identity is the problem, but it does not appear to be general.
I looked into Maximilians comment that the polygon with FID 2883 does not have any vertices in common with the polygon it intersects in the intersection layer. I do not understand what you mean here - since the land cover layer (prior to intersection) is clipped with the NationalBoundaries layer it results in identical outer boundaries (i.e. identical vertices) in the two layers used in the intersection.
I do not know how the intersection algorithm works but I wonder if partly identical boundaries is a special case that the algorithm does not manage to handle, and if so, whether it is important to fix it. It is probably rare to first clip and then intersect with the same layer as you can do the intersection directly. But as I have done it this way it must mean that other users may do it too...
#7 Updated by Claas Leiner over 4 years ago
- File beispiel_intersect.7z added
The problem is in QGIS 2.14.1 still exists. In QGIS 2.8.7 it works.
Intersect form ftools (Vector > Geoprocessing > intersect) works in QGIS 2.14 fine. Intersect from the Toolbox lost Objekts.
The same is, when I use Clip.
Intersect and clip in the toolbox provides the sample data incomplete Ertgebnisse.
#10 Updated by Sandro Santilli about 4 years ago
About "Returns a Geometry representing the points shared by this Geometry and other." -- the GEOS statement has nothing to do with vertexes. Rather "points" is meant in the point-set theory.
strk=# select ST_AsText(ST_Intersection('POLYGON((0 5, 10 0, 10 10, 0 5))'::geometry, 'POLYGON((10 5, 20 10, 20 5, 10 5))'::geometry)); st_astext ------------- POINT(10 5) (1 row)
Are you sure your expectancies are correct ?