qgsrasteranalysisutils.cpp

Go to the documentation of this file.

/***************************************************************************
  qgsrasteranalysisutils.cpp
  ---------------------
  Date                 : June 2018
  Copyright            : (C) 2018 by Nyall Dawson
  Email                : nyall dot dawson 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 "qgsrasteranalysisutils.h"
 
#include <cmath>
#include <map>
#include <unordered_map>
#include <unordered_set>
 
#include "qgsfeedback.h"
#include "qgsgeos.h"
#include "qgsprocessingparameters.h"
#include "qgsrasterblock.h"
#include "qgsrasteriterator.h"
 
#include <QString>
 
using namespace Qt::StringLiterals;

 
void QgsRasterAnalysisUtils::cellInfoForBBox(
  const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY, int &nCellsX, int &nCellsY, int rasterWidth, int rasterHeight, QgsRectangle &rasterBlockExtent
)
{
  //get intersecting bbox
  const QgsRectangle intersectBox = rasterBBox.intersect( featureBBox );
  if ( intersectBox.isEmpty() )
  {
    nCellsX = 0;
    nCellsY = 0;
    rasterBlockExtent = QgsRectangle();
    return;
  }
 
  //get offset in pixels in x- and y- direction
  const int offsetX = static_cast<int>( std::floor( ( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX ) );
  const int offsetY = static_cast<int>( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY ) );
 
  const int maxColumn = static_cast<int>( std::floor( ( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) ) + 1;
  const int maxRow = static_cast<int>( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) ) + 1;
 
  nCellsX = maxColumn - offsetX;
  nCellsY = maxRow - offsetY;
 
  //avoid access to cells outside of the raster (may occur because of rounding)
  nCellsX = std::min( offsetX + nCellsX, rasterWidth ) - offsetX;
  nCellsY = std::min( offsetY + nCellsY, rasterHeight ) - offsetY;
 
  rasterBlockExtent
    = QgsRectangle( rasterBBox.xMinimum() + offsetX * cellSizeX, rasterBBox.yMaximum() - offsetY * cellSizeY, rasterBBox.xMinimum() + ( nCellsX + offsetX ) * cellSizeX, rasterBBox.yMaximum() - ( nCellsY + offsetY ) * cellSizeY );
}
 
void QgsRasterAnalysisUtils::statisticsFromMiddlePointTest(
  QgsRasterInterface *rasterInterface,
  int rasterBand,
  const QgsGeometry &poly,
  int nCellsX,
  int nCellsY,
  double cellSizeX,
  double cellSizeY,
  const QgsRectangle &rasterBBox,
  const std::function<void( double, const QgsPointXY & )> &addValue,
  bool skipNodata
)
{
  if ( !poly.constGet() )
    return;
 
  auto polyEngine = std::make_unique<QgsGeos>( poly.constGet() );
  polyEngine->prepareGeometry();
 
  QgsRasterIterator iter( rasterInterface );
  iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
 
  std::unique_ptr<QgsRasterBlock> block;
  int iterLeft = 0;
  int iterTop = 0;
  int iterCols = 0;
  int iterRows = 0;
  QgsRectangle blockExtent;
  bool isNoData = false;
  while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
  {
    double cellCenterY = blockExtent.yMaximum() - 0.5 * cellSizeY;
 
    for ( int row = 0; row < iterRows; ++row )
    {
      double cellCenterX = blockExtent.xMinimum() + 0.5 * cellSizeX;
      for ( int col = 0; col < iterCols; ++col )
      {
        const double pixelValue = block->valueAndNoData( row, col, isNoData );
        if ( validPixel( pixelValue ) && ( !skipNodata || !isNoData ) )
        {
          if ( polyEngine->contains( cellCenterX, cellCenterY ) )
          {
            addValue( pixelValue, QgsPointXY( cellCenterX, cellCenterY ) );
          }
        }
        cellCenterX += cellSizeX;
      }
      cellCenterY -= cellSizeY;
    }
  }
}
 
void QgsRasterAnalysisUtils::statisticsFromPreciseIntersection(
  QgsRasterInterface *rasterInterface,
  int rasterBand,
  const QgsGeometry &poly,
  int nCellsX,
  int nCellsY,
  double cellSizeX,
  double cellSizeY,
  const QgsRectangle &rasterBBox,
  const std::function<void( double, double, const QgsPointXY & )> &addValue,
  bool skipNodata
)
{
  QgsGeometry pixelRectGeometry;
 
  const double hCellSizeX = cellSizeX / 2.0;
  const double hCellSizeY = cellSizeY / 2.0;
  const double pixelArea = cellSizeX * cellSizeY;
  double weight = 0;
 
  std::unique_ptr<QgsGeometryEngine> polyEngine( QgsGeometry::createGeometryEngine( poly.constGet() ) );
  if ( !polyEngine )
  {
    return;
  }
  polyEngine->prepareGeometry();
 
  QgsRasterIterator iter( rasterInterface );
  iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
 
  std::unique_ptr<QgsRasterBlock> block;
  int iterLeft = 0;
  int iterTop = 0;
  int iterCols = 0;
  int iterRows = 0;
  QgsRectangle blockExtent;
  bool isNoData = false;
  while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
  {
    double currentY = blockExtent.yMaximum() - 0.5 * cellSizeY;
    for ( int row = 0; row < iterRows; ++row )
    {
      double currentX = blockExtent.xMinimum() + 0.5 * cellSizeX;
      for ( int col = 0; col < iterCols; ++col )
      {
        const double pixelValue = block->valueAndNoData( row, col, isNoData );
        if ( validPixel( pixelValue ) && ( !skipNodata || !isNoData ) )
        {
          pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
          // GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
          // so we first test to see if there IS an intersection before doing the actual calculation
          if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
          {
            //intersection
            const QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
            if ( !intersectGeometry.isEmpty() )
            {
              const double intersectionArea = intersectGeometry.area();
              if ( intersectionArea > 0.0 )
              {
                weight = intersectionArea / pixelArea;
                addValue( pixelValue, weight, QgsPointXY( currentX, currentY ) );
              }
            }
          }
        }
        currentX += cellSizeX;
      }
      currentY -= cellSizeY;
    }
  }
}
 
bool QgsRasterAnalysisUtils::validPixel( double value )
{
  return !std::isnan( value );
}
 
void QgsRasterAnalysisUtils::mapToPixel( const double x, const double y, const QgsRectangle bounds, const double unitsPerPixelX, const double unitsPerPixelY, int &px, int &py )
{
  px = trunc( ( x - bounds.xMinimum() ) / unitsPerPixelX );
  py = trunc( ( y - bounds.yMaximum() ) / -unitsPerPixelY );
}
 
void QgsRasterAnalysisUtils::pixelToMap( const int px, const int py, const QgsRectangle bounds, const double unitsPerPixelX, const double unitsPerPixelY, double &x, double &y )
{
  x = bounds.xMinimum() + ( px + 0.5 ) * unitsPerPixelX;
  y = bounds.yMaximum() - ( py + 0.5 ) * unitsPerPixelY;
}
 
static QVector<QPair<QString, Qgis::DataType>> sDataTypes;
 
void populateDataTypes()
{
  if ( sDataTypes.empty() )
  {
    sDataTypes.append( qMakePair( u"Byte"_s, Qgis::DataType::Byte ) );
    sDataTypes.append( qMakePair( u"Int16"_s, Qgis::DataType::Int16 ) );
    sDataTypes.append( qMakePair( u"UInt16"_s, Qgis::DataType::UInt16 ) );
    sDataTypes.append( qMakePair( u"Int32"_s, Qgis::DataType::Int32 ) );
    sDataTypes.append( qMakePair( u"UInt32"_s, Qgis::DataType::UInt32 ) );
    sDataTypes.append( qMakePair( u"Float32"_s, Qgis::DataType::Float32 ) );
    sDataTypes.append( qMakePair( u"Float64"_s, Qgis::DataType::Float64 ) );
    sDataTypes.append( qMakePair( u"CInt16"_s, Qgis::DataType::CInt16 ) );
    sDataTypes.append( qMakePair( u"CInt32"_s, Qgis::DataType::CInt32 ) );
    sDataTypes.append( qMakePair( u"CFloat32"_s, Qgis::DataType::CFloat32 ) );
    sDataTypes.append( qMakePair( u"CFloat64"_s, Qgis::DataType::CFloat64 ) );
    sDataTypes.append( qMakePair( u"Int8"_s, Qgis::DataType::Int8 ) );
  }
}
 
std::unique_ptr<QgsProcessingParameterDefinition> QgsRasterAnalysisUtils::createRasterTypeParameter( const QString &name, const QString &description, Qgis::DataType defaultType )
{
  populateDataTypes();
 
  QStringList names;
  int defaultChoice = 0;
  int i = 0;
  for ( auto it = sDataTypes.constBegin(); it != sDataTypes.constEnd(); ++it )
  {
    names.append( it->first );
    if ( it->second == defaultType )
      defaultChoice = i;
    i++;
  }
 
  return std::make_unique<QgsProcessingParameterEnum>( name, description, names, false, defaultChoice );
}
 
Qgis::DataType QgsRasterAnalysisUtils::rasterTypeChoiceToDataType( int choice )
{
  if ( choice < 0 || choice >= sDataTypes.count() )
    return Qgis::DataType::Float32;
 
  return sDataTypes.value( choice ).second;
}
 
void QgsRasterAnalysisUtils::applyRasterLogicOperator(
  const std::vector<QgsRasterAnalysisUtils::RasterLogicInput> &inputs,
  std::unique_ptr<QgsRasterDataProvider> destinationRaster,
  double outputNoDataValue,
  const bool treatNoDataAsFalse,
  int width,
  int height,
  const QgsRectangle &extent,
  QgsFeedback *feedback,
  std::function<void( const std::vector<std::unique_ptr<QgsRasterBlock>> &, bool &, bool &, int, int, bool )> &applyLogicFunc,
  qgssize &noDataCount,
  qgssize &trueCount,
  qgssize &falseCount
)
{
  destinationRaster->setEditable( true );
  QgsRasterIterator outputIter( destinationRaster.get() );
  outputIter.startRasterRead( 1, width, height, extent );
 
  const bool hasReportsDuringClose = destinationRaster->hasReportsDuringClose();
  const double maxProgressDuringBlockWriting = hasReportsDuringClose ? 50.0 : 100.0;
 
  int iterLeft = 0;
  int iterTop = 0;
  int iterCols = 0;
  int iterRows = 0;
  QgsRectangle blockExtent;
  std::unique_ptr<QgsRasterBlock> outputBlock;
  while ( outputIter.readNextRasterPart( 1, iterCols, iterRows, outputBlock, iterLeft, iterTop, &blockExtent ) )
  {
    std::vector<std::unique_ptr<QgsRasterBlock>> inputBlocks;
    for ( const QgsRasterAnalysisUtils::RasterLogicInput &i : inputs )
    {
      for ( const int band : i.bands )
      {
        std::unique_ptr<QgsRasterBlock> b( i.interface->block( band, blockExtent, iterCols, iterRows ) );
        inputBlocks.emplace_back( std::move( b ) );
      }
    }
 
    feedback->setProgress( maxProgressDuringBlockWriting * outputIter.progress( 1 ) );
    for ( int row = 0; row < iterRows; row++ )
    {
      if ( feedback->isCanceled() )
        break;
 
      for ( int column = 0; column < iterCols; column++ )
      {
        bool res = false;
        bool resIsNoData = false;
        applyLogicFunc( inputBlocks, res, resIsNoData, row, column, treatNoDataAsFalse );
        if ( resIsNoData )
          noDataCount++;
        else if ( res )
          trueCount++;
        else
          falseCount++;
 
        outputBlock->setValue( row, column, resIsNoData ? outputNoDataValue : ( res ? 1 : 0 ) );
      }
    }
    if ( !destinationRaster->writeBlock( outputBlock.get(), 1, iterLeft, iterTop ) )
    {
      throw QgsProcessingException( QObject::tr( "Could not write raster block: %1" ).arg( destinationRaster->error().summary() ) );
    }
  }
  destinationRaster->setEditable( false );
 
  if ( hasReportsDuringClose )
  {
    std::unique_ptr<QgsFeedback> scaledFeedback( QgsFeedback::createScaledFeedback( feedback, maxProgressDuringBlockWriting, 100.0 ) );
    if ( !destinationRaster->closeWithProgress( scaledFeedback.get() ) )
    {
      if ( feedback->isCanceled() )
        return;
      throw QgsProcessingException( QObject::tr( "Could not write raster dataset" ) );
    }
  }
}
 
std::vector<double> QgsRasterAnalysisUtils::getCellValuesFromBlockStack( const std::vector<std::unique_ptr<QgsRasterBlock>> &inputBlocks, int &row, int &col, bool &noDataInStack )
{
  //get all values from inputBlocks
  std::vector<double> cellValues;
  bool hasNoData = false;
  cellValues.reserve( inputBlocks.size() );
 
  for ( auto &block : inputBlocks )
  {
    double value = 0;
    if ( !block || !block->isValid() )
    {
      noDataInStack = true;
      break;
    }
    else
    {
      value = block->valueAndNoData( row, col, hasNoData );
      if ( hasNoData )
      {
        noDataInStack = true;
        continue; //NoData is not included in the cell value vector
      }
      else
      {
        cellValues.push_back( value );
      }
    }
  }
  return cellValues;
}
 
double QgsRasterAnalysisUtils::meanFromCellValues( std::vector<double> &cellValues, int stackSize )
{
  const double sum = std::accumulate( cellValues.begin(), cellValues.end(), 0.0 );
  const double mean = sum / static_cast<double>( stackSize );
  return mean;
}
 
double QgsRasterAnalysisUtils::medianFromCellValues( std::vector<double> &cellValues, int stackSize )
{
  std::sort( cellValues.begin(), cellValues.end() );
  const bool even = ( stackSize % 2 ) < 1;
  if ( even )
  {
    return ( cellValues[stackSize / 2 - 1] + cellValues[stackSize / 2] ) / 2.0;
  }
  else //odd
  {
    return cellValues[( stackSize + 1 ) / 2 - 1];
  }
}
 
 
double QgsRasterAnalysisUtils::stddevFromCellValues( std::vector<double> &cellValues, int stackSize )
{
  const double variance = varianceFromCellValues( cellValues, stackSize );
  const double stddev = std::sqrt( variance );
  return stddev;
}
 
double QgsRasterAnalysisUtils::varianceFromCellValues( std::vector<double> &cellValues, int stackSize )
{
  const double mean = meanFromCellValues( cellValues, stackSize );
  double accum = 0.0;
  for ( int i = 0; i < stackSize; i++ )
  {
    accum += std::pow( ( cellValues.at( i ) - mean ), 2.0 );
  }
  const double variance = accum / static_cast<double>( stackSize );
  return variance;
}
 
double QgsRasterAnalysisUtils::maximumFromCellValues( std::vector<double> &cellValues )
{
  return *std::max_element( cellValues.begin(), cellValues.end() );
}
 
double QgsRasterAnalysisUtils::minimumFromCellValues( std::vector<double> &cellValues )
{
  return *std::min_element( cellValues.begin(), cellValues.end() );
}
 
double QgsRasterAnalysisUtils::majorityFromCellValues( std::vector<double> &cellValues, const double noDataValue, int stackSize )
{
  if ( stackSize == 1 )
  {
    //output will be same as input if only one layer is entered
    return cellValues[0];
  }
  else if ( stackSize == 2 )
  {
    //if only two layers are input, return NoData if values are not the same (eg. no Majority could  be found)
    return ( qgsDoubleNear( cellValues[0], cellValues[1] ) ) ? cellValues[0] : noDataValue;
  }
  else if ( std::adjacent_find( cellValues.begin(), cellValues.end(), std::not_equal_to<double>() ) == cellValues.end() )
  {
    //check if all values in cellValues are equal
    //output will be same as input if all cellValues of the stack are the same
    return cellValues[0];
  }
  else
  {
    //search for majority using hash map [O(n)]
    std::unordered_map<double, int> map;
 
    for ( int i = 0; i < stackSize; i++ )
    {
      map[cellValues[i]]++;
    }
 
    int maxCount = 0;
    bool multipleMajorities = false;
    double result = noDataValue;
    for ( const auto &pair : std::as_const( map ) )
    {
      if ( maxCount < pair.second )
      {
        result = pair.first;
        maxCount = pair.second;
        multipleMajorities = false;
      }
      else if ( maxCount == pair.second )
      {
        multipleMajorities = true;
      }
    }
    return multipleMajorities ? noDataValue : result;
  }
}
 
double QgsRasterAnalysisUtils::minorityFromCellValues( std::vector<double> &cellValues, const double noDataValue, int stackSize )
{
  if ( stackSize == 1 )
  {
    //output will be same as input if only one layer is entered
    return cellValues[0];
  }
  else if ( stackSize == 2 )
  {
    //if only two layers are input, return NoData if values are not the same (eg. no minority could  be found)
    return ( qgsDoubleNear( cellValues[0], cellValues[1] ) ) ? cellValues[0] : noDataValue;
  }
  else if ( std::adjacent_find( cellValues.begin(), cellValues.end(), std::not_equal_to<double>() ) == cellValues.end() )
  {
    //check if all values in cellValues are equal
    //output will be same as input if all cellValues of the stack are the same
    return cellValues[0];
  }
  else
  {
    //search for minority using hash map [O(n)]
    std::unordered_map<double, int> map;
 
    for ( int i = 0; i < stackSize; i++ )
    {
      map[cellValues[i]]++;
    }
 
    int minCount = stackSize;
    bool multipleMinorities = false;
    double result = noDataValue; //result will stay NoData if no minority value exists
    for ( const auto &pair : std::as_const( map ) )
    {
      if ( minCount > pair.second )
      {
        result = pair.first;
        minCount = pair.second;
        multipleMinorities = false;
      }
      else if ( minCount == pair.second )
      {
        multipleMinorities = true;
      }
    }
    return multipleMinorities ? noDataValue : result;
  }
}
 
double QgsRasterAnalysisUtils::rangeFromCellValues( std::vector<double> &cellValues )
{
  const double max = *std::max_element( cellValues.begin(), cellValues.end() );
  const double min = *std::min_element( cellValues.begin(), cellValues.end() );
  return max - min;
}
 
double QgsRasterAnalysisUtils::varietyFromCellValues( std::vector<double> &cellValues )
{
  const std::unordered_set<double> uniqueValues( cellValues.begin(), cellValues.end() );
  return uniqueValues.size();
}
 
double QgsRasterAnalysisUtils::nearestRankPercentile( std::vector<double> &cellValues, int stackSize, double percentile )
{
  //if percentile equals 0 -> pick the first element of the ordered list
  std::sort( cellValues.begin(), cellValues.end() );
 
  int i = 0;
  if ( percentile > 0 )
  {
    i = std::ceil( percentile * static_cast<double>( stackSize ) ) - 1;
  }
 
  return cellValues[i];
}
 
double QgsRasterAnalysisUtils::interpolatedPercentileInc( std::vector<double> &cellValues, int stackSize, double percentile )
{
  std::sort( cellValues.begin(), cellValues.end() );
 
  if ( qgsDoubleNear( percentile, 1.0 ) )
  {
    return cellValues[stackSize - 1];
  }
  else if ( qgsDoubleNear( percentile, 0.0 ) )
  {
    return cellValues[0];
  }
 
  const double x = ( percentile * ( stackSize - 1 ) );
 
  const int i = static_cast<int>( std::floor( x ) );
  const double xFraction = std::fmod( x, 1 );
 
  if ( stackSize == 1 )
  {
    return cellValues[0];
  }
  else if ( stackSize == 2 )
  {
    return cellValues[0] + ( cellValues[1] - cellValues[0] ) * percentile;
  }
  else
  {
    return cellValues[i] + ( cellValues[i + 1] - cellValues[i] ) * xFraction;
  }
}
 
double QgsRasterAnalysisUtils::interpolatedPercentileExc( std::vector<double> &cellValues, int stackSize, double percentile, double noDataValue )
{
  std::sort( cellValues.begin(), cellValues.end() );
  const double x = ( percentile * ( stackSize + 1 ) );
 
  const int i = static_cast<int>( std::floor( x ) ) - 1;
  const double xFraction = std::fmod( x, 1 );
  const double lowerExcValue = 1.0 / ( static_cast<double>( stackSize ) + 1.0 );
  const double upperExcValue = static_cast<double>( stackSize ) / ( static_cast<double>( stackSize ) + 1.0 );
 
  if ( stackSize < 2 || ( ( percentile < lowerExcValue || percentile > upperExcValue ) ) )
  {
    return noDataValue;
  }
  else
  {
    return cellValues[i] + ( cellValues[i + 1] - cellValues[i] ) * xFraction;
  }
}
 
double QgsRasterAnalysisUtils::interpolatedPercentRankInc( std::vector<double> &cellValues, int stackSize, double value, double noDataValue )
{
  std::sort( cellValues.begin(), cellValues.end() );
 
  if ( value < cellValues[0] || value > cellValues[stackSize - 1] )
  {
    return noDataValue;
  }
  else
  {
    for ( int i = 0; i < stackSize - 1; i++ )
    {
      if ( cellValues[i] <= value && cellValues[i + 1] >= value )
      {
        double fraction = 0.0;
 
        //make sure that next number in the distribution is not the same to prevent NaN fractions
        if ( !qgsDoubleNear( cellValues[i], cellValues[i + 1] ) )
          fraction = ( value - cellValues[i] ) / ( cellValues[i + 1] - cellValues[i] );
 
        return ( fraction + i ) / ( stackSize - 1 );
      }
    }
    return noDataValue;
  }
}
 
double QgsRasterAnalysisUtils::interpolatedPercentRankExc( std::vector<double> &cellValues, int stackSize, double value, double noDataValue )
{
  std::sort( cellValues.begin(), cellValues.end() );
 
  if ( value < cellValues[0] || value > cellValues[stackSize - 1] )
  {
    return noDataValue;
  }
  else
  {
    for ( int i = 0; i < stackSize - 1; i++ )
    {
      if ( cellValues[i] <= value && cellValues[i + 1] >= value )
      {
        double fraction = 0.0;
 
        //make sure that next number in the distribution is not the same to prevent NaN fractions
        if ( !qgsDoubleNear( cellValues[i], cellValues[i + 1] ) )
          fraction = ( value - cellValues[i] ) / ( cellValues[i + 1] - cellValues[i] );
 
        return ( ( i + 1 ) + fraction ) / ( stackSize + 1 );
      }
    }
    return noDataValue;
  }
}
 
double QgsRasterAnalysisUtils::neighborCellDistance( int direction, double cellSizeX, double cellSizeY )
{
  switch ( direction )
  {
    case 0: // north
    case 4: // south
      return cellSizeY;
    case 2: // east
    case 6: // west
      return cellSizeX;
    default: // diagonal
      return std::hypot( cellSizeX, cellSizeY );
  }
  BUILTIN_UNREACHABLE
}
 
bool QgsRasterAnalysisUtils::neighborCellCoordinates( int direction, int row, int column, int &neighborRow, int &neighborColumn, int rows, int columns )
{
  static constexpr std::array<int, 8> COL_DIRECTION_OFFSETS { 0, 1, 1, 1, 0, -1, -1, -1 };
  static constexpr std::array<int, 8> ROW_DIRECTION_OFFSETS { -1, -1, 0, 1, 1, 1, 0, -1 };
 
  neighborColumn = column + COL_DIRECTION_OFFSETS[direction];
  neighborRow = row + ROW_DIRECTION_OFFSETS[direction];
  return ( neighborColumn >= 0 && neighborColumn < columns && neighborRow >= 0 && neighborRow < rows );
}
 
int QgsRasterAnalysisUtils::steepestGradientDirection( const QgsRasterBlock *demBlock, int row, int column, double cellSizeX, double cellSizeY, bool down, bool noEdges )
{
  const int rows = demBlock->height();
  const int columns = demBlock->width();
 
  bool isCenterNoData = false;
  const double z = demBlock->valueAndNoData( row, column, isCenterNoData );
  if ( isCenterNoData )
    return -3;
 
  double maxGradient = 0.0;
  int steepestDir = -1;
 
  bool isNeighborNodata = false;
  bool isFlat = true;
  for ( int dir = 0; dir < 8; ++dir )
  {
    int nCol = 0;
    int nRow = 0;
    if ( neighborCellCoordinates( dir, row, column, nRow, nCol, rows, columns ) )
    {
      const double nZ = demBlock->valueAndNoData( nRow, nCol, isNeighborNodata );
      if ( !isNeighborNodata )
      {
        const double grad = ( z - nZ ) / neighborCellDistance( dir, cellSizeX, cellSizeY );
        isFlat &= qgsDoubleNear( grad, 0 );
        if ( ( !down || grad > 0 ) && ( steepestDir == -1 || grad > maxGradient ) )
        {
          maxGradient = grad;
          steepestDir = dir;
        }
      }
      else if ( noEdges )
      {
        return -3;
      }
    }
    else if ( noEdges )
    {
      return -3;
    }
  }
  return isFlat ? -2 : steepestDir;
}