Skip to content
Permalink
Browse files

Merge pull request #39392 from olivierdalang/geometry_fixer_fixes

Geometry fixer fixes
  • Loading branch information
m-kuhn committed Nov 8, 2020
2 parents 7d012ee + 4d11d79 commit 63d875b084551c5c4309ca985ecabcf06b77fd60
@@ -367,20 +367,33 @@ bool QgsGeometryGapCheck::mergeWithNeighbor( const QMap<QString, QgsFeaturePool
return false;
}

QgsSpatialIndex neighbourIndex( QgsSpatialIndex::Flag::FlagStoreFeatureGeometries );
neighbourIndex.addFeatures( neighbours );
// Create an index of all neighbouring vertices
QgsSpatialIndex neighbourVerticesIndex( QgsSpatialIndex::Flag::FlagStoreFeatureGeometries );
int id = 0;
for ( const QgsFeature &neighbour : neighbours )
{
QgsVertexIterator vit = neighbour.geometry().vertices();
while ( vit.hasNext() )
{
QgsPoint pt = vit.next();
QgsFeature f;
f.setId( id ); // required for SpatialIndex to return the correct result
f.setGeometry( QgsGeometry( pt.clone() ) );
neighbourVerticesIndex.addFeature( f );
id++;
}
}

// Snap to the closest vertex
QgsPolyline snappedRing;
QgsVertexIterator iterator = errGeometry->vertices();
while ( iterator.hasNext() )
{
QgsPoint pt = iterator.next();
QgsVertexId id;
QgsGeometry closestGeom = neighbourIndex.geometry( neighbourIndex.nearestNeighbor( QgsPointXY( pt ) ).first() );
QgsGeometry closestGeom = neighbourVerticesIndex.geometry( neighbourVerticesIndex.nearestNeighbor( QgsPointXY( pt ) ).first() );
if ( !closestGeom.isEmpty() )
{
QgsPoint closestPoint = QgsGeometryUtils::closestVertex( *closestGeom.constGet(), pt, id );
snappedRing.append( closestPoint );
snappedRing.append( QgsPoint( closestGeom.vertexAt( 0 ) ) );
}
}

