qgsalgorithmzonalhistogram.cpp

Go to the documentation of this file.

/***************************************************************************
                         qgsalgorithmzonalhistogram.cpp
                         ---------------------
    begin                : May, 2018
    copyright            : (C) 2018 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 "qgsalgorithmzonalhistogram.h"
 
#include "qgslogger.h"
#include "qgsrasteranalysisutils.h"
 
#include <QString>
 
using namespace Qt::StringLiterals;

 
QString QgsZonalHistogramAlgorithm::name() const
{
  return u"zonalhistogram"_s;
}
 
QString QgsZonalHistogramAlgorithm::displayName() const
{
  return QObject::tr( "Zonal histogram" );
}
 
QStringList QgsZonalHistogramAlgorithm::tags() const
{
  return QObject::tr( "raster,unique,values,count,area,statistics" ).split( ',' );
}
 
QString QgsZonalHistogramAlgorithm::group() const
{
  return QObject::tr( "Raster analysis" );
}
 
QString QgsZonalHistogramAlgorithm::groupId() const
{
  return u"rasteranalysis"_s;
}
 
void QgsZonalHistogramAlgorithm::initAlgorithm( const QVariantMap & )
{
  addParameter( new QgsProcessingParameterRasterLayer( u"INPUT_RASTER"_s, QObject::tr( "Raster layer" ) ) );
  addParameter( new QgsProcessingParameterBand( u"RASTER_BAND"_s, QObject::tr( "Band number" ), 1, u"INPUT_RASTER"_s ) );
 
  addParameter(
    new QgsProcessingParameterFeatureSource( u"INPUT_VECTOR"_s, QObject::tr( "Vector layer containing zones" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorPolygon ) )
  );
 
  addParameter( new QgsProcessingParameterString( u"COLUMN_PREFIX"_s, QObject::tr( "Output column prefix" ), u"HISTO_"_s, false, true ) );
 
  addParameter( new QgsProcessingParameterFeatureSink( u"OUTPUT"_s, QObject::tr( "Output zones" ), Qgis::ProcessingSourceType::VectorPolygon ) );
}
 
QString QgsZonalHistogramAlgorithm::shortHelpString() const
{
  return QObject::tr( "This algorithm appends fields representing counts of each unique value from a raster layer contained within zones defined as polygons." );
}
 
QString QgsZonalHistogramAlgorithm::shortDescription() const
{
  return QObject::tr( "Appends fields representing counts of each unique value from a raster layer contained within zones defined as polygons." );
}
 
QgsZonalHistogramAlgorithm *QgsZonalHistogramAlgorithm::createInstance() const
{
  return new QgsZonalHistogramAlgorithm();
}
 
bool QgsZonalHistogramAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
  QgsRasterLayer *layer = parameterAsRasterLayer( parameters, u"INPUT_RASTER"_s, context );
  if ( !layer )
    throw QgsProcessingException( invalidRasterError( parameters, u"INPUT_RASTER"_s ) );
 
  mRasterBand = parameterAsInt( parameters, u"RASTER_BAND"_s, context );
  mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( mRasterBand );
  mNodataValue = layer->dataProvider()->sourceNoDataValue( mRasterBand );
  mRasterInterface.reset( layer->dataProvider()->clone() );
  mRasterExtent = layer->extent();
  mCrs = layer->crs();
  mCellSizeX = std::abs( layer->rasterUnitsPerPixelX() );
  mCellSizeY = std::abs( layer->rasterUnitsPerPixelX() );
  mNbCellsXProvider = mRasterInterface->xSize();
  mNbCellsYProvider = mRasterInterface->ySize();
  Qgis::DataType dataType = mRasterInterface->dataType( mRasterBand );
 
  switch ( dataType )
  {
    case Qgis::DataType::Byte:
    case Qgis::DataType::Int16:
    case Qgis::DataType::UInt16:
    case Qgis::DataType::Int32:
    case Qgis::DataType::UInt32:
      break;
    default:
      feedback->pushWarning(
        QObject::tr(
          "The input raster is a floating-point raster. Such rasters are not suitable for use with zonal histogram algorithm.\n"
          "Please use Round raster or Reclassify by table tools to reduce number of decimal places or define histogram bins."
        )
      );
      break;
  }
 
  return true;
}
 
QVariantMap QgsZonalHistogramAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
{
  QGS_MARK_ALGORITHM_SOURCE
 
  std::unique_ptr<QgsFeatureSource> zones( parameterAsSource( parameters, u"INPUT_VECTOR"_s, context ) );
  if ( !zones )
    throw QgsProcessingException( invalidSourceError( parameters, u"INPUT_VECTOR"_s ) );
 
  const long count = zones->featureCount();
  const double step = count > 0 ? 100.0 / count : 1;
  long current = 0;
 
  QList<double> uniqueValues;
  QMap<QgsFeatureId, QHash<double, qgssize>> featuresUniqueValues;
 
  // First loop through the zones to build up a list of unique values across all zones to determine sink fields list
  QgsFeatureRequest request;
  request.setNoAttributes();
  if ( zones->sourceCrs() != mCrs )
  {
    request.setDestinationCrs( mCrs, context.transformContext() );
  }
  QgsFeatureIterator it = zones->getFeatures( request );
  QgsFeature f;
  while ( it.nextFeature( f ) )
  {
    if ( feedback->isCanceled() )
    {
      break;
    }
    feedback->setProgress( current * step );
 
    if ( !f.hasGeometry() )
    {
      current++;
      continue;
    }
 
    const QgsGeometry featureGeometry = f.geometry();
    const QgsRectangle featureRect = featureGeometry.boundingBox().intersect( mRasterExtent );
    if ( featureRect.isEmpty() )
    {
      current++;
      continue;
    }
 
    int nCellsX, nCellsY;
    QgsRectangle rasterBlockExtent;
    QgsRasterAnalysisUtils::cellInfoForBBox( mRasterExtent, featureRect, mCellSizeX, mCellSizeY, nCellsX, nCellsY, mNbCellsXProvider, mNbCellsYProvider, rasterBlockExtent );
 
    QHash<double, qgssize> fUniqueValues;
    QgsRasterAnalysisUtils::statisticsFromMiddlePointTest(
      mRasterInterface.get(), mRasterBand, featureGeometry, nCellsX, nCellsY, mCellSizeX, mCellSizeY, rasterBlockExtent, [&fUniqueValues]( double value, const QgsPointXY & ) { fUniqueValues[value]++; }, false
    );
 
    if ( fUniqueValues.count() < 1 )
    {
      // The cell resolution is probably larger than the polygon area. We switch to slower precise pixel - polygon intersection in this case
      // TODO: eventually deal with weight if needed
      QgsRasterAnalysisUtils::statisticsFromPreciseIntersection(
        mRasterInterface.get(),
        mRasterBand,
        featureGeometry,
        nCellsX,
        nCellsY,
        mCellSizeX,
        mCellSizeY,
        rasterBlockExtent,
        [&fUniqueValues]( double value, double, const QgsPointXY & ) { fUniqueValues[value]++; },
        false
      );
    }
 
    for ( auto it = fUniqueValues.constBegin(); it != fUniqueValues.constEnd(); ++it )
    {
      if ( uniqueValues.indexOf( it.key() ) == -1 )
      {
        uniqueValues << it.key();
      }
      featuresUniqueValues[f.id()][it.key()] += it.value();
    }
 
    current++;
  }
 
  std::sort( uniqueValues.begin(), uniqueValues.end() );
 
  const QString fieldPrefix = parameterAsString( parameters, u"COLUMN_PREFIX"_s, context );
  QgsFields newFields;
  for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
  {
    newFields.append( QgsField( u"%1%2"_s.arg( fieldPrefix, mHasNoDataValue && *it == mNodataValue ? u"NODATA"_s : QString::number( *it ) ), QMetaType::Type::LongLong, QString(), -1, 0 ) );
  }
  const QgsFields fields = QgsProcessingUtils::combineFields( zones->fields(), newFields );
 
  QString dest;
  std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, u"OUTPUT"_s, context, dest, fields, zones->wkbType(), zones->sourceCrs() ) );
  if ( !sink )
    throw QgsProcessingException( invalidSinkError( parameters, u"OUTPUT"_s ) );
 
  it = zones->getFeatures( QgsFeatureRequest() );
  while ( it.nextFeature( f ) )
  {
    QgsAttributes attributes = f.attributes();
    const QHash<double, qgssize> fUniqueValues = featuresUniqueValues.value( f.id() );
    for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it )
    {
      attributes += fUniqueValues.value( *it, 0 );
    }
 
    QgsFeature outputFeature;
    outputFeature.setGeometry( f.geometry() );
    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;
}