qgsalgorithmrastersampling.cpp

Go to the documentation of this file.

/***************************************************************************
                         qgsalgorithmrastersampling.cpp
                         --------------------------
    begin                : August 2020
    copyright            : (C) 2020 by Mathieu Pellerin
    email                : nirvn dot asia at gmail dot com
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "qgsalgorithmrastersampling.h"
 
#include "qgsmultipoint.h"
 
#include <QString>
 
using namespace Qt::StringLiterals;

 
QString QgsRasterSamplingAlgorithm::name() const
{
  return u"rastersampling"_s;
}
 
QString QgsRasterSamplingAlgorithm::displayName() const
{
  return QObject::tr( "Sample raster values" );
}
 
QStringList QgsRasterSamplingAlgorithm::tags() const
{
  return QObject::tr( "extract,point,pixel,sample,value" ).split( ',' );
}
 
QString QgsRasterSamplingAlgorithm::group() const
{
  return QObject::tr( "Raster analysis" );
}
 
QString QgsRasterSamplingAlgorithm::groupId() const
{
  return u"rasteranalysis"_s;
}
 
QString QgsRasterSamplingAlgorithm::shortDescription() const
{
  return QObject::tr( "Samples raster values under a set of points." );
}
 
QString QgsRasterSamplingAlgorithm::shortHelpString() const
{
  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"
    "If the raster layer has more than one band, all the band values are sampled."
  );
}
 
QgsRasterSamplingAlgorithm *QgsRasterSamplingAlgorithm::createInstance() const
{
  return new QgsRasterSamplingAlgorithm();
}
 
void QgsRasterSamplingAlgorithm::initAlgorithm( const QVariantMap & )
{
  addParameter( new QgsProcessingParameterFeatureSource( u"INPUT"_s, QObject::tr( "Input layer" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorPoint ) ) );
  addParameter( new QgsProcessingParameterRasterLayer( u"RASTERCOPY"_s, QObject::tr( "Raster layer" ) ) );
 
  addParameter( new QgsProcessingParameterString( u"COLUMN_PREFIX"_s, QObject::tr( "Output column prefix" ), u"SAMPLE_"_s, false, true ) );
  addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT"_s, QObject::tr( "Sampled" ), Qgis::ProcessingSourceType::VectorPoint ) );
}
 
bool QgsRasterSamplingAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
{
  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, u"RASTERCOPY"_s, context );
  if ( !layer )
    throw QgsProcessingException( invalidRasterError( parameters, u"RASTERCOPY"_s ) );
 
  mBandCount = layer->bandCount();
  mCrs = layer->crs();
  mDataProvider.reset( static_cast<QgsRasterDataProvider *>( layer->dataProvider()->clone() ) );
 
  return true;
}
 
QVariantMap QgsRasterSamplingAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
  QGS_MARK_ALGORITHM_SOURCE
 
  std::unique_ptr<QgsFeatureSource> source( parameterAsSource( parameters, u"INPUT"_s, context ) );
  if ( !source )
    throw QgsProcessingException( invalidSourceError( parameters, u"INPUT"_s ) );
 
  const QString fieldPrefix = parameterAsString( parameters, u"COLUMN_PREFIX"_s, context );
  QgsFields newFields;
  QgsAttributes emptySampleAttributes;
  for ( int band = 1; band <= mBandCount; band++ )
  {
    const Qgis::DataType dataType = mDataProvider->dataType( band );
    const bool intSafe
      = ( dataType == Qgis::DataType::Byte || dataType == Qgis::DataType::UInt16 || dataType == Qgis::DataType::Int16 || dataType == Qgis::DataType::UInt32 || dataType == Qgis::DataType::Int32 || dataType == Qgis::DataType::CInt16 || dataType == Qgis::DataType::CInt32 );
 
    newFields.append( QgsField( u"%1%2"_s.arg( fieldPrefix, QString::number( band ) ), intSafe ? QMetaType::Type::Int : QMetaType::Type::Double ) );
    emptySampleAttributes += QVariant();
  }
  const QgsFields fields = QgsProcessingUtils::combineFields( source->fields(), newFields );
 
  QString dest;
  std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, u"OUTPUT"_s, context, dest, fields, source->wkbType(), source->sourceCrs() ) );
  if ( !sink )
    throw QgsProcessingException( invalidSinkError( parameters, u"OUTPUT"_s ) );
 
  const long count = source->featureCount();
  const double step = count > 0 ? 100.0 / count : 1;
  long current = 0;
 
  const QgsCoordinateTransform ct( source->sourceCrs(), mCrs, context.transformContext() );
  QgsFeatureIterator it = source->getFeatures( QgsFeatureRequest() );
  QgsFeature feature;
  while ( it.nextFeature( feature ) )
  {
    if ( feedback->isCanceled() )
    {
      break;
    }
    feedback->setProgress( current * step );
    current++;
 
    QgsAttributes attributes = feature.attributes();
    QgsFeature outputFeature( feature );
    if ( !feature.hasGeometry() )
    {
      attributes += emptySampleAttributes;
      outputFeature.setAttributes( attributes );
      if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
        throw QgsProcessingException( writeFeatureError( sink.get(), parameters, u"OUTPUT"_s ) );
      else
        feedback->featureAddedToSink( u"OUTPUT"_s );
      feedback->reportError( QObject::tr( "No geometry attached to feature %1." ).arg( feature.id() ) );
      continue;
    }
 
    QgsGeometry geometry = feature.geometry();
    if ( geometry.isMultipart() && geometry.get()->partCount() != 1 )
    {
      attributes += emptySampleAttributes;
      outputFeature.setAttributes( attributes );
      if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
        throw QgsProcessingException( writeFeatureError( sink.get(), parameters, u"OUTPUT"_s ) );
      else
        feedback->featureAddedToSink( u"OUTPUT"_s );
      feedback->reportError( QObject::tr( "Impossible to sample data of multipart feature %1." ).arg( feature.id() ) );
      continue;
    }
    QgsPointXY point(
      *( geometry.isMultipart() ? qgsgeometry_cast<const QgsPoint *>( qgsgeometry_cast<const QgsMultiPoint *>( geometry.constGet() )->geometryN( 0 ) )
                                : qgsgeometry_cast<const QgsPoint *>( geometry.constGet() ) )
    );
    try
    {
      point = ct.transform( point );
    }
    catch ( const QgsException & )
    {
      attributes += emptySampleAttributes;
      outputFeature.setAttributes( attributes );
      if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
        throw QgsProcessingException( writeFeatureError( sink.get(), parameters, u"OUTPUT"_s ) );
      else
        feedback->featureAddedToSink( u"OUTPUT"_s );
      feedback->reportError( QObject::tr( "Could not reproject feature %1 to raster CRS." ).arg( feature.id() ) );
      continue;
    }
 
    for ( int band = 1; band <= mBandCount; band++ )
    {
      bool ok = false;
      const double value = mDataProvider->sample( point, band, &ok );
      attributes += ok ? value : QVariant();
    }
    outputFeature.setAttributes( attributes );
    if ( !sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ) )
      throw QgsProcessingException( writeFeatureError( sink.get(), parameters, u"OUTPUT"_s ) );
    else
      feedback->featureAddedToSink( u"OUTPUT"_s );
  }
 
  sink->finalize();
  feedback->featureSinkFinalized( u"OUTPUT"_s );
 
  QVariantMap outputs;
  outputs.insert( u"OUTPUT"_s, dest );
  return outputs;
}