"Difference" produce incorrect result when using a point layer and a polygon overlay layer

[jfbourdon]

What is the bug or the crash?

Using points as an input layer and a polygon as an overlay layer with the "Difference" processing tool produces an incorrect result which is that points inside the polygon aren't shown on the map but are still present in the attribute table.

Probably related to:

• Union fails when using a point and a polygon layer #44801
• Algorithm "Difference" failing to execute properly #49291
• Fix difference algorithm when differencing a point results in a null geometry #49300

Steps to reproduce the issue

# Create point layer
pts = QgsVectorLayer("Point?crs=EPSG:3857", "pts", "memory")
for wkt in ('Point (0 0)', 'Point (0 2)'):
    geom = QgsGeometry.fromWkt(wkt)
    fet = QgsFeature()
    fet.setGeometry(geom)
    pts.dataProvider().addFeature(fet)


# Create polygon layer
poly = QgsVectorLayer("Polygon?crs=EPSG:3857", "polygon", "memory")
wkt = 'Polygon ((-1 1, 1 1, 1 -1, -1 -1, -1 1))'
geom = QgsGeometry.fromWkt(wkt)
fet = QgsFeature()
fet.setGeometry(geom)
poly.dataProvider().addFeature(fet)


# Keep point outside polygon using "Difference"
pts_diff = processing.run("native:difference", {
    'INPUT':pts,
    'OVERLAY':poly,
    'OUTPUT':'TEMPORARY_OUTPUT'
    })["OUTPUT"]


# Keep point outside polygon using "Extract by Location" with "disjoint" predicate
pts_extract = processing.run("native:extractbylocation", {
    'INPUT':pts,
    'PREDICATE':[2],
    'INTERSECT':poly,
    'OUTPUT':'TEMPORARY_OUTPUT'
    })["OUTPUT"]


# Check remaining points
pts_diff.dataProvider().featureCount()     # 2 (bad)
pts_extract.dataProvider().featureCount()  # 1 (ok)

QgsProject.instance().addMapLayer(poly)
QgsProject.instance().addMapLayer(pts_diff)  # show one point on map but two points in attribute table
QgsProject.instance().addMapLayer(pts_extract)  # show one point on map and in attribute table

Versions

QGIS 3.40.4
GEOS 3.13.0-CAPI-1.19.0

Windows 11 Enterprise 23H2

Supported QGIS version

• I'm running a supported QGIS version according to the roadmap.

New profile

• I tried with a new QGIS profile

Additional context

No response

[adansereau]

I came across the same problem (QGIS 3.40.1, GEOS 3.13.0-CAPI-1.19.0).

It also occurs in the test data at https://github.com/qgis/QGIS/blob/master/python/plugins/processing/tests/testdata/expected/difference_points_vs_polys.gml, where we can see all 9 features from the input are kept in the output layer, but there are ony 3 points on the map.

---

Based on the title of the PR ("Fix difference algorithm when differencing a point results in a null geometry"), and more specifically on #49300 (comment), I think this is caused by NULL geometries in the output.

These NULL geometries also appear in other functions that rely on the Difference, such as Union or Symmetrical difference.

