Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
[FEATURE] Union algorithm for single layer
Resolves all overlapping geometries just like GRASS or Arc do.

So now we have two variants of union:
- union(A) - does union within geometries of one layer
- union(A,B) - does union between geometries of two layers

For union(A,B) algorithm if there are overlaps among geometries of layer A or among geometries of layer B,
these are not resolved: one needs to do union(union(A,B)) to resolve all overlaps, i.e. run single layer
union(X) on the produced result X=union(A,B)

This should also address issues raised in #17131
  • Loading branch information
wonder-sk committed May 10, 2018
1 parent c738bcf commit 64b8c72
Show file tree
Hide file tree
Showing 3 changed files with 262 additions and 31 deletions.
15 changes: 10 additions & 5 deletions src/analysis/processing/qgsalgorithmunion.cpp
Expand Up @@ -53,25 +53,23 @@ QgsProcessingAlgorithm *QgsUnionAlgorithm::createInstance() const
void QgsUnionAlgorithm::initAlgorithm( const QVariantMap & )
{
addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ) ) );
addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "OVERLAY" ), QObject::tr( "Union layer" ) ) );
addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "OVERLAY" ), QObject::tr( "Union layer" ), QList< int >(), QVariant(), true ) );

addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Union" ) ) );
}


QVariantMap QgsUnionAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
std::unique_ptr< QgsFeatureSource > sourceA( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
if ( !sourceA )
throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );

std::unique_ptr< QgsFeatureSource > sourceB( parameterAsSource( parameters, QStringLiteral( "OVERLAY" ), context ) );
if ( !sourceB )
throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "OVERLAY" ) ) );
// sourceB is optional so we do not throw an error if it is not a valid source

QgsWkbTypes::Type geomType = QgsWkbTypes::multiType( sourceA->wkbType() );

QgsFields fields = QgsProcessingUtils::combineFields( sourceA->fields(), sourceB->fields() );
QgsFields fields = sourceB ? QgsProcessingUtils::combineFields( sourceA->fields(), sourceB->fields() ) : sourceA->fields();

