add new raster stack position algorithms

qgis · Aug 5, 2020 · 06f9840 · 06f9840
1 parent a56604a
commit 06f9840
Show file tree

Hide file tree

Showing 7 changed files with 378 additions and 0 deletions.
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
@@ -132,6 +132,7 @@ SET(QGIS_ANALYSIS_SRCS
   processing/qgsalgorithmrasterlayeruniquevalues.cpp
   processing/qgsalgorithmrasterlogicalop.cpp
   processing/qgsalgorithmrasterize.cpp
+  processing/qgsalgorithmrasterstackposition.cpp
   processing/qgsalgorithmrasterstatistics.cpp
   processing/qgsalgorithmrastersurfacevolume.cpp
   processing/qgsalgorithmrasterzonalstats.cpp

diff --git a/src/analysis/processing/qgsalgorithmrasterstackposition.cpp b/src/analysis/processing/qgsalgorithmrasterstackposition.cpp
@@ -0,0 +1,288 @@
+/***************************************************************************
+                         qgsrasterstackposition.cpp
+                         ---------------------
+    begin                : July 2020
+    copyright            : (C) 2020 by Clemens Raffler
+    email                : clemens dot raffler 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 "qgsalgorithmrasterstackposition.h"
+#include "qgsrasterprojector.h"
+#include "qgsrasterfilewriter.h"
+#include "qgsrasteranalysisutils.h"
+
+///@cond PRIVATE
+
+//
+//QgsRasterFrequencyByComparisonOperatorBase
+//
+
+QString QgsRasterStackPositionAlgorithmBase::group() const
+{
+  return QObject::tr( "Raster analysis" );
+}
+
+QString QgsRasterStackPositionAlgorithmBase::groupId() const
+{
+  return QStringLiteral( "rasteranalysis" );
+}
+
+void QgsRasterStackPositionAlgorithmBase::initAlgorithm( const QVariantMap & )
+{
+  addParameter( new QgsProcessingParameterMultipleLayers( QStringLiteral( "INPUT_RASTERS" ),
+                QObject::tr( "Input raster layers" ), QgsProcessing::TypeRaster ) );
+
+  addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "REFERENCE_LAYER" ), QObject::tr( "Reference layer" ) ) );
+
+  addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "IGNORE_NODATA" ), QObject::tr( "Ignore NoData values" ), false ) );
+
+  std::unique_ptr< QgsProcessingParameterNumber > output_nodata_parameter = qgis::make_unique< QgsProcessingParameterNumber >( QStringLiteral( "OUTPUT_NODATA_VALUE" ), QObject::tr( "Output NoData value" ), QgsProcessingParameterNumber::Double, -9999, true );
+  output_nodata_parameter->setFlags( output_nodata_parameter->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
+  addParameter( output_nodata_parameter.release() );
+
+  addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ),
+                QObject::tr( "Output layer" ) ) );
+  addOutput( new QgsProcessingOutputString( QStringLiteral( "EXTENT" ), QObject::tr( "Extent" ) ) );
+  addOutput( new QgsProcessingOutputString( QStringLiteral( "CRS_AUTHID" ), QObject::tr( "CRS authority identifier" ) ) );
+  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "WIDTH_IN_PIXELS" ), QObject::tr( "Width in pixels" ) ) );
+  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "HEIGHT_IN_PIXELS" ), QObject::tr( "Height in pixels" ) ) );
+  addOutput( new QgsProcessingOutputNumber( QStringLiteral( "TOTAL_PIXEL_COUNT" ), QObject::tr( "Total pixel count" ) ) );
+}
+
+bool QgsRasterStackPositionAlgorithmBase::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  QgsRasterLayer *referenceLayer = parameterAsRasterLayer( parameters, QStringLiteral( "REFERENCE_LAYER" ), context );
+  if ( !referenceLayer )
+    throw QgsProcessingException( invalidRasterError( parameters, QStringLiteral( "REFERENCE_LAYER" ) ) );
+
+  mIgnoreNoData = parameterAsBool( parameters, QStringLiteral( "IGNORE_NODATA" ), context );
+  mNoDataValue = parameterAsDouble( parameters, QStringLiteral( "OUTPUT_NODATA_VALUE" ), context );
+
+  mCrs = referenceLayer->crs();
+  mRasterUnitsPerPixelX = referenceLayer->rasterUnitsPerPixelX();
+  mRasterUnitsPerPixelY = referenceLayer->rasterUnitsPerPixelY();
+  mLayerWidth = referenceLayer->width();
+  mLayerHeight = referenceLayer->height();
+  mExtent = referenceLayer->extent();
+
+  const QList< QgsMapLayer * > layers = parameterAsLayerList( parameters, QStringLiteral( "INPUT_RASTERS" ), context );
+  QList< QgsRasterLayer * > rasterLayers;
+  rasterLayers.reserve( layers.count() );
+  for ( QgsMapLayer *l : layers )
+  {
+    if ( feedback->isCanceled() )
+      break; //in case some slow data sources are loaded
+
+    if ( l->type() == QgsMapLayerType::RasterLayer )
+    {
+      QgsRasterLayer *layer = qobject_cast< QgsRasterLayer * >( l );
+      QgsRasterAnalysisUtils::RasterLogicInput input;
+      const int band = 1; //could be made dynamic
+      input.hasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( band );
+      input.sourceDataProvider.reset( layer->dataProvider()->clone() );
+      input.interface = input.sourceDataProvider.get();
+      // add projector if necessary
+      if ( layer->crs() != mCrs )
+      {
+        input.projector = qgis::make_unique< QgsRasterProjector >();
+        input.projector->setInput( input.sourceDataProvider.get() );
+        input.projector->setCrs( layer->crs(), mCrs, context.transformContext() );
+        input.interface = input.projector.get();
+      }
+      mInputs.emplace_back( std::move( input ) );
+    }
+  }
+
+  return true;
+}
+
+QVariantMap QgsRasterStackPositionAlgorithmBase::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
+  QFileInfo fi( outputFile );
+  const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
+
+  std::unique_ptr< QgsRasterFileWriter > writer = qgis::make_unique< QgsRasterFileWriter >( outputFile );
+  writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
+  writer->setOutputFormat( outputFormat );
+  std::unique_ptr<QgsRasterDataProvider > provider( writer->createOneBandRaster( Qgis::Int32, mLayerWidth, mLayerHeight, mExtent, mCrs ) );
+  if ( !provider )
+    throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
+  if ( !provider->isValid() )
+    throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
+
+  provider->setNoDataValue( 1, mNoDataValue );
+  qgssize layerSize = static_cast< qgssize >( mLayerWidth ) * static_cast< qgssize >( mLayerHeight );
+
+  int maxWidth = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_WIDTH;
+  int maxHeight = QgsRasterIterator::DEFAULT_MAXIMUM_TILE_HEIGHT;
+  int nbBlocksWidth = static_cast< int>( std::ceil( 1.0 * mLayerWidth / maxWidth ) );
+  int nbBlocksHeight = static_cast< int >( std::ceil( 1.0 * mLayerHeight / maxHeight ) );
+  int nbBlocks = nbBlocksWidth * nbBlocksHeight;
+  provider->setEditable( true );
+
+  QgsRasterIterator iter( provider.get() );
+  iter.startRasterRead( 1, mLayerWidth, mLayerHeight, mExtent );
+  int iterLeft = 0;
+  int iterTop = 0;
+  int iterCols = 0;
+  int iterRows = 0;
+  QgsRectangle blockExtent;
+
+  std::unique_ptr< QgsRasterBlock > outputBlock;
+  while ( iter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
+  {
+    std::vector< std::unique_ptr< QgsRasterBlock > > inputBlocks;
+    for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : mInputs )
+    {
+      if ( feedback->isCanceled() )
+        break; //in case some slow data sources are loaded
+      for ( int band : i.bands )
+      {
+        if ( feedback->isCanceled() )
+          break; //in case some slow data sources are loaded
+        std::unique_ptr< QgsRasterBlock > b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
+        inputBlocks.emplace_back( std::move( b ) );
+      }
+    }
+
+    feedback->setProgress( 100 * ( ( iterTop / maxHeight * nbBlocksWidth ) + iterLeft / maxWidth ) / nbBlocks );
+    for ( int row = 0; row < iterRows; row++ )
+    {
+      if ( feedback->isCanceled() )
+        break;
+
+      for ( int col = 0; col < iterCols; col++ )
+      {
+        bool noDataInStack = false;
+        //rewrite getCellValuesFromBlock
+        std::vector<double> cellValues = QgsRasterAnalysisUtils::getCellValuesFromBlockStack( inputBlocks, row, col, noDataInStack );
+
+        if ( noDataInStack && !mIgnoreNoData )
+        {
+          //output cell will always be NoData if NoData occurs the cellValueStack and NoData is not ignored
+          //this saves unnecessary iterations on the cellValueStack
+          outputBlock->setValue( row, col, mNoDataValue );
+        }
+        else
+        {
+          int position = getPosition( cellValues );
+          outputBlock->setValue( row, col, position );
+        }
+      }
+    }
+    provider->writeBlock( outputBlock.get(), 1, iterLeft, iterTop );
+  }
+  provider->setEditable( false );
+
+  QVariantMap outputs;
+  outputs.insert( QStringLiteral( "EXTENT" ), mExtent.toString() );
+  outputs.insert( QStringLiteral( "CRS_AUTHID" ), mCrs.authid() );
+  outputs.insert( QStringLiteral( "WIDTH_IN_PIXELS" ), mLayerWidth );
+  outputs.insert( QStringLiteral( "HEIGHT_IN_PIXELS" ), mLayerHeight );
+  outputs.insert( QStringLiteral( "TOTAL_PIXEL_COUNT" ), layerSize );
+  outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
+
+  return outputs;
+}
+
+//
+// QgsRasterStackLowestPositionAlgorithm
+//
+QString QgsRasterStackLowestPositionAlgorithm::displayName() const
+{
+  return QObject::tr( "Lowest position in raster stack" );
+}
+
+QString QgsRasterStackLowestPositionAlgorithm::name() const
+{
+  return QObject::tr( "lowestpositioninrasterstack" );
+}
+
+QStringList QgsRasterStackLowestPositionAlgorithm::tags() const
+{
+  return QObject::tr( "cell,lowest,position,pixel,stack" ).split( ',' );
+}
+
+QString QgsRasterStackLowestPositionAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "The Lowest position algorithm evaluates on a cell-by-cell basis the position "
+                      "of the raster with the highest value in a stack of rasters. Position counts start "
+                      "with 1 and range to the total number input rasters. The order of the input "
+                      "rasters is relevant for the algorithm. If multiple rasters feature the lowest value, "
+                      "the first raster will be used for the position value.\n "
+                      "If multiband rasters are used in the data raster stack, the algorithm will always "
+                      "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
+                      "Any NoData cells in the raster layer stack will result in a NoData cell "
+                      "in the output raster unless the ignore NoData parameter is checked. "
+                      "The output NoData value can be set manually. The output rasters extent and resolution "
+                      "is defined by a reference raster layer and is always of int32 type." );
+}
+
+QgsRasterStackLowestPositionAlgorithm *QgsRasterStackLowestPositionAlgorithm::createInstance() const
+{
+  return new QgsRasterStackLowestPositionAlgorithm();
+}
+
+int QgsRasterStackLowestPositionAlgorithm::getPosition( std::vector<QgsRasterBlock> inputBlocks )
+{
+  //the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
+  return static_cast<int>( std::min_element( cellValueStack.begin(), cellValueStack.end() ) - cellValueStack.begin() ) + 1; //therefore... index + 1
+}
+
+//
+// QgsRasterStackHighestPositionAlgorithmAlgorithm
+//
+
+QString QgsRasterStackHighestPositionAlgorithm::displayName() const
+{
+  return QObject::tr( "Highest position in raster stack" );
+}
+
+QString QgsRasterStackHighestPositionAlgorithm::name() const
+{
+  return QObject::tr( "highestpositioninrasterstack" );
+}
+
+QStringList QgsRasterStackHighestPositionAlgorithm::tags() const
+{
+  return QObject::tr( "cell,highest,position,pixel,stack" ).split( ',' );
+}
+
+QString QgsRasterStackHighestPositionAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "The Highest position algorithm evaluates on a cell-by-cell basis the position "
+                      "of the raster with the highest value in a stack of rasters. Position counts start "
+                      "with 1 and range to the total number input rasters. The order of the input "
+                      "rasters is relevant for the algorithm. If multiple rasters feature the highest value, "
+                      "the first raster will be used for the position value.\n "
+                      "If multiband rasters are used in the data raster stack, the algorithm will always "
+                      "perform the analysis on the first band of the rasters - use GDAL to use other bands in the analysis. "
+                      "Any NoData cells in the raster layer stack will result in a NoData cell "
+                      "in the output raster unless the ignore NoData parameter is checked. "
+                      "The output NoData value can be set manually. The output rasters extent and resolution "
+                      "is defined by a reference raster layer and is always of int32 type." );
+}
+
+QgsRasterStackHighestPositionAlgorithm *QgsRasterStackHighestPositionAlgorithm::createInstance() const
+{
+  return new QgsRasterStackHighestPositionAlgorithm();
+}
+
+int QgsRasterStackHighestPositionAlgorithm::getPosition( std::vector<QgsRasterBlock> inputBlocks )
+{
+  //the ArcGIS implementation uses 1 for first position value instead of 0 as in standard c++
+  return static_cast<int>( std::max_element( cellValueStack.begin(), cellValueStack.end() ) - cellValueStack.begin() ) + 1; //therefore... index + 1
+}
+
+///@endcond
+
diff --git a/src/analysis/processing/qgsalgorithmrasterstackposition.h b/src/analysis/processing/qgsalgorithmrasterstackposition.h
@@ -0,0 +1,86 @@
+/***************************************************************************
+                         qgsalgorithmrasterstackposition.h
+                         ---------------------
+    begin                : July 2020
+    copyright            : (C) 2020 by Clemens Raffler
+    email                : clemens dot raffler 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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSALGORITHMRASTERSTACKPOSITION_H
+#define QGSALGORITHMRASTERSTACKPOSITION_H
+
+#define SIP_NO_FILE
+
+#include "qgis_sip.h"
+#include "qgsapplication.h"
+#include "qgsprocessingalgorithm.h"
+#include "qgsrasterprojector.h"
+#include "qgsrasteranalysisutils.h"
+
+///@cond PRIVATE
+
+class QgsRasterStackPositionAlgorithmBase : public QgsProcessingAlgorithm
+{
+  public:
+    QgsRasterStackPositionAlgorithmBase() = default;
+    void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
+    QString group() const override;
+    QString groupId() const override;
+
+  protected:
+    bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+    QVariantMap processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+    virtual int getPosition( std::vector<QgsRasterBlock>inputBlocks ) = 0;
+
+  private:
+    std::vector< QgsRasterAnalysisUtils::RasterLogicInput > mInputs;
+    bool mIgnoreNoData;
+    double mNoDataValue = -9999;
+    int mLayerWidth;
+    int mLayerHeight;
+    QgsRectangle mExtent;
+    QgsCoordinateReferenceSystem mCrs;
+    double mRasterUnitsPerPixelX;
+    double mRasterUnitsPerPixelY;
+};
+
+class QgsRasterStackLowestPositionAlgorithm : public QgsRasterStackPositionAlgorithmBase
+{
+  public:
+    QgsRasterStackLowestPositionAlgorithm() = default;
+    QString name() const override;
+    QString displayName() const override;
+    QStringList tags() const override;
+    QString shortHelpString() const override;
+    QgsRasterStackLowestPositionAlgorithm *createInstance() const override SIP_FACTORY;
+
+  protected:
+    int getPosition( std::vector<QgsRasterBlock>inputBlocks ) override;
+};
+
+class QgsRasterStackHighestPositionAlgorithm : public QgsRasterStackPositionAlgorithmBase
+{
+  public:
+    QgsRasterStackHighestPositionAlgorithm() = default;
+    QString name() const override;
+    QString displayName() const override;
+    QStringList tags() const override;
+    QString shortHelpString() const override;
+    QgsRasterStackHighestPositionAlgorithm *createInstance() const override SIP_FACTORY;
+
+  protected:
+    int getPosition( std::vector<QgsRasterBlock>inputBlocks ) override;
+};
+
+///@endcond PRIVATE
+
+#endif // QGSALGORITHMRASTERSTACKPOSITION_H