For example, when applying these operations on a single layer containing a single feature, for Point, Line, and Polygon, where the layer is both the input and the overlay (ie let's say we have two layers that happen to have the same feature), if we look at the geometry in the output layer:

	Point (0 0)	LineString (0 0, 0 1)	Polygon ((0 0, 1 0, 1 1, 0 1, 0 0)
INTERSECTION (ok)	[<QgsGeometry: MultiPoint ((0 0))>]	[<QgsGeometry: MultiLineString ((0 0, 0 1))>]	[<QgsGeometry: MultiPolygon (((0 0, 0 1, 1 1, 1 0, 0 0)))>]
DIFFERENCE (not ok)	[<QgsGeometry: null>]	[]	[]
SYMMETRICAL DIFFERENCE (not ok)	[<QgsGeometry: null>, <QgsGeometry: null>]	[]	[]
UNION (not ok)	[<QgsGeometry: MultiPoint ((0 0))>, <QgsGeometry: null>, <QgsGeometry: null>]	[<QgsGeometry: MultiLineString ((0 0, 0 1))>]	[<QgsGeometry: MultiPolygon (((0 0, 0 1, 1 1, 1 0, 0 0)))>]

code example

# small function to shorten layer creation
def layerFromWKTList(vectorLayer, wkt_strings):
  for wkt in wkt_strings:
    geom = QgsGeometry.fromWkt(wkt)
    fet = QgsFeature()
    fet.setGeometry(geom)
    vectorLayer.dataProvider().addFeature(fet)
  return vectorLayer

# create the 3 layers to test

pt_layer = layerFromWKTList(
  QgsVectorLayer("Point?crs=EPSG:3857", "pt", "memory"),
  ('Point (0 0)',)
)

line_layer = layerFromWKTList(
  QgsVectorLayer("LineString?crs=EPSG:3857", "line", "memory"),
  ('LineString (0 0, 0 1)',)
)

poly_layer = layerFromWKTList(
  QgsVectorLayer("Polygon?crs=EPSG:3857", "poly", "memory"),
  ('Polygon ((0 0, 1 0, 1 1, 0 1, 0 0)',)
)

# which operations to try
# (intersection is just to help see where the union result comes from)
algos = ('intersection', 'difference', 'symmetricaldifference', 'union')

layers = (pt_layer, line_layer, poly_layer,)

# loop on the operations and the layers
for alg in algos:
  print()
  print(alg.upper())
  
  for layer in layers:
      # print geometry type
      print(QgsWkbTypes.displayString(layer.getGeometry(1).wkbType()), end=":")
      
      # apply operation on the layer
      result = processing.run(f"native:{alg}", {
          'INPUT':layer,
          'OVERLAY':layer,
          'OUTPUT':'TEMPORARY_OUTPUT'
      })["OUTPUT"]
      
      # print the geometry content
      print([f.geometry() for f in result.getFeatures()])


# printed result

# INTERSECTION (ok)
# Point:[<QgsGeometry: MultiPoint ((0 0))>]
# LineString:[<QgsGeometry: MultiLineString ((0 0, 0 1))>]
# Polygon:[<QgsGeometry: MultiPolygon (((0 0, 0 1, 1 1, 1 0, 0 0)))>]
#
# DIFFERENCE (not ok for Point)
# Point:[<QgsGeometry: null>]
# LineString:[]
# Polygon:[]
#
# SYMMETRICALDIFFERENCE (not ok for Point)
# Point:[<QgsGeometry: null>, <QgsGeometry: null>]
# LineString:[]
# Polygon:[]
#
# UNION (not ok for Point)
# Point:[<QgsGeometry: MultiPoint ((0 0))>, <QgsGeometry: null>, <QgsGeometry: null>]
# LineString:[<QgsGeometry: MultiLineString ((0 0, 0 1))>]
# Polygon:[<QgsGeometry: MultiPolygon (((0 0, 0 1, 1 1, 1 0, 0 0)))>]

---

Following the original reprex, we can see that the output layer contains a feature with <QgsGeometry: null>:

# geometry in input layer
print("pts:", [f.geometry() for f in pts.getFeatures()])
# pts: [<QgsGeometry: Point (0 0)>, <QgsGeometry: Point (0 2)>]
# --> two geometries in the input layer

# geometry in output layer
print("pts_diff:", [f.geometry() for f in pts_diff.getFeatures()])
# pts_diff: [<QgsGeometry: null>, <QgsGeometry: Point (0 2)>]
# --> the first feature has a null geometry
# (this would explain why it appears in the attribute table but not on the map)

print([f.geometry().isNull() for f in pts_diff.getFeatures()])  # [True, False]
print([f.geometry().isEmpty() for f in pts_diff.getFeatures()])  # [True, False]
print([f.hasGeometry() for f in pts_diff.getFeatures()])  # [False, True]
print([f.geometry().wkbType() for f in pts_diff.getFeatures()])  # [<WkbType.Unknown: 0>, <WkbType.Point: 1>]
print([f.geometry().type() for f in pts_diff.getFeatures()])  # [<GeometryType.Unknown: 3>, <GeometryType.Point: 0>]

# geometry in expected layer
print("pts_extract:", [f.geometry() for f in pts_extract.getFeatures()])
# pts_extract: [<QgsGeometry: Point (0 2)>]

A workaround I found is to apply the "Remove null geometries" tool on the output layer of "Difference" tool to get rid of the unwanted features.

pts_without_null = processing.run("native:removenullgeometries",  {
    'INPUT':pts_diff,
    'OUTPUT':'TEMPORARY_OUTPUT'
    })["OUTPUT"]

print("pts_without_null:", [f.geometry() for f in pts_without_null.getFeatures()])
# pts_without_null: [<QgsGeometry: Point (0 2)>]
# --> OK

[adansereau]

While I am neither familiar with QGIS development nor with C++ code, it seems to me that replacing the following line in QgsOverlayUtils::difference

QGIS/src/analysis/processing/qgsoverlayutils.cpp

Lines 195 to 196 in d8195bb

	if ( !geom.isNull() && !sanitizeDifferenceResult( geom, geometryType, flags ) )
	continue;

with

-      if ( !geom.isNull() && !sanitizeDifferenceResult( geom, geometryType, flags ) )
+      if ( geom.isNull() || !sanitizeDifferenceResult( geom, geometryType, flags ) )
        continue;

could solve the problem. (unless there are other side effects that I am not aware of)

Indeed, in both cases, only a non-null geometry will be passed to sanitizeDifferenceResult() (because we want to avoid the error GEOS geoprocessing error: difference failed thrown at https://github.com/qgis/QGIS/blob/d8195bb589f9eed95e35418d7ca31086f9d2dfba/src/analysis/processing/qgsoverlayutils.cpp#L64C1-L71C4)

• with !geom.isNull() &&, a null geometry is false for the left side of the condition and the right side will not be evaluated (because "AND" is already false)
• with geom.isNull() ||, a null geometry is true for the left side and the right side will not be evaluated (because "OR" is already true)

In the first case, only EMPTY geometry features can be skipped, while in the second case a NULL geometry feature can also reach the continue statement (and thus not be added to the output layer).