QString dest;
std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, geomType, sourceA->sourceCrs() ) );
Expand All @@ -81,6 +79,13 @@ QVariantMap QgsUnionAlgorithm::processAlgorithm( const QVariantMap &parameters,
QVariantMap outputs;
outputs.insert( QStringLiteral( "OUTPUT" ), dest );

if ( !sourceB )
{
// we are doing single layer union
QgsOverlayUtils::resolveOverlaps( *sourceA.get(), *sink.get(), feedback );
return outputs;
}

QList<int> fieldIndicesA = QgsProcessingUtils::fieldNamesToIndices( QStringList(), sourceA->fields() );
QList<int> fieldIndicesB = QgsProcessingUtils::fieldNamesToIndices( QStringList(), sourceB->fields() );

Expand Down
268 changes: 242 additions & 26 deletions src/analysis/processing/qgsoverlayutils.cpp
Expand Up @@ -20,6 +20,61 @@

///@cond PRIVATE

//! Makes sure that what came out from intersection of two geometries is good to be used in the output
static bool sanitizeIntersectionResult( QgsGeometry &geom, QgsWkbTypes::GeometryType geometryType )
{
if ( geom.isNull() )
{
// TODO: not sure if this ever happens - if it does, that means GEOS failed badly - would be good to have a test for such situation
throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: intersection failed." ), geom.lastError() ) );
}

// Intersection of geometries may give use also geometries we do not want in our results.
// For example, two square polygons touching at the corner have a point as the intersection, but no area.
// In other cases we may get a mixture of geometries in the output - we want to keep only the expected types.
if ( QgsWkbTypes::flatType( geom.wkbType() ) == QgsWkbTypes::GeometryCollection )
{
// try to filter out irrelevant parts with different geometry type than what we want
geom.convertGeometryCollectionToSubclass( geometryType );
if ( geom.isEmpty() )
return false;
}

if ( QgsWkbTypes::geometryType( geom.wkbType() ) != geometryType )
{
// we can't make use of this resulting geometry
return false;
}

// some data providers are picky about the geometries we pass to them: we can't add single-part geometries
// when we promised multi-part geometries, so ensure we have the right type
geom.convertToMultiType();

return true;
}


//! Makes sure that what came out from difference of two geometries is good to be used in the output
static bool sanitizeDifferenceResult( QgsGeometry &geom )
{
if ( geom.isNull() )
{
// TODO: not sure if this ever happens - if it does, that means GEOS failed badly - would be good to have a test for such situation
throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: difference failed." ), geom.lastError() ) );
}

// if geomB covers the whole source geometry, we get an empty geometry collection
if ( geom.isEmpty() )
return false;

// some data providers are picky about the geometries we pass to them: we can't add single-part geometries
// when we promised multi-part geometries, so ensure we have the right type
geom.convertToMultiType();

return true;
}


void QgsOverlayUtils::difference( const QgsFeatureSource &sourceA, const QgsFeatureSource &sourceB, QgsFeatureSink &sink, QgsProcessingContext &context, QgsProcessingFeedback *feedback, int &count, int totalCount, QgsOverlayUtils::DifferenceOutput outputAttrs )
{
QgsFeatureRequest requestB;
Expand Down Expand Up @@ -83,8 +138,7 @@ void QgsOverlayUtils::difference( const QgsFeatureSource &sourceA, const QgsFeat
geom = geom.difference( geomB );
}

// if geomB covers the whole source geometry, we get an empty geometry collection
if ( geom.isEmpty() )
if ( !sanitizeDifferenceResult( geom ) )
continue;

const QgsAttributes attrsA( featA.attributes() );
Expand All @@ -104,7 +158,6 @@ void QgsOverlayUtils::difference( const QgsFeatureSource &sourceA, const QgsFeat
}

QgsFeature outFeat;
geom.convertToMultiType();
outFeat.setGeometry( geom );
outFeat.setAttributes( attrs );
sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
Expand Down Expand Up @@ -179,35 +232,13 @@ void QgsOverlayUtils::intersection( const QgsFeatureSource &sourceA, const QgsFe
continue;

QgsGeometry intGeom = geom.intersection( tmpGeom );

if ( intGeom.isNull() )
{
// TODO: not sure if this ever happens - if it does, that means GEOS failed badly - would be good to have a test for such situation
throw QgsProcessingException( QStringLiteral( "%1\n\n%2" ).arg( QObject::tr( "GEOS geoprocessing error: intersection failed." ), intGeom.lastError() ) );
}

// Intersection of geometries may give use also geometries we do not want in our results.
// For example, two square polygons touching at the corner have a point as the intersection, but no area.
// In other cases we may get a mixture of geometries in the output - we want to keep only the expected types.
if ( QgsWkbTypes::flatType( intGeom.wkbType() ) == QgsWkbTypes::GeometryCollection )
{
// try to filter out irrelevant parts with different geometry type than what we want
intGeom.convertGeometryCollectionToSubclass( geometryType );
if ( intGeom.isEmpty() )
continue;
}

if ( QgsWkbTypes::geometryType( intGeom.wkbType() ) != geometryType )
{
// we can't make use of this resulting geometry
if ( !sanitizeIntersectionResult( intGeom, geometryType ) )
continue;
}

const QgsAttributes attrsB( featB.attributes() );
for ( int i = 0; i < fieldIndicesB.count(); ++i )
outAttributes[fieldIndicesA.count() + i] = attrsB[fieldIndicesB[i]];

intGeom.convertToMultiType();
outFeat.setGeometry( intGeom );
outFeat.setAttributes( outAttributes );
sink.addFeature( outFeat, QgsFeatureSink::FastInsert );
Expand All @@ -218,4 +249,189 @@ void QgsOverlayUtils::intersection( const QgsFeatureSource &sourceA, const QgsFe
}
}

void QgsOverlayUtils::resolveOverlaps( const QgsFeatureSource &source, QgsFeatureSink &sink, QgsProcessingFeedback *feedback )
{
int count = 0;
int totalCount = source.featureCount();
if ( totalCount == 0 )
return; // nothing to do here

QgsFeatureId newFid = -1;

QgsWkbTypes::GeometryType geometryType = QgsWkbTypes::geometryType( QgsWkbTypes::multiType( source.wkbType() ) );

QgsFeatureRequest requestOnlyGeoms;
requestOnlyGeoms.setSubsetOfAttributes( QgsAttributeList() );

QgsFeatureRequest requestOnlyAttrs;
requestOnlyAttrs.setFlags( QgsFeatureRequest::NoGeometry );

QgsFeatureRequest requestOnlyIds;
requestOnlyIds.setFlags( QgsFeatureRequest::NoGeometry );
requestOnlyIds.setSubsetOfAttributes( QgsAttributeList() );

// make a set of used feature IDs so they we do not try to reuse them for newly added features
QgsFeature f;
QSet<QgsFeatureId> fids;
QgsFeatureIterator it = source.getFeatures( requestOnlyIds );
while ( it.nextFeature( f ) )
{
if ( feedback->isCanceled() )
return;

fids.insert( f.id() );
}

QHash<QgsFeatureId, QgsGeometry> geometries;
QgsSpatialIndex index;
QHash<QgsFeatureId, QList<QgsFeatureId> > intersectingIds; // which features overlap a particular area

// resolve intersections

it = source.getFeatures( requestOnlyGeoms );
while ( it.nextFeature( f ) )
{
if ( feedback->isCanceled() )
return;

QgsFeatureId fid1 = f.id();
QgsGeometry g1 = f.geometry();

geometries.insert( fid1, g1 );
index.insertFeature( f );

QgsRectangle bbox( f.geometry().boundingBox() );
const QList<QgsFeatureId> ids = index.intersects( bbox );
for ( QgsFeatureId fid2 : ids )
{
if ( fid1 == fid2 )
continue;

QgsGeometry g2 = geometries.value( fid2 );
if ( !g1.intersects( g2 ) )
continue;

QgsGeometry geomIntersection = g1.intersection( g2 );
if ( !sanitizeIntersectionResult( geomIntersection, geometryType ) )
continue;

//
// add intersection geometry
//

// figure out new fid
while ( fids.contains( newFid ) )
--newFid;
fids.insert( newFid );

geometries.insert( newFid, geomIntersection );
QgsFeature fx( newFid );
fx.setGeometry( geomIntersection );

index.insertFeature( fx );

// figure out which feature IDs belong to this intersection. Some of the IDs can be of the newly
// created geometries - in such case we need to retrieve original IDs
QList<QgsFeatureId> lst;
if ( intersectingIds.contains( fid1 ) )
lst << intersectingIds.value( fid1 );
else
lst << fid1;
if ( intersectingIds.contains( fid2 ) )
lst << intersectingIds.value( fid2 );
else
lst << fid2;
intersectingIds.insert( newFid, lst );

//
// update f1
//

QgsGeometry g12 = g1.difference( g2 );

index.deleteFeature( f );
geometries.remove( fid1 );

if ( sanitizeDifferenceResult( g12 ) )
{
geometries.insert( fid1, g12 );

QgsFeature f1x( fid1 );
f1x.setGeometry( g12 );
index.insertFeature( f1x );
}

//
// update f2
//

QgsGeometry g21 = g2.difference( g1 );

QgsFeature f2old( fid2 );
f2old.setGeometry( g2 );
index.deleteFeature( f2old );

geometries.remove( fid2 );

if ( sanitizeDifferenceResult( g21 ) )
{
geometries.insert( fid2, g21 );

QgsFeature f2x( fid2 );
f2x.setGeometry( g21 );
index.insertFeature( f2x );
}

// update our temporary copy of the geometry to what is left from it
g1 = g12;
}

++count;
feedback->setProgress( count / ( double ) totalCount * 100. );
}

// release some memory of structures we don't need anymore

fids.clear();
index = QgsSpatialIndex();

// load attributes

QHash<QgsFeatureId, QgsAttributes> attributesHash;
it = source.getFeatures( requestOnlyAttrs );
while ( it.nextFeature( f ) )
{
if ( feedback->isCanceled() )
return;

attributesHash.insert( f.id(), f.attributes() );
}

// store stuff in the sink

for ( auto i = geometries.constBegin(); i != geometries.constEnd(); ++i )
{
if ( feedback->isCanceled() )
return;

QgsFeature outFeature( i.key() );
outFeature.setGeometry( i.value() );

if ( intersectingIds.contains( i.key() ) )
{
const QList<QgsFeatureId> ids = intersectingIds.value( i.key() );
for ( QgsFeatureId id : ids )
{
outFeature.setAttributes( attributesHash.value( id ) );
sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
}
}
else
{
outFeature.setAttributes( attributesHash.value( i.key() ) );
sink.addFeature( outFeature, QgsFeatureSink::FastInsert );
}
}
}

///@endcond PRIVATE
10 changes: 10 additions & 0 deletions src/analysis/processing/qgsoverlayutils.h
Expand Up @@ -43,6 +43,16 @@ namespace QgsOverlayUtils

void intersection( const QgsFeatureSource &sourceA, const QgsFeatureSource &sourceB, QgsFeatureSink &sink, QgsProcessingContext &context, QgsProcessingFeedback *feedback, int &count, int totalCount, const QList<int> &fieldIndicesA, const QList<int> &fieldIndicesB );

/**
* Copies features from the source to the sink and resolves overlaps: for each pair of overlapping features A and B
* it will produce:
* 1. a feature with geometry A - B with A's attributes
* 2. a feature with geometry B - A with B's attributes
* 3. two features with geometry intersection(A, B) - one with A's attributes, one with B's attributes.
*
* As a result, for all pairs of features in the output, a pair either havs no common interior or their interior is the same.
*/
void resolveOverlaps( const QgsFeatureSource &source, QgsFeatureSink &sink, QgsProcessingFeedback *feedback );
}

///@endcond PRIVATE
Expand Down

0 comments on commit 64b8c72

Please sign in to comment.