Branch data     Line data    Source code

       1                 :            : /***************************************************************************
       2                 :            :                          qgsalgorithmcalculateoverlaps.cpp
       3                 :            :                          ------------------
       4                 :            :     begin                : May 2019
       5                 :            :     copyright            : (C) 2019 by Nyall Dawson
       6                 :            :     email                : nyall dot dawson at gmail dot com
       7                 :            :  ***************************************************************************/
       8                 :            : 
       9                 :            : /***************************************************************************
      10                 :            :  *                                                                         *
      11                 :            :  *   This program is free software; you can redistribute it and/or modify  *
      12                 :            :  *   it under the terms of the GNU General Public License as published by  *
      13                 :            :  *   the Free Software Foundation; either version 2 of the License, or     *
      14                 :            :  *   (at your option) any later version.                                   *
      15                 :            :  *                                                                         *
      16                 :            :  ***************************************************************************/
      17                 :            : 
      18                 :            : #include "qgsalgorithmcalculateoverlaps.h"
      19                 :            : #include "qgsvectorlayer.h"
      20                 :            : #include "qgsgeometryengine.h"
      21                 :            : 
      22                 :            : ///@cond PRIVATE
      23                 :            : 
      24                 :          0 : QString QgsCalculateVectorOverlapsAlgorithm::name() const
      25                 :            : {
      26                 :          0 :   return QStringLiteral( "calculatevectoroverlaps" );
      27                 :            : }
      28                 :            : 
      29                 :          0 : QString QgsCalculateVectorOverlapsAlgorithm::displayName() const
      30                 :            : {
      31                 :          0 :   return QObject::tr( "Overlap analysis" );
      32                 :            : }
      33                 :            : 
      34                 :          0 : QStringList QgsCalculateVectorOverlapsAlgorithm::tags() const
      35                 :            : {
      36                 :          0 :   return QObject::tr( "vector,overlay,area,percentage,intersection" ).split( ',' );
      37                 :          0 : }
      38                 :            : 
      39                 :          0 : QString QgsCalculateVectorOverlapsAlgorithm::group() const
      40                 :            : {
      41                 :          0 :   return QObject::tr( "Vector analysis" );
      42                 :            : }
      43                 :            : 
      44                 :          0 : QString QgsCalculateVectorOverlapsAlgorithm::groupId() const
      45                 :            : {
      46                 :          0 :   return QStringLiteral( "vectoranalysis" );
      47                 :            : }
      48                 :            : 
      49                 :          0 : void QgsCalculateVectorOverlapsAlgorithm::initAlgorithm( const QVariantMap & )
      50                 :            : {
      51                 :          0 :   addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList< int >() << QgsProcessing::TypeVectorPolygon ) );
      52                 :          0 :   addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "LAYERS" ), QObject::tr( "Overlay layers" ), QgsProcessing::TypeVectorPolygon ) );
      53                 :          0 :   addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Overlap" ) ) );
      54                 :          0 : }
      55                 :            : 
      56                 :          0 : QIcon QgsCalculateVectorOverlapsAlgorithm::icon() const
      57                 :            : {
      58                 :          0 :   return QgsApplication::getThemeIcon( QStringLiteral( "/algorithms/mAlgorithmClip.svg" ) );
      59                 :          0 : }
      60                 :            : 
      61                 :          0 : QString QgsCalculateVectorOverlapsAlgorithm::svgIconPath() const
      62                 :            : {
      63                 :          0 :   return QgsApplication::iconPath( QStringLiteral( "/algorithms/mAlgorithmClip.svg" ) );
      64                 :          0 : }
      65                 :            : 
      66                 :          0 : QString QgsCalculateVectorOverlapsAlgorithm::shortHelpString() const
      67                 :            : {
      68                 :          0 :   return QObject::tr( "This algorithm calculates the area and percentage cover by which features from an input layer "
      69                 :            :                       "are overlapped by features from a selection of overlay layers.\n\n"
      70                 :            :                       "New attributes are added to the output layer reporting the total area of overlap and percentage of the input feature overlapped "
      71                 :            :                       "by each of the selected overlay layers." );
      72                 :            : }
      73                 :            : 
      74                 :          0 : QgsCalculateVectorOverlapsAlgorithm *QgsCalculateVectorOverlapsAlgorithm::createInstance() const
      75                 :            : {
      76                 :          0 :   return new QgsCalculateVectorOverlapsAlgorithm();
      77                 :          0 : }
      78                 :            : 
      79                 :          0 : bool QgsCalculateVectorOverlapsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
      80                 :            : {
      81                 :          0 :   mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
      82                 :          0 :   if ( !mSource )
      83                 :          0 :     throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
      84                 :            : 
      85                 :          0 :   mOutputFields = mSource->fields();
      86                 :            : 
      87                 :          0 :   const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "LAYERS" ), context );
      88                 :          0 :   mOverlayerSources.reserve( layers.size() );
      89                 :          0 :   mLayerNames.reserve( layers.size() );
      90                 :          0 :   for ( QgsMapLayer *layer : layers )
      91                 :            :   {
      92                 :          0 :     if ( QgsVectorLayer *vl = qobject_cast< QgsVectorLayer * >( layer ) )
      93                 :            :     {
      94                 :          0 :       mLayerNames << layer->name();
      95                 :          0 :       mOverlayerSources.emplace_back( std::make_unique< QgsVectorLayerFeatureSource >( vl ) );
      96                 :          0 :       mOutputFields.append( QgsField( QStringLiteral( "%1_area" ).arg( vl->name() ), QVariant::Double ) );
      97                 :          0 :       mOutputFields.append( QgsField( QStringLiteral( "%1_pc" ).arg( vl->name() ), QVariant::Double ) );
      98                 :          0 :     }
      99                 :            :   }
     100                 :            : 
     101                 :          0 :   mOutputType = mSource->wkbType();
     102                 :          0 :   mCrs = mSource->sourceCrs();
     103                 :          0 :   mInputCount = mSource->featureCount();
     104                 :          0 :   mInputFeatures = mSource->getFeatures();
     105                 :            :   return true;
     106                 :          0 : }
     107                 :            : 
     108                 :          0 : QVariantMap QgsCalculateVectorOverlapsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
     109                 :            : {
     110                 :          0 :   QString destId;
     111                 :          0 :   std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, destId, mOutputFields,
     112                 :          0 :                                           mOutputType, mCrs ) );
     113                 :          0 :   if ( !sink )
     114                 :          0 :     throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
     115                 :            : 
     116                 :            :   // build a spatial index for each constraint layer for speed. We also store input constraint geometries here,
     117                 :            :   // to avoid refetching and projecting them later
     118                 :          0 :   QList< QgsSpatialIndex > spatialIndices;
     119                 :          0 :   spatialIndices.reserve( mLayerNames.size() );
     120                 :          0 :   auto nameIt = mLayerNames.constBegin();
     121                 :          0 :   for ( auto sourceIt = mOverlayerSources.begin(); sourceIt != mOverlayerSources.end(); ++sourceIt, ++nameIt )
     122                 :            :   {
     123                 :          0 :     feedback->pushInfo( QObject::tr( "Preparing %1" ).arg( *nameIt ) );
     124                 :          0 :     QgsFeatureIterator featureIt = ( *sourceIt )->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() ).setDestinationCrs( mCrs, context.transformContext() ).setInvalidGeometryCheck( context.invalidGeometryCheck() ).setInvalidGeometryCallback( context.invalidGeometryCallback() ) );
     125                 :          0 :     spatialIndices << QgsSpatialIndex( featureIt, feedback, QgsSpatialIndex::FlagStoreFeatureGeometries );
     126                 :          0 :   }
     127                 :            : 
     128                 :          0 :   QgsDistanceArea da;
     129                 :          0 :   da.setSourceCrs( mCrs, context.transformContext() );
     130                 :          0 :   da.setEllipsoid( context.ellipsoid() );
     131                 :            : 
     132                 :            :   // loop through input
     133                 :          0 :   double step = mInputCount > 0 ? 100.0 / mInputCount : 0;
     134                 :          0 :   long i = 0;
     135                 :          0 :   QgsFeature feature;
     136                 :          0 :   while ( mInputFeatures.nextFeature( feature ) )
     137                 :            :   {
     138                 :          0 :     if ( feedback->isCanceled() )
     139                 :          0 :       break;
     140                 :            : 
     141                 :          0 :     QgsAttributes outAttributes = feature.attributes();
     142                 :          0 :     if ( feature.hasGeometry() && !qgsDoubleNear( feature.geometry().area(), 0.0 ) )
     143                 :            :     {
     144                 :          0 :       const QgsGeometry inputGeom = feature.geometry();
     145                 :          0 :       const double inputArea = da.measureArea( inputGeom );
     146                 :            : 
     147                 :            :       // prepare for lots of intersection tests (for speed)
     148                 :          0 :       std::unique_ptr< QgsGeometryEngine > bufferGeomEngine( QgsGeometry::createGeometryEngine( inputGeom.constGet() ) );
     149                 :          0 :       bufferGeomEngine->prepareGeometry();
     150                 :            : 
     151                 :            :       // calculate overlap attributes
     152                 :          0 :       auto spatialIteratorIt = spatialIndices.begin();
     153                 :          0 :       for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it, ++spatialIteratorIt )
     154                 :            :       {
     155                 :          0 :         if ( feedback->isCanceled() )
     156                 :          0 :           break;
     157                 :            : 
     158                 :          0 :         QgsSpatialIndex &index = *spatialIteratorIt;
     159                 :          0 :         const QList<QgsFeatureId> matches = index.intersects( inputGeom.boundingBox() );
     160                 :          0 :         QVector< QgsGeometry > intersectingGeoms;
     161                 :          0 :         intersectingGeoms.reserve( matches.count() );
     162                 :          0 :         for ( QgsFeatureId match : matches )
     163                 :            :         {
     164                 :          0 :           if ( feedback->isCanceled() )
     165                 :          0 :             break;
     166                 :            : 
     167                 :          0 :           const QgsGeometry overlayGeometry = index.geometry( match );
     168                 :          0 :           if ( bufferGeomEngine->intersects( overlayGeometry.constGet() ) )
     169                 :            :           {
     170                 :          0 :             intersectingGeoms.append( overlayGeometry );
     171                 :          0 :           }
     172                 :          0 :         }
     173                 :            : 
     174                 :          0 :         if ( feedback->isCanceled() )
     175                 :          0 :           break;
     176                 :            : 
     177                 :            :         // dissolve intersecting features, calculate total area of them within our buffer
     178                 :          0 :         QgsGeometry overlayDissolved = QgsGeometry::unaryUnion( intersectingGeoms );
     179                 :            : 
     180                 :          0 :         if ( feedback->isCanceled() )
     181                 :          0 :           break;
     182                 :            : 
     183                 :          0 :         QgsGeometry overlayIntersection = inputGeom.intersection( overlayDissolved );
     184                 :            : 
     185                 :          0 :         const double overlayArea = da.measureArea( overlayIntersection );
     186                 :          0 :         outAttributes.append( overlayArea );
     187                 :          0 :         outAttributes.append( 100 * overlayArea / inputArea );
     188                 :          0 :       }
     189                 :          0 :     }
     190                 :            :     else
     191                 :            :     {
     192                 :            :       // input feature has no geometry
     193                 :          0 :       for ( auto it = mOverlayerSources.begin(); it != mOverlayerSources.end(); ++ it )
     194                 :            :       {
     195                 :          0 :         outAttributes.append( QVariant() );
     196                 :          0 :         outAttributes.append( QVariant() );
     197                 :          0 :       }
     198                 :            :     }
     199                 :            : 
     200                 :          0 :     feature.setAttributes( outAttributes );
     201                 :          0 :     sink->addFeature( feature, QgsFeatureSink::FastInsert );
     202                 :            : 
     203                 :          0 :     i++;
     204                 :          0 :     feedback->setProgress( i * step );
     205                 :          0 :   }
     206                 :            : 
     207                 :          0 :   QVariantMap outputs;
     208                 :          0 :   outputs.insert( QStringLiteral( "OUTPUT" ), destId );
     209                 :          0 :   return outputs;
     210                 :          0 : }
     211                 :            : 
     212                 :            : ///@endcond