Skip to content
Permalink
Browse files

Move internal intersection to overlay utils, code de-duplication

  • Loading branch information
wonder-sk committed Apr 25, 2018
1 parent d0130d2 commit 85b2efae2c6a05e9859097a52d44e35c6c8b49b4
@@ -17,6 +17,7 @@

#include "qgsgeometrycollection.h"
#include "qgsgeometryengine.h"
#include "qgsoverlayutils.h"

///@cond PRIVATE

@@ -87,147 +88,23 @@ QVariantMap QgsIntersectionAlgorithm::processAlgorithm( const QVariantMap &param
const QStringList fieldsA = parameterAsFields( parameters, QStringLiteral( "INPUT_FIELDS" ), context );
const QStringList fieldsB = parameterAsFields( parameters, QStringLiteral( "OVERLAY_FIELDS" ), context );

QgsFields fieldListA;
QList<int> fieldIndicesA;
if ( !fieldsA.isEmpty() )
{
for ( const QString &f : fieldsA )
{
int idxA = sourceA->fields().lookupField( f );
if ( idxA >= 0 )
{
fieldIndicesA.append( idxA );
fieldListA.append( sourceA->fields()[idxA] );
}
}
}
else
{
fieldListA = sourceA->fields();
for ( int i = 0; i < fieldListA.count(); ++i )
fieldIndicesA.append( i );
}

QgsFields fieldListB;
QList<int> fieldIndicesB;
if ( !fieldsB.isEmpty() )
{
for ( const QString &f : fieldsB )
{
int idxB = sourceB->fields().lookupField( f );
if ( idxB >= 0 )
{
fieldIndicesB.append( idxB );
fieldListB.append( sourceB->fields()[idxB] );
}
}
}
else
{
fieldListB = sourceB->fields();
for ( int i = 0; i < fieldListB.count(); ++i )
fieldIndicesB.append( i );
}


QgsFields outputFields = QgsProcessingUtils::combineFields( fieldListA, fieldListB );
QList<int> fieldIndicesA = QgsOverlayUtils::fieldNamesToIndices( fieldsA, sourceA->fields() );
QList<int> fieldIndicesB = QgsOverlayUtils::fieldNamesToIndices( fieldsB, sourceB->fields() );

QgsFields outputFields = QgsProcessingUtils::combineFields(
QgsOverlayUtils::indicesToFields( fieldIndicesA, sourceA->fields() ),
QgsOverlayUtils::indicesToFields( fieldIndicesB, sourceB->fields() ) );

QString dest;
std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, outputFields, geomType, sourceA->sourceCrs() ) );

QVariantMap outputs;
outputs.insert( QStringLiteral( "OUTPUT" ), dest );

QgsFeatureRequest request;
request.setSubsetOfAttributes( QgsAttributeList() );
request.setDestinationCrs( sourceA->sourceCrs(), context.transformContext() );

QgsFeature outFeat;
QgsSpatialIndex indexB( sourceB->getFeatures( request ), feedback );

double total = 100.0 / ( sourceA->featureCount() ? sourceA->featureCount() : 1 );
int count = 0;
int total = sourceA->featureCount();

QgsFeature featA;
QgsFeatureIterator fitA = sourceA->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( fieldIndicesA ) );
while ( fitA.nextFeature( featA ) )
{
if ( feedback->isCanceled() )
break;

if ( !featA.hasGeometry() )
continue;

QgsGeometry geom( featA.geometry() );
QgsFeatureIds intersects = indexB.intersects( geom.boundingBox() ).toSet();

QgsFeatureRequest request;
request.setFilterFids( intersects );
request.setDestinationCrs( sourceA->sourceCrs(), context.transformContext() );
request.setSubsetOfAttributes( fieldIndicesB );

std::unique_ptr< QgsGeometryEngine > engine;
if ( !intersects.isEmpty() )
{
// use prepared geometries for faster intersection tests
engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
engine->prepareGeometry();
}

QgsAttributes outAttributes( outputFields.count() );
const QgsAttributes attrsA( featA.attributes() );
for ( int i = 0; i < fieldIndicesA.count(); ++i )
outAttributes[i] = attrsA[fieldIndicesA[i]];

QgsFeature featB;
QgsFeatureIterator fitB = sourceB->getFeatures( request );
while ( fitB.nextFeature( featB ) )
{
if ( feedback->isCanceled() )
break;

QgsGeometry tmpGeom( featB.geometry() );
if ( !engine->intersects( tmpGeom.constGet() ) )
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( QgsWkbTypes::geometryType( geomType ) );
if ( intGeom.isEmpty() )
continue;
}

if ( QgsWkbTypes::geometryType( intGeom.wkbType() ) != QgsWkbTypes::geometryType( geomType ) )
{
// we can't make use of this resulting geometry
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 );
}

