qgszonalstatistics.cpp

Go to the documentation of this file.

/***************************************************************************
                          qgszonalstatistics.cpp  -  description
                          ----------------------------
    begin                : August 29th, 2009
    copyright            : (C) 2009 by Marco Hugentobler
    email                : marco at hugis dot net
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   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 "qgszonalstatistics.h"
 
#include "processing/qgsrasteranalysisutils.h"
#include "qgsfeatureiterator.h"
#include "qgsfeedback.h"
#include "qgsgeometry.h"
#include "qgsmessagelog.h"
#include "qgsproject.h"
#include "qgsrasterlayer.h"
#include "qgsvectordataprovider.h"
#include "qgsvectorlayer.h"
 
#include <QFile>
#include <QString>
 
using namespace Qt::StringLiterals;
 
QgsZonalStatistics::QgsZonalStatistics( QgsVectorLayer *polygonLayer, QgsRasterLayer *rasterLayer, const QString &attributePrefix, int rasterBand, Qgis::ZonalStatistics stats, bool overwriteExistingFields )
  : QgsZonalStatistics(
      polygonLayer,
      rasterLayer ? rasterLayer->dataProvider() : nullptr,
      rasterLayer ? rasterLayer->crs() : QgsCoordinateReferenceSystem(),
      rasterLayer ? rasterLayer->rasterUnitsPerPixelX() : 0,
      rasterLayer ? rasterLayer->rasterUnitsPerPixelY() : 0,
      attributePrefix,
      rasterBand,
      stats,
      overwriteExistingFields
    )
{}
 
QgsZonalStatistics::QgsZonalStatistics(
  QgsVectorLayer *polygonLayer,
  QgsRasterInterface *rasterInterface,
  const QgsCoordinateReferenceSystem &rasterCrs,
  double rasterUnitsPerPixelX,
  double rasterUnitsPerPixelY,
  const QString &attributePrefix,
  int rasterBand,
  Qgis::ZonalStatistics stats,
  bool overwriteExistingFields
)
  : mRasterInterface( rasterInterface )
  , mRasterCrs( rasterCrs )
  , mCellSizeX( std::fabs( rasterUnitsPerPixelX ) )
  , mCellSizeY( std::fabs( rasterUnitsPerPixelY ) )
  , mRasterBand( rasterBand )
  , mPolygonLayer( polygonLayer )
  , mAttributePrefix( attributePrefix )
  , mStatistics( stats )
  , mOverwriteExistingFields( overwriteExistingFields )
{}
 
Qgis::ZonalStatisticResult QgsZonalStatistics::calculateStatistics( QgsFeedback *feedback )
{
  if ( !mRasterInterface )
  {
    return Qgis::ZonalStatisticResult::RasterInvalid;
  }
 
  if ( mRasterInterface->bandCount() < mRasterBand )
  {
    return Qgis::ZonalStatisticResult::RasterBandInvalid;
  }
 
  if ( !mPolygonLayer || mPolygonLayer->geometryType() != Qgis::GeometryType::Polygon )
  {
    return Qgis::ZonalStatisticResult::LayerTypeWrong;
  }
 
  QgsVectorDataProvider *vectorProvider = mPolygonLayer->dataProvider();
  if ( !vectorProvider )
  {
    return Qgis::ZonalStatisticResult::LayerInvalid;
  }
 
  QMap<Qgis::ZonalStatistic, int> statFieldIndexes;
 
  //add the new fields to the provider
  int oldFieldCount = vectorProvider->fields().count();
  QList<QgsField> newFieldList;
  for ( Qgis::ZonalStatistic stat :
        { Qgis::ZonalStatistic::Count,
          Qgis::ZonalStatistic::Sum,
          Qgis::ZonalStatistic::Mean,
          Qgis::ZonalStatistic::Median,
          Qgis::ZonalStatistic::StDev,
          Qgis::ZonalStatistic::Min,
          Qgis::ZonalStatistic::Max,
          Qgis::ZonalStatistic::Range,
          Qgis::ZonalStatistic::Minority,
          Qgis::ZonalStatistic::Majority,
          Qgis::ZonalStatistic::Variety,
          Qgis::ZonalStatistic::Variance } )
  {
    if ( mStatistics & stat )
    {
      QString targetFieldName = mAttributePrefix + QgsZonalStatistics::shortName( stat );
      int targetIndex = vectorProvider->fields().indexFromName( targetFieldName );
 
      if ( mOverwriteExistingFields && targetIndex != -1 )
      {
        statFieldIndexes.insert( stat, targetIndex );
      }
      else
      {
        if ( mOverwriteExistingFields && targetIndex == -1 )
        {
          QgsMessageLog::logMessage( QObject::tr( "Could not find field '%1'. A new field will be created instead." ).arg( targetFieldName ), QObject::tr( "Zonal Statistics" ), Qgis::MessageLevel::Warning );
        }
 
        QString fieldName = getUniqueFieldName( targetFieldName, newFieldList );
        QgsField field( fieldName, QMetaType::Type::Double, u"double precision"_s );
        newFieldList.push_back( field );
        statFieldIndexes.insert( stat, oldFieldCount + newFieldList.count() - 1 );
      }
    }
  }
 
  vectorProvider->addAttributes( newFieldList );
 
  long featureCount = vectorProvider->featureCount();
 
  QgsFeatureRequest request;
  request.setNoAttributes();
 
  request.setDestinationCrs( mRasterCrs, QgsProject::instance()->transformContext() );
  QgsFeatureIterator fi = vectorProvider->getFeatures( request );
  QgsFeature feature;
 
  int featureCounter = 0;
 
  QgsChangedAttributesMap changeMap;
  while ( fi.nextFeature( feature ) )
  {
    ++featureCounter;
    if ( feedback && feedback->isCanceled() )
    {
      break;
    }
 
    if ( feedback )
    {
      feedback->setProgress( 100.0 * static_cast<double>( featureCounter ) / featureCount );
    }
 
    QgsGeometry featureGeometry = feature.geometry();
 
    QMap<Qgis::ZonalStatistic, QVariant> results = calculateStatistics( mRasterInterface, featureGeometry, mCellSizeX, mCellSizeY, mRasterBand, mStatistics );
 
    if ( results.empty() )
      continue;
 
    QgsAttributeMap changeAttributeMap;
    for ( const auto &result : results.toStdMap() )
    {
      changeAttributeMap.insert( statFieldIndexes.value( result.first ), result.second );
    }
 
    changeMap.insert( feature.id(), changeAttributeMap );
  }
 
  vectorProvider->changeAttributeValues( changeMap );
  mPolygonLayer->updateFields();
 
  if ( feedback )
  {
    if ( feedback->isCanceled() )
      return Qgis::ZonalStatisticResult::Canceled;
 
    feedback->setProgress( 100 );
  }
 
  return Qgis::ZonalStatisticResult::Success;
}
 
QString QgsZonalStatistics::getUniqueFieldName( const QString &fieldName, const QList<QgsField> &newFields )
{
  QgsVectorDataProvider *provider = mPolygonLayer->dataProvider();
  const bool isShapefile = provider->storageType().contains( "ESRI Shapefile"_L1 );
 
  QList<QgsField> allFields = provider->fields().toList();
  allFields.append( newFields );
 
  QString baseName = isShapefile ? fieldName.mid( 0, 10 ) : fieldName;
 
  auto fieldExists = [&allFields]( const QString &name ) -> bool {
    for ( const QgsField &field : std::as_const( allFields ) )
    {
      if ( field.name().compare( name, Qt::CaseSensitivity::CaseSensitive ) == 0 )
      {
        return true;
      }
    }
    return false;
  };
 
  if ( !fieldExists( baseName ) )
  {
    return baseName;
  }
 
  int n = 1;
  QString newName;
  while ( true )
  {
    if ( isShapefile )
    {
      newName = u"%1_%2"_s.arg( fieldName.mid( 0, n < 9 ? 8 : 7 ) ).arg( n );
    }
    else
    {
      newName = u"%1_%2"_s.arg( fieldName ).arg( n );
    }
 
    if ( !fieldExists( newName ) )
      return newName;
 
    n++;
  }
}
 
QString QgsZonalStatistics::displayName( Qgis::ZonalStatistic statistic )
{
  switch ( statistic )
  {
    case Qgis::ZonalStatistic::Count:
      return QObject::tr( "Count" );
    case Qgis::ZonalStatistic::Sum:
      return QObject::tr( "Sum" );
    case Qgis::ZonalStatistic::Mean:
      return QObject::tr( "Mean" );
    case Qgis::ZonalStatistic::Median:
      return QObject::tr( "Median" );
    case Qgis::ZonalStatistic::StDev:
      return QObject::tr( "St dev" );
    case Qgis::ZonalStatistic::Min:
      return QObject::tr( "Minimum" );
    case Qgis::ZonalStatistic::Max:
      return QObject::tr( "Maximum" );
    case Qgis::ZonalStatistic::MinimumPoint:
      return QObject::tr( "Minimum point" );
    case Qgis::ZonalStatistic::MaximumPoint:
      return QObject::tr( "Maximum point" );
    case Qgis::ZonalStatistic::Range:
      return QObject::tr( "Range" );
    case Qgis::ZonalStatistic::Minority:
      return QObject::tr( "Minority" );
    case Qgis::ZonalStatistic::Majority:
      return QObject::tr( "Majority" );
    case Qgis::ZonalStatistic::Variety:
      return QObject::tr( "Variety" );
    case Qgis::ZonalStatistic::Variance:
      return QObject::tr( "Variance" );
    case Qgis::ZonalStatistic::All:
    case Qgis::ZonalStatistic::Default:
      return QString();
  }
  return QString();
}
 
QString QgsZonalStatistics::shortName( Qgis::ZonalStatistic statistic )
{
  switch ( statistic )
  {
    case Qgis::ZonalStatistic::Count:
      return u"count"_s;
    case Qgis::ZonalStatistic::Sum:
      return u"sum"_s;
    case Qgis::ZonalStatistic::Mean:
      return u"mean"_s;
    case Qgis::ZonalStatistic::Median:
      return u"median"_s;
    case Qgis::ZonalStatistic::StDev:
      return u"stdev"_s;
    case Qgis::ZonalStatistic::Min:
      return u"min"_s;
    case Qgis::ZonalStatistic::Max:
      return u"max"_s;
    case Qgis::ZonalStatistic::MinimumPoint:
      return u"minpoint"_s;
    case Qgis::ZonalStatistic::MaximumPoint:
      return u"maxpoint"_s;
    case Qgis::ZonalStatistic::Range:
      return u"range"_s;
    case Qgis::ZonalStatistic::Minority:
      return u"minority"_s;
    case Qgis::ZonalStatistic::Majority:
      return u"majority"_s;
    case Qgis::ZonalStatistic::Variety:
      return u"variety"_s;
    case Qgis::ZonalStatistic::Variance:
      return u"variance"_s;
    case Qgis::ZonalStatistic::All:
    case Qgis::ZonalStatistic::Default:
      return QString();
  }
  return QString();
}

QMap<int, QVariant> QgsZonalStatistics::calculateStatisticsInt(
  QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, double cellSizeX, double cellSizeY, int rasterBand, Qgis::ZonalStatistics statistics
)
{
  const auto result { QgsZonalStatistics::calculateStatistics( rasterInterface, geometry, cellSizeX, cellSizeY, rasterBand, statistics ) };
  QMap<int, QVariant> pyResult;
  for ( auto it = result.constBegin(); it != result.constEnd(); ++it )
  {
    pyResult.insert( static_cast<int>( it.key() ), it.value() );
  }
  return pyResult;
}
 
QMap<Qgis::ZonalStatistic, QVariant> QgsZonalStatistics::calculateStatistics(
  QgsRasterInterface *rasterInterface, const QgsGeometry &geometry, double cellSizeX, double cellSizeY, int rasterBand, Qgis::ZonalStatistics statistics
)
{
  QMap<Qgis::ZonalStatistic, QVariant> results;
 
  if ( geometry.isEmpty() )
    return results;
 
  QgsRectangle rasterBBox = rasterInterface->extent();
 
  QgsRectangle featureRect = geometry.boundingBox().intersect( rasterBBox );
 
  if ( featureRect.isEmpty() )
    return results;
 
  bool statsStoreValues = ( statistics & Qgis::ZonalStatistic::Median ) || ( statistics & Qgis::ZonalStatistic::StDev ) || ( statistics & Qgis::ZonalStatistic::Variance );
  bool statsStoreValueCount = ( statistics & Qgis::ZonalStatistic::Minority ) || ( statistics & Qgis::ZonalStatistic::Majority );
 
  FeatureStats featureStats( statsStoreValues, statsStoreValueCount );
 
  int nCellsXProvider = rasterInterface->xSize();
  int nCellsYProvider = rasterInterface->ySize();
 
  int nCellsX, nCellsY;
  QgsRectangle rasterBlockExtent;
  QgsRasterAnalysisUtils::cellInfoForBBox( rasterBBox, featureRect, cellSizeX, cellSizeY, nCellsX, nCellsY, nCellsXProvider, nCellsYProvider, rasterBlockExtent );
 
  featureStats.reset();
  QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [&featureStats]( double value, const QgsPointXY &point ) {
    featureStats.addValue( value, point );
  } );
 
  if ( featureStats.count <= 1 )
  {
    //the cell resolution is probably larger than the polygon area. We switch to precise pixel - polygon intersection in this case
    featureStats.reset();
    QgsRasterAnalysisUtils::
      statisticsFromPreciseIntersection( rasterInterface, rasterBand, geometry, nCellsX, nCellsY, cellSizeX, cellSizeY, rasterBlockExtent, [&featureStats]( double value, double weight, const QgsPointXY &point ) {
        featureStats.addValue( value, point, weight );
      } );
  }
 
  // calculate the statistics
 
  if ( statistics & Qgis::ZonalStatistic::Count )
    results.insert( Qgis::ZonalStatistic::Count, QVariant( featureStats.count ) );
  if ( statistics & Qgis::ZonalStatistic::Sum )
    results.insert( Qgis::ZonalStatistic::Sum, QVariant( featureStats.sum ) );
  if ( featureStats.count > 0 )
  {
    double mean = featureStats.sum / featureStats.count;
    if ( statistics & Qgis::ZonalStatistic::Mean )
      results.insert( Qgis::ZonalStatistic::Mean, QVariant( mean ) );
    if ( statistics & Qgis::ZonalStatistic::Median )
    {
      std::sort( featureStats.values.begin(), featureStats.values.end() );
      int size = featureStats.values.count();
      bool even = ( size % 2 ) < 1;
      double medianValue;
      if ( even )
      {
        medianValue = ( featureStats.values.at( size / 2 - 1 ) + featureStats.values.at( size / 2 ) ) / 2;
      }
      else //odd
      {
        medianValue = featureStats.values.at( ( size + 1 ) / 2 - 1 );
      }
      results.insert( Qgis::ZonalStatistic::Median, QVariant( medianValue ) );
    }
    if ( statistics & Qgis::ZonalStatistic::StDev || statistics & Qgis::ZonalStatistic::Variance )
    {
      double sumSquared = 0;
      for ( int i = 0; i < featureStats.values.count(); ++i )
      {
        double diff = featureStats.values.at( i ) - mean;
        sumSquared += diff * diff;
      }
      double variance = sumSquared / featureStats.values.count();
      if ( statistics & Qgis::ZonalStatistic::StDev )
      {
        double stdev = std::pow( variance, 0.5 );
        results.insert( Qgis::ZonalStatistic::StDev, QVariant( stdev ) );
      }
      if ( statistics & Qgis::ZonalStatistic::Variance )
        results.insert( Qgis::ZonalStatistic::Variance, QVariant( variance ) );
    }
    if ( statistics & Qgis::ZonalStatistic::Min )
      results.insert( Qgis::ZonalStatistic::Min, QVariant( featureStats.min ) );
    if ( statistics & Qgis::ZonalStatistic::Max )
      results.insert( Qgis::ZonalStatistic::Max, QVariant( featureStats.max ) );
    if ( statistics & Qgis::ZonalStatistic::MinimumPoint )
      results.insert( Qgis::ZonalStatistic::MinimumPoint, QVariant( featureStats.minPoint ) );
    if ( statistics & Qgis::ZonalStatistic::MaximumPoint )
      results.insert( Qgis::ZonalStatistic::MaximumPoint, QVariant( featureStats.maxPoint ) );
    if ( statistics & Qgis::ZonalStatistic::Range )
      results.insert( Qgis::ZonalStatistic::Range, QVariant( featureStats.max - featureStats.min ) );
    if ( statistics & Qgis::ZonalStatistic::Minority || statistics & Qgis::ZonalStatistic::Majority )
    {
      int64_t lowestCount = std::numeric_limits<int64_t>::max();
      int64_t highestCount = 0;
      double valueWithLowestCount = 0;
      double valueWithHighestCount = 0;
      for ( auto [value, count] : std::as_const( featureStats.valueCount ) )
      {
        if ( count < lowestCount )
        {
          lowestCount = count;
          valueWithLowestCount = value;
        }
        if ( count > highestCount )
        {
          highestCount = count;
          valueWithHighestCount = value;
        }
      }
      if ( statistics & Qgis::ZonalStatistic::Minority )
      {
        results.insert( Qgis::ZonalStatistic::Minority, QVariant( valueWithLowestCount ) );
      }
      if ( statistics & Qgis::ZonalStatistic::Majority )
      {
        results.insert( Qgis::ZonalStatistic::Majority, QVariant( valueWithHighestCount ) );
      }
    }
    if ( statistics & Qgis::ZonalStatistic::Variety )
      results.insert( Qgis::ZonalStatistic::Variety, QVariant( static_cast<qlonglong>( featureStats.valueCount.size() ) ) );
  }
 
  return results;
}