qgis / QGIS

QGIS is a free, open source, cross platform (lin/win/mac) geographical information system (GIS)
https://qgis.org
GNU General Public License v2.0
10.63k stars 3.01k forks source link

The QGIS intersection tool yields incomplete results #22477

Closed qgib closed 6 years ago

qgib commented 8 years ago

Author Name: Mattias Lindman (Mattias Lindman) Original Redmine Issue: 14504 Affected QGIS version: 2.12.3 Redmine category:processing/qgis


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 #21307 may be related.



Related issue(s): #20192 (relates), #21307 (relates), #22572 (relates) Redmine related issue(s): 11986, 13246, 14606


qgib commented 8 years ago

Author Name: Maximilian Krambach (Maximilian Krambach)


Related to #22469.

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."[1] 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"

[1]http:// geos.refractions.net/ro/doxygen_docs/html/classgeos_1_1geom_1_1Geometry.html#a44

qgib commented 8 years ago

Author Name: Maximilian Krambach (Maximilian Krambach)


Geos has such a function, geos::algorithm::LineIntersector:computeIntersection.

https://github.com/qgis/QGIS/blob/fab8dc21762b8ca0a9a48ffffe343788fddabc24/src/core/geometry/qgsgeometry.cpp#L1381 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?

qgib commented 8 years ago

Author Name: Jürgen Fischer (@jef-n)


Maximilian Krambach wrote:

Geos has such a function, geos::algorithm::LineIntersector:computeIntersection.

is it exposed in the C-API?

qgib commented 8 years ago

Author Name: Maximilian Krambach (Maximilian Krambach)


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.

qgib commented 8 years ago

Author Name: Maximilian Krambach (Maximilian Krambach)


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.

qgib commented 8 years ago

Author Name: Mattias Lindman (Mattias Lindman)


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

/Mattias

qgib commented 8 years ago

Author Name: Claas Leiner (@claasleiner)


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.

Best wishes

Claas



qgib commented 8 years ago

Author Name: Claas Leiner (@claasleiner)


When I use the files clip.py and intersection.py from the processing 2.12 (Https://plugins.qgis.org/plugins/processing/version/2.12.2/) in the processing 2.14 everything works well. (copy in the folder */python/plugins/processing/algs/qgis).

Claas

qgib commented 8 years ago

Author Name: Giovanni Manghi (@gioman)


see also #20192 for more intersection problems.


qgib commented 8 years ago

Author Name: Sandro Santilli (@strk)


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.

Example:

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 ?

qgib commented 7 years ago

Author Name: Giovanni Manghi (@gioman)


qgib commented 6 years ago

Author Name: Alexander Bruy (@alexbruy)


Should be fixed in master/3.2. Please check if it works now.


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 #21307 may be related. to 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 #21307 may be related.

qgib commented 6 years ago

Author Name: Nyall Dawson (@nyalldawson)


Issue is fixed in master, please open a new ticket if new issues are encountered