@@ -112,7 +112,7 @@ void QgsGeometryOverlapCheck::fixError( const QMap<QString, QgsFeaturePool *> &f
QgsGeometryCheckerUtils::LayerFeature layerFeatureA( featurePoolA, featureA, mContext, true );
QgsGeometryCheckerUtils::LayerFeature layerFeatureB( featurePoolB, featureB, mContext, true );
const QgsGeometry geometryA = layerFeatureA.geometry();
std::unique_ptr< QgsGeometryEngine > geomEngineA = QgsGeometryCheckerUtils::createGeomEngine( geometryA.constGet(), mContext->reducedTolerance );
std::unique_ptr< QgsGeometryEngine > geomEngineA = QgsGeometryCheckerUtils::createGeomEngine( geometryA.constGet(), mContext->tolerance );
geomEngineA->prepareGeometry();

const QgsGeometry geometryB = layerFeatureB.geometry();
@@ -153,7 +153,8 @@ void QgsGeometryOverlapCheck::fixError( const QMap<QString, QgsFeaturePool *> &f
}
else if ( method == Subtract )
{
std::unique_ptr< QgsAbstractGeometry > diff1( geomEngineA->difference( interPart, &errMsg ) );
std::unique_ptr< QgsGeometryEngine > geomEngineDiffA = QgsGeometryCheckerUtils::createGeomEngine( geometryA.constGet(), 0 );
std::unique_ptr< QgsAbstractGeometry > diff1( geomEngineDiffA->difference( interPart, &errMsg ) );
if ( !diff1 || diff1->isEmpty() )
{
diff1.reset();
@@ -162,9 +163,8 @@ void QgsGeometryOverlapCheck::fixError( const QMap<QString, QgsFeaturePool *> &f
{
QgsGeometryCheckerUtils::filter1DTypes( diff1.get() );
}
std::unique_ptr< QgsGeometryEngine > geomEngineB = QgsGeometryCheckerUtils::createGeomEngine( geometryB.constGet(), mContext->reducedTolerance );
geomEngineB->prepareGeometry();
std::unique_ptr< QgsAbstractGeometry > diff2( geomEngineB->difference( interPart, &errMsg ) );
std::unique_ptr< QgsGeometryEngine > geomEngineDiffB = QgsGeometryCheckerUtils::createGeomEngine( geometryB.constGet(), 0 );
std::unique_ptr< QgsAbstractGeometry > diff2( geomEngineDiffB->difference( interPart, &errMsg ) );
if ( !diff2 || diff2->isEmpty() )
{
diff2.reset();
@@ -98,6 +98,8 @@ class TestQgsGeometryChecks: public QObject
void testSelfContactCheck();
void testSelfIntersectionCheck();
void testSliverPolygonCheck();
void testGapCheckPointInPoly();
void testOverlapCheckToleranceBug();
};

void TestQgsGeometryChecks::initTestCase()
@@ -1156,6 +1158,116 @@ 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 );
}

void TestQgsGeometryChecks::testOverlapCheckToleranceBug()
{
// The overlap (intersection) was computed with a different tolerance when collecting errors
// than when fixing them, leading to failures to fix the issue esp. with big coordinates.
//
// Also, it used to offset unaffected points (far from the actual overlap) on the affected
// feature, leading to both unwanted shifts and remaining slivers.

QTemporaryDir dir;
QMap<QString, QString> layers;
layers.insert( "overlap_layer_tolerance_bug.shp", "" );
auto testContext = createTestContext( dir, layers );

// Test detection
QList<QgsGeometryCheckError *> checkErrors;
QStringList messages;

QVariantMap configuration;
configuration.insert( "gapThreshold", 1000.0 );

QgsProject::instance()->setCrs( QgsCoordinateReferenceSystem::fromEpsgId( 2056 ) );

QgsGeometryOverlapCheck check( testContext.first, configuration );
QgsFeedback feedback;
check.collectErrors( testContext.second, checkErrors, messages, &feedback );
listErrors( checkErrors, messages );

QCOMPARE( checkErrors.size(), 1 );

QgsGeometryCheckError *error = checkErrors.first();

// Test fixes
QgsFeature f;
testContext.second[layers["overlap_layer_tolerance_bug.shp"]]->getFeature( 0, f );
double areaOld = f.geometry().area();
QgsPoint pointOld_1 = f.geometry().vertexAt( 1 );
QgsPoint pointOld_2 = f.geometry().vertexAt( 2 );

// Just making sure we've got the right feature/point
QCOMPARE( areaOld, 10442.710061549426 );
QGSCOMPARENEARPOINT( pointOld_1, QgsPoint( 2537221.53079314017668366, 1152360.02460834058001637 ), 0.00001 );
QGSCOMPARENEARPOINT( pointOld_2, QgsPoint( 2537366.84566075634211302, 1152360.28978145681321621 ), 0.00001 );

QgsGeometryCheck::Changes changes;
QMap<QString, int> mergeAttrs;
error->check()->fixError( testContext.second, error, QgsGeometryOverlapCheck::Subtract, mergeAttrs, changes );

// Ensure it worked
QCOMPARE( error->status(), QgsGeometryCheckError::StatusFixed );

// Ensure it actually worked
testContext.second[layers["overlap_layer_tolerance_bug.shp"]]->getFeature( 0, f );
QVERIFY( f.geometry().area() < areaOld );
// And that we don't have unexpected changes on unaffected points
QCOMPARE( f.geometry().vertexAt( 1 ), pointOld_1 );
QCOMPARE( f.geometry().vertexAt( 2 ), pointOld_2 );

cleanupTestContext( testContext );
}

///////////////////////////////////////////////////////////////////////////////

double TestQgsGeometryChecks::layerToMapUnits( const QgsMapLayer *layer, const QgsCoordinateReferenceSystem &mapCrs ) const
@@ -0,0 +1 @@
UTF-8
Binary file not shown.
@@ -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]]
Binary file not shown.
Binary file not shown.
@@ -0,0 +1 @@
UTF-8
Binary file not shown.
@@ -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]]
Binary file not shown.
Binary file not shown.

0 comments on commit 63d875b

Please sign in to comment.
You can’t perform that action at this time.