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.48k stars 2.99k forks source link

Overlay union: errors and data loss #31552

Open aborruso opened 5 years ago

aborruso commented 5 years ago

Hi, this issue to make a summary of two others (#31528 and #31541), opened by @pigreco

These are the files I have used to make my tests. Essentially two shapefiles (input_01 and input02), both projected in EPSG:4326 and EPSG:3004, so in total four shapefiles.

My goal is to apply to them Overlay/Union, without using overlay layer and producing a temp output (as below).

image

input_01 set

If I run Overlay/Union to this set I have:

input_02 set

If I run Overlay/Union to this set I have:

You can see all these steps in this video https://www.youtube.com/watch?v=y7EI9aIJX6U

cc @wonder-sk @gioman

Thank you

My QGIS

QGIS version
3.8.2-Zanzibar
QGIS code revision
4470baa1a3
Compiled against Qt
5.11.2
Running against Qt
5.11.2
Compiled against GDAL/OGR
2.4.1
Running against GDAL/OGR
2.4.1
Compiled against GEOS
3.7.2-CAPI-1.11.0
Running against GEOS
3.7.2-CAPI-1.11.0 b55d2125
PostgreSQL Client Version
10.8
SpatiaLite Version
4.3.0
QWT Version
6.1.3
QScintilla2 Version
2.10.8
Compiled against PROJ
5.2.0
Running against PROJ
Rel. 5.2.0, September 15th, 2018
OS Version
Windows 10 (10.0)
gioman commented 5 years ago

input_01 set

If I run Overlay/Union to this set I have:

* for 4326 one

  * error message
  * data loss

* for 3004 one

  * no error message
  * no data loss

input_02 set

If I run Overlay/Union to this set I have:

* for 4326 one

  * error message
  * no data loss (mine is only a visual check)

* for 3004 one

  * an error message, but of a different kind
  * no data loss (mine is only a visual check)

my observations (on 3.8.2/Linux) are slightly different:

in case of input_01 in 4326, error message, no data loss.

in case of input_02 in 4326, no error message, no data loss.

in case of the input_02 the "error message, but of a different kind" is

GEOS geoprocessing error: difference failed.
TopologyException: found non-noded intersection between LINESTRING (2.36891e+06 4.22927e+06, 2.36892e+06 4.22927e+06) and LINESTRING (2.36892e+06 4.22928e+06, 2.3689e+06 4.22926e+06) at 2368909.6404680978 4229266.2896102117
Execution failed after 0.14 seconds

but I guess this could be very specific for this dataset, and not a general problem.

wonder-sk commented 5 years ago

If there is a topology exception emitted from GEOS, that means some XY coordinate values are very very close to each other. Because of that, some math operations give wrong results.

The proper solution to this problem is to run "Snap geometries to layer" processing algorithm in "Snap to anchor nodes (single layer only)" mode with a small tolerance (note it is in map units!). Afterwards run "Fix geometries" algorithm to make sure that if snapping produced invalid geometries, they will get fixed.

A good question is whether QGIS should try harder to help users in these scenarios:

aborruso commented 5 years ago

The proper solution to this problem is to run "Snap geometries to layer" processing algorithm in "Snap to anchor nodes (single layer only)" mode with a small tolerance (note it is in map units!). Afterwards run "Fix geometries" algorithm to make sure that if snapping produced invalid geometries, they will get fixed.

Hi @wonder-sk I have tested it with input_01 set and it works, thank you.

For input_02-3004 shapefile instead, applying a snap coarse tolerance (1 meter) and than "Fix geometries", I still have the same GEOS error.

I agree with you: when there is a GEOS error, the user should have a note (with a URL for more detailed descriptions). Your note is perfect.

gioman commented 5 years ago

For input_02-3004 shapefile instead, applying a snap coarse tolerance (1 meter) and than "Fix geometries", I still have the same GEOS error.

"fun" fact, here doing the snapping and fix operations produced an input that when using it to Union in single layer mode does not show anymore the GEOS error, instead is shows an error like

Feature could not be written to Union_a52b9760_1f5a_42ae_b38d_e23910ed1223
Execution completed in 0.13 seconds
Results:
{'OUTPUT': 'Union_a52b9760_1f5a_42ae_b38d_e23910ed1223'}

and the output is created anyway, but with missing parts.

We should try to get to the bottom of this mess.