Branch data     Line data    Source code

       1                 :            : /***************************************************************************
       2                 :            :                          qgsalgorithmrastersampling.cpp
       3                 :            :                          --------------------------
       4                 :            :     begin                : August 2020
       5                 :            :     copyright            : (C) 2020 by Mathieu Pellerin
       6                 :            :     email                : nirvn dot asia 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 "qgsalgorithmrastersampling.h"
      19                 :            : #include "qgsmultipoint.h"
      20                 :            : 
      21                 :            : ///@cond PRIVATE
      22                 :            : 
      23                 :          0 : QString QgsRasterSamplingAlgorithm::name() const
      24                 :            : {
      25                 :          0 :   return QStringLiteral( "rastersampling" );
      26                 :            : }
      27                 :            : 
      28                 :          0 : QString QgsRasterSamplingAlgorithm::displayName() const
      29                 :            : {
      30                 :          0 :   return QObject::tr( "Sample raster values" );
      31                 :            : }
      32                 :            : 
      33                 :          0 : QStringList QgsRasterSamplingAlgorithm::tags() const
      34                 :            : {
      35                 :          0 :   return QObject::tr( "extract,point,pixel,value" ).split( ',' );
      36                 :          0 : }
      37                 :            : 
      38                 :          0 : QString QgsRasterSamplingAlgorithm::group() const
      39                 :            : {
      40                 :          0 :   return QObject::tr( "Raster analysis" );
      41                 :            : }
      42                 :            : 
      43                 :          0 : QString QgsRasterSamplingAlgorithm::groupId() const
      44                 :            : {
      45                 :          0 :   return QStringLiteral( "rasteranalysis" );
      46                 :            : }
      47                 :            : 
      48                 :          0 : QString QgsRasterSamplingAlgorithm::shortDescription() const
      49                 :            : {
      50                 :          0 :   return QObject::tr( "Samples raster values under a set of points." );
      51                 :            : }
      52                 :            : 
      53                 :          0 : QString QgsRasterSamplingAlgorithm::shortHelpString() const
      54                 :            : {
      55                 :          0 :   return QObject::tr( "This algorithm creates a new vector layer with the same attributes of the input layer and the raster values corresponding on the point location.\n\n"
      56                 :            :                       "If the raster layer has more than one band, all the band values are sampled." );
      57                 :            : }
      58                 :            : 
      59                 :          0 : QgsRasterSamplingAlgorithm *QgsRasterSamplingAlgorithm::createInstance() const
      60                 :            : {
      61                 :          0 :   return new QgsRasterSamplingAlgorithm();
      62                 :          0 : }
      63                 :            : 
      64                 :          0 : void QgsRasterSamplingAlgorithm::initAlgorithm( const QVariantMap & )
      65                 :            : {
      66                 :          0 :   addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList< int >() << QgsProcessing::TypeVectorPoint ) );
      67                 :          0 :   addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "RASTERCOPY" ), QObject::tr( "Raster layer" ) ) );
      68                 :            : 
      69                 :          0 :   addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "SAMPLE_" ), false, true ) );
      70                 :          0 :   addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Sampled" ), QgsProcessing::TypeVectorPoint ) );
      71                 :          0 : }
      72                 :            : 
      73                 :          0 : bool QgsRasterSamplingAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
      74                 :            : {
      75                 :          0 :   QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "RASTERCOPY" ), context );
      76                 :          0 :   if ( !layer )
      77                 :          0 :     throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "RASTERCOPY" ) ) );
      78                 :            : 
      79                 :          0 :   mBandCount = layer->bandCount();
      80                 :          0 :   mCrs = layer->crs();
      81                 :          0 :   mDataProvider.reset( static_cast< QgsRasterDataProvider *>( layer->dataProvider()->clone() ) );
      82                 :            : 
      83                 :          0 :   return true;
      84                 :          0 : }
      85                 :            : 
      86                 :          0 : QVariantMap QgsRasterSamplingAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
      87                 :            : {
      88                 :          0 :   std::unique_ptr< QgsFeatureSource > source( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
      89                 :          0 :   if ( !source )
      90                 :          0 :     throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
      91                 :            : 
      92                 :          0 :   QString fieldPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context );
      93                 :          0 :   QgsFields newFields;
      94                 :          0 :   QgsAttributes emptySampleAttributes;
      95                 :          0 :   for ( int band = 1; band <= mBandCount; band++ )
      96                 :            :   {
      97                 :          0 :     Qgis::DataType dataType = mDataProvider->dataType( band );
      98                 :          0 :     bool intSafe = ( dataType == Qgis::Byte || dataType == Qgis::UInt16 || dataType == Qgis::Int16 || dataType == Qgis::UInt32 ||
      99                 :          0 :                      dataType == Qgis::Int32 || dataType == Qgis::CInt16 || dataType == Qgis::CInt32 );
     100                 :            : 
     101                 :          0 :     newFields.append( QgsField( QStringLiteral( "%1%2" ).arg( fieldPrefix, QString::number( band ) ), intSafe ? QVariant::Int : QVariant::Double ) );
     102                 :          0 :     emptySampleAttributes += QVariant();
     103                 :          0 :   }
     104                 :          0 :   QgsFields fields = QgsProcessingUtils::combineFields( source->fields(), newFields );
     105                 :            : 
     106                 :          0 :   QString dest;
     107                 :          0 :   std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields,
     108                 :          0 :                                           source->wkbType(), source->sourceCrs() ) );
     109                 :          0 :   if ( !sink )
     110                 :          0 :     throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
     111                 :            : 
     112                 :          0 :   const long count = source->featureCount();
     113                 :          0 :   const double step = count > 0 ? 100.0 / count : 1;
     114                 :          0 :   long current = 0;
     115                 :            : 
     116                 :          0 :   QgsCoordinateTransform ct( source->sourceCrs(), mCrs, context.transformContext() );
     117                 :          0 :   QgsFeatureIterator it = source->getFeatures( QgsFeatureRequest() );
     118                 :          0 :   QgsFeature feature;
     119                 :          0 :   while ( it.nextFeature( feature ) )
     120                 :            :   {
     121                 :          0 :     if ( feedback->isCanceled() )
     122                 :            :     {
     123                 :          0 :       break;
     124                 :            :     }
     125                 :          0 :     feedback->setProgress( current * step );
     126                 :          0 :     current++;
     127                 :            : 
     128                 :          0 :     QgsAttributes attributes = feature.attributes();
     129                 :          0 :     QgsFeature outputFeature( feature );
     130                 :          0 :     if ( !feature.hasGeometry() )
     131                 :            :     {
     132                 :          0 :       attributes += emptySampleAttributes;
     133                 :          0 :       outputFeature.setAttributes( attributes );
     134                 :          0 :       sink->addFeature( outputFeature, QgsFeatureSink::FastInsert );
     135                 :          0 :       feedback->reportError( QObject::tr( "No geometry attached to feature %1." ).arg( feature.id() ) );
     136                 :          0 :       continue;
     137                 :            :     }
     138                 :            : 
     139                 :          0 :     QgsGeometry geometry = feature.geometry();
     140                 :          0 :     if ( geometry.isMultipart() && geometry.get()->partCount() != 1 )
     141                 :            :     {
     142                 :          0 :       attributes += emptySampleAttributes;
     143                 :          0 :       outputFeature.setAttributes( attributes );
     144                 :          0 :       sink->addFeature( outputFeature, QgsFeatureSink::FastInsert );
     145                 :          0 :       feedback->reportError( QObject::tr( "Impossible to sample data of multipart feature %1." ).arg( feature.id() ) );
     146                 :          0 :       continue;
     147                 :            :     }
     148                 :          0 :     QgsPointXY point( *( geometry.isMultipart() ? qgsgeometry_cast< const QgsPoint * >( qgsgeometry_cast< const QgsMultiPoint * >( geometry.constGet() )->geometryN( 0 ) ) :
     149                 :          0 :                          qgsgeometry_cast< const QgsPoint * >( geometry.constGet() ) ) );
     150                 :            :     try
     151                 :            :     {
     152                 :          0 :       point = ct.transform( point );
     153                 :          0 :     }
     154                 :            :     catch ( const QgsException & )
     155                 :            :     {
     156                 :          0 :       attributes += emptySampleAttributes;
     157                 :          0 :       outputFeature.setAttributes( attributes );
     158                 :          0 :       sink->addFeature( outputFeature, QgsFeatureSink::FastInsert );
     159                 :          0 :       feedback->reportError( QObject::tr( "Could not reproject feature %1 to raster CRS." ).arg( feature.id() ) );
     160                 :            :       continue;
     161                 :          0 :     }
     162                 :            : 
     163                 :          0 :     for ( int band = 1; band <= mBandCount; band ++ )
     164                 :            :     {
     165                 :          0 :       bool ok = false;
     166                 :          0 :       double value = mDataProvider->sample( point, band, &ok );
     167                 :          0 :       attributes += ok ? value : QVariant();
     168                 :          0 :     }
     169                 :          0 :     outputFeature.setAttributes( attributes );
     170                 :          0 :     sink->addFeature( outputFeature, QgsFeatureSink::FastInsert );
     171                 :          0 :   }
     172                 :            : 
     173                 :          0 :   QVariantMap outputs;
     174                 :          0 :   outputs.insert( QStringLiteral( "OUTPUT" ), dest );
     175                 :          0 :   return outputs;
     176                 :          0 : }
     177                 :            : 
     178                 :            : ///@endcond