++count;
feedback->setProgress( int( count * total ) );
}
QgsOverlayUtils::intersection( *sourceA.get(), *sourceB.get(), *sink.get(), context, feedback, count, total, fieldIndicesA, fieldIndicesB );

return outputs;
}
@@ -120,4 +120,131 @@ void QgsOverlayUtils::difference( const QgsFeatureSource &sourceA, const QgsFeat
}
}


void QgsOverlayUtils::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 )
{
QgsWkbTypes::GeometryType geometryType = QgsWkbTypes::geometryType( QgsWkbTypes::multiType( sourceA.wkbType() ) );
int attrCount = fieldIndicesA.count() + fieldIndicesB.count();

QgsFeatureRequest request;
request.setSubsetOfAttributes( QgsAttributeList() );
request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );

QgsFeature outFeat;
QgsSpatialIndex indexB( sourceB.getFeatures( request ), feedback );

if ( totalCount == 0 )
totalCount = 1; // avoid division by zero

QgsFeature featA;
QgsFeatureIterator fitA = sourceA.getFeatures( QgsFeatureRequest().setSubsetOfAttributes( fieldIndicesA ) );
while ( fitA.nextFeature( featA ) )
{
if ( feedback->isCanceled() )
break;

if ( !featA.hasGeometry() )
continue;

QgsGeometry geom( featA.geometry() );
QgsFeatureIds intersects = indexB.intersects( geom.boundingBox() ).toSet();

QgsFeatureRequest request;
request.setFilterFids( intersects );
request.setDestinationCrs( sourceA.sourceCrs(), context.transformContext() );
request.setSubsetOfAttributes( fieldIndicesB );

std::unique_ptr< QgsGeometryEngine > engine;
if ( !intersects.isEmpty() )
{
// use prepared geometries for faster intersection tests
engine.reset( QgsGeometry::createGeometryEngine( geom.constGet() ) );
engine->prepareGeometry();
}

QgsAttributes outAttributes( attrCount );
const QgsAttributes attrsA( featA.attributes() );
for ( int i = 0; i < fieldIndicesA.count(); ++i )
outAttributes[i] = attrsA[fieldIndicesA[i]];

QgsFeature featB;
QgsFeatureIterator fitB = sourceB.getFeatures( request );
while ( fitB.nextFeature( featB ) )
{
if ( feedback->isCanceled() )
break;

QgsGeometry tmpGeom( featB.geometry() );
if ( !engine->intersects( tmpGeom.constGet() ) )
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
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 );
}

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


QList<int> QgsOverlayUtils::fieldNamesToIndices( const QStringList &fieldNames, const QgsFields &fields )
{
QList<int> indices;
if ( !fieldNames.isEmpty() )
{
for ( const QString &f : fieldNames )
{
int idx = fields.lookupField( f );
if ( idx >= 0 )
indices.append( idx );
}
}
else
{
for ( int i = 0; i < fields.count(); ++i )
indices.append( i );
}
return indices;
}

QgsFields QgsOverlayUtils::indicesToFields( const QList<int> &indices, const QgsFields &fields )
{
QgsFields fieldsSubset;
for ( int i : indices )
fieldsSubset.append( fields.at( i ) );
return fieldsSubset;
}

///@endcond PRIVATE
@@ -16,12 +16,15 @@
#ifndef QGSOVERLAYUTILS_H
#define QGSOVERLAYUTILS_H

#include <QList>

#define SIP_NO_FILE

///@cond PRIVATE

class QgsFeatureSource;
class QgsFeatureSink;
class QgsFields;
class QgsProcessingContext;
class QgsProcessingFeedback;

@@ -38,6 +41,11 @@ namespace QgsOverlayUtils

void difference( const QgsFeatureSource &sourceA, const QgsFeatureSource &sourceB, QgsFeatureSink &sink, QgsProcessingContext &context, QgsProcessingFeedback *feedback, int &count, int totalCount, DifferenceOutput outputAttrs );

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 );

QList<int> fieldNamesToIndices( const QStringList &fieldNames, const QgsFields &fields );

QgsFields indicesToFields( const QList<int> &indices, const QgsFields &fields );
}

///@endcond PRIVATE

0 comments on commit 85b2efa

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