[geometry fixer] fix gap check geometry fix edge case

The case where the gap was containing a point that was lying inside (or on the edge) of a neighbouring polygon used to fail, as we were using that neighbour's points to snap to (since it was also at distance 0). This was leading to flaky wrong results.
qgis · Oct 15, 2020 · 279b051 · 279b051
1 parent 5ddf9aa
commit 279b051
Show file tree

Hide file tree

Showing 7 changed files with 78 additions and 5 deletions.
diff --git a/src/analysis/vector/geometry_checker/qgsgeometrygapcheck.cpp b/src/analysis/vector/geometry_checker/qgsgeometrygapcheck.cpp
@@ -367,19 +367,39 @@ bool QgsGeometryGapCheck::mergeWithNeighbor( const QMap<QString, QgsFeaturePool
     return false;
   }
 
-  QgsSpatialIndex neighbourIndex( QgsSpatialIndex::Flag::FlagStoreFeatureGeometries );
-  neighbourIndex.addFeatures( neighbours );
 
   QgsPolyline snappedRing;
   QgsVertexIterator iterator = errGeometry->vertices();
   while ( iterator.hasNext() )
   {
     QgsPoint pt = iterator.next();
     QgsVertexId id;
-    QgsGeometry closestGeom = neighbourIndex.geometry( neighbourIndex.nearestNeighbor( QgsPointXY( pt ) ).first() );
-    if ( !closestGeom.isEmpty() )
+
+    // We iterate through all neighbours points to find the closest one.
+    // Note that we cannot use a nearestNeighbour as it's not guaranteed to return the feature
+    // with the closest vertex (distance can be 0 if we are lying inside the feature, even if far
+    // from any vertex)
+    double nearestDist = -1;
+    QgsPoint closestPoint;
+    QListIterator<QgsFeature> neighboursIt( neighbours );
+    while ( neighboursIt.hasNext() )
+    {
+      QgsGeometry neighbourGeom = neighboursIt.next().geometry();
+      QgsVertexIterator neighbourIterator = neighbourGeom.vertices();
+      while ( neighbourIterator.hasNext() )
+      {
+        QgsPoint neighbourPt = neighbourIterator.next();
+        double dist = pt.distance( neighbourPt );
+        if ( nearestDist == -1 || dist < nearestDist )
+        {
+          nearestDist = dist;
+          closestPoint = neighbourPt;
+        }
+      }
+    }
+
+    if ( nearestDist != -1 )
     {
-      QgsPoint closestPoint = QgsGeometryUtils::closestVertex( *closestGeom.constGet(), pt, id );
       snappedRing.append( closestPoint );
     }
   }

diff --git a/tests/src/geometry_checker/testqgsgeometrychecks.cpp b/tests/src/geometry_checker/testqgsgeometrychecks.cpp
@@ -98,6 +98,7 @@ class TestQgsGeometryChecks: public QObject
     void testSelfContactCheck();
     void testSelfIntersectionCheck();
     void testSliverPolygonCheck();
+    void testGapCheckPointInPoly();
 };
 
 void TestQgsGeometryChecks::initTestCase()
@@ -1156,6 +1157,56 @@ void TestQgsGeometryChecks::testSliverPolygonCheck()
   cleanupTestContext( testContext );
 }
 
+void TestQgsGeometryChecks::testGapCheckPointInPoly()
+{
+  // The case where the gap was containing a point that was lying inside (or on the edge)
+  // of a neighbouring polygon used to fail, as we were using that neighbour's points to
+  // snap to (since it was also at distance 0). This was leading to flaky wrong results.
+
+  QTemporaryDir dir;
+  QMap<QString, QString> layers;
+  layers.insert( "gap_layer_point_in_poly.shp", "" );
+  auto testContext = createTestContext( dir, layers );
+
+  // Test detection
+  QList<QgsGeometryCheckError *> checkErrors;
+  QStringList messages;
+
+  QVariantMap configuration;
+
+  QgsProject::instance()->setCrs( QgsCoordinateReferenceSystem::fromEpsgId( 2056 ) );
+
+  QgsGeometryGapCheck check( testContext.first, configuration );
+  QgsFeedback feedback;
+  check.collectErrors( testContext.second, checkErrors, messages, &feedback );
+  listErrors( checkErrors, messages );
+
+  QCOMPARE( checkErrors.size(), 1 );
+
+  QgsGeometryCheckError *error = checkErrors.first();
+  QCOMPARE( error->contextBoundingBox().snappedToGrid( 100.0 ), QgsRectangle( 2.5372e+06, 1.1522e+06, 2.5375e+06, 1.1524e+06 ) );
+  QCOMPARE( error->affectedAreaBBox().snappedToGrid( 100.0 ), QgsRectangle( 2.5373e+06, 1.1523e+06, 2.5375e+06, 1.1523e+06 ) );
+
+  // Test fixes
+  QgsFeature f;
+  testContext.second[layers["gap_layer_point_in_poly.shp"]]->getFeature( 1, f );
+  double areaOld = f.geometry().area();
+  QCOMPARE( areaOld, 19913.135772452362 );
+
+  QgsGeometryCheck::Changes changes;
+  QMap<QString, int> mergeAttrs;
+  error->check()->fixError( testContext.second, error, QgsGeometryGapCheck::MergeLongestEdge, mergeAttrs, changes );
+
+  // Ensure it worked
+  QCOMPARE( error->status(), QgsGeometryCheckError::StatusFixed );
+
+  // Ensure it worked on geom
+  testContext.second[layers["gap_layer_point_in_poly.shp"]]->getFeature( 1, f );
+  QVERIFY( f.geometry().area() > areaOld );
+
+  cleanupTestContext( testContext );
+}
+
 ///////////////////////////////////////////////////////////////////////////////
 
 double TestQgsGeometryChecks::layerToMapUnits( const QgsMapLayer *layer, const QgsCoordinateReferenceSystem &mapCrs ) const

diff --git a/tests/testdata/geometry_checker/gap_layer_point_in_poly.cpg b/tests/testdata/geometry_checker/gap_layer_point_in_poly.cpg
@@ -0,0 +1 @@
+UTF-8
diff --git a/tests/testdata/geometry_checker/gap_layer_point_in_poly.dbf b/tests/testdata/geometry_checker/gap_layer_point_in_poly.dbf
diff --git a/tests/testdata/geometry_checker/gap_layer_point_in_poly.prj b/tests/testdata/geometry_checker/gap_layer_point_in_poly.prj
@@ -0,0 +1 @@
+PROJCS["CH1903+_LV95",GEOGCS["GCS_CH1903+",DATUM["D_CH1903+",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],PARAMETER["False_Easting",2600000.0],PARAMETER["False_Northing",1200000.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Azimuth",90.0],PARAMETER["Longitude_Of_Center",7.43958333333333],PARAMETER["Latitude_Of_Center",46.9524055555556],UNIT["Meter",1.0]]
diff --git a/tests/testdata/geometry_checker/gap_layer_point_in_poly.shp b/tests/testdata/geometry_checker/gap_layer_point_in_poly.shp
diff --git a/tests/testdata/geometry_checker/gap_layer_point_in_poly.shx b/tests/testdata/geometry_checker/gap_layer_point_in_poly.shx