From fdaa57a273a98aab286115d159eafa5e7a6cda09 Mon Sep 17 00:00:00 2001 From: nirvn Date: Tue, 1 May 2018 11:42:31 +0700 Subject: [PATCH] [FEATURE][processing] Zonal histogram algorithm --- .../testdata/custom/circular_strings.gpkg | Bin 131072 -> 131072 bytes .../testdata/expected/zones_histogram.gml | 34 ++ .../testdata/expected/zones_histogram.xsd | 58 ++++ .../tests/testdata/qgis_algorithm_tests.yaml | 16 + .../processing/tests/testdata/zones.gml | 24 ++ .../processing/tests/testdata/zones.xsd | 23 ++ src/analysis/CMakeLists.txt | 1 + .../processing/qgsalgorithmzonalhistogram.cpp | 314 ++++++++++++++++++ .../processing/qgsalgorithmzonalhistogram.h | 77 +++++ .../processing/qgsnativealgorithms.cpp | 2 + 10 files changed, 549 insertions(+) create mode 100644 python/plugins/processing/tests/testdata/expected/zones_histogram.gml create mode 100644 python/plugins/processing/tests/testdata/expected/zones_histogram.xsd create mode 100644 python/plugins/processing/tests/testdata/zones.gml create mode 100644 python/plugins/processing/tests/testdata/zones.xsd create mode 100644 src/analysis/processing/qgsalgorithmzonalhistogram.cpp create mode 100644 src/analysis/processing/qgsalgorithmzonalhistogram.h diff --git a/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg b/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg index d1b1e702ca1da2bd5454a8492bd98651253cada5..601c38a7fec5a6c137b673a63f0ef58080ad6d9f 100644 GIT binary patch delta 26 hcmZo@;Am*zm>|uVKT*b+HJ?FGW|uVJW))dA${{d~02>Ad2 diff --git a/python/plugins/processing/tests/testdata/expected/zones_histogram.gml b/python/plugins/processing/tests/testdata/expected/zones_histogram.gml new file mode 100644 index 00000000000..7067a1369cf --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/zones_histogram.gml @@ -0,0 +1,34 @@ + + + + + 270743.9470328064458959.32367457 + 270848.2436236554459031.472508 + + + + + + 270744.399022365,4459030.56852888 270797.959785163,4459031.472508 270797.959785163,4458978.58972954 270743.947032806,4458979.0417191 270744.399022365,4459030.56852888 + 4 + 6 + 6 + 9 + 0 + + + + + 270818.41231273,4458978.87222301 270818.18631795,4458959.54966935 270848.243623655,4458959.32367457 270848.017628875,4458978.87222301 270818.41231273,4458978.87222301 + 0 + 0 + 0 + 0 + 6 + + + diff --git a/python/plugins/processing/tests/testdata/expected/zones_histogram.xsd b/python/plugins/processing/tests/testdata/expected/zones_histogram.xsd new file mode 100644 index 00000000000..71481a18fe5 --- /dev/null +++ b/python/plugins/processing/tests/testdata/expected/zones_histogram.xsd @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml b/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml index 0d46c3b2b8d..b2eae7eddce 100644 --- a/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml +++ b/python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml @@ -3460,6 +3460,22 @@ tests: geometry: precision: 5 + - algorithm: qgis:zonalhistogram + name: zonal histogram + params: + COLUMN_PREFIX: HISTO_ + INPUT_RASTER: + name: raster.tif + type: raster + INPUT_VECTOR: + name: zones.gml + type: vector + RASTER_BAND: 1 + results: + OUTPUT: + name: expected/zones_histogram.gml + type: vector + - algorithm: native:fixgeometries name: Fix geometries params: diff --git a/python/plugins/processing/tests/testdata/zones.gml b/python/plugins/processing/tests/testdata/zones.gml new file mode 100644 index 00000000000..929c35ad135 --- /dev/null +++ b/python/plugins/processing/tests/testdata/zones.gml @@ -0,0 +1,24 @@ + + + + + 270743.94703280574458959.323674565 + 270848.24362365524459031.472508 + + + + + + 270744.399022365,4459030.56852888 270797.959785163,4459031.472508 270797.959785163,4458978.58972954 270743.947032806,4458979.0417191 270744.399022365,4459030.56852888 + + + + + 270818.41231273,4458978.87222301 270818.18631795,4458959.54966935 270848.243623655,4458959.32367457 270848.017628875,4458978.87222301 270818.41231273,4458978.87222301 + + + diff --git a/python/plugins/processing/tests/testdata/zones.xsd b/python/plugins/processing/tests/testdata/zones.xsd new file mode 100644 index 00000000000..dfdd253975c --- /dev/null +++ b/python/plugins/processing/tests/testdata/zones.xsd @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index aee1374d44a..83d8c35f621 100755 --- a/src/analysis/CMakeLists.txt +++ b/src/analysis/CMakeLists.txt @@ -82,6 +82,7 @@ SET(QGIS_ANALYSIS_SRCS processing/qgsalgorithmunion.cpp processing/qgsalgorithmuniquevalueindex.cpp processing/qgsalgorithmwedgebuffers.cpp + processing/qgsalgorithmzonalhistogram.cpp processing/qgsnativealgorithms.cpp processing/qgsoverlayutils.cpp diff --git a/src/analysis/processing/qgsalgorithmzonalhistogram.cpp b/src/analysis/processing/qgsalgorithmzonalhistogram.cpp new file mode 100644 index 00000000000..e2719c08e8c --- /dev/null +++ b/src/analysis/processing/qgsalgorithmzonalhistogram.cpp @@ -0,0 +1,314 @@ +/*************************************************************************** + 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 "qgsgeos.h" +#include "qgslogger.h" + +///@cond PRIVATE + +QString QgsZonalHistogramAlgorithm::name() const +{ + return QStringLiteral( "zonalhistogram" ); +} + +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 QStringLiteral( "rasteranalysis" ); +} + +void QgsZonalHistogramAlgorithm::initAlgorithm( const QVariantMap & ) +{ + addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ), + QObject::tr( "Raster layer" ) ) ); + addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ), + QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT_RASTER" ) ) ); + + addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT_VECTOR" ), + QObject::tr( "Vector layer containing zones" ), QList< int >() << QgsProcessing::TypeVectorPolygon ) ); + + addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "HISTO_" ), false, true ) ); + + addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Output zones" ), QgsProcessing::TypeVectorPolygon ) ); +} + +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." ); +} + +QgsZonalHistogramAlgorithm *QgsZonalHistogramAlgorithm::createInstance() const +{ + return new QgsZonalHistogramAlgorithm(); +} + +bool QgsZonalHistogramAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback * ) +{ + QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_RASTER" ), context ); + mBand = parameterAsInt( parameters, QStringLiteral( "RASTER_BAND" ), context ); + + if ( !layer ) + throw QgsProcessingException( QObject::tr( "Could not load raster layer" ) ); + + mInterface.reset( layer->dataProvider()->clone() ); + mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( mBand ); + mNodataValue = layer->dataProvider()->sourceNoDataValue( mBand ); + mExtent = layer->extent(); + mCrs = layer->crs(); + mRasterUnitsPerPixelX = std::abs( layer->rasterUnitsPerPixelX() ); + mRasterUnitsPerPixelY = std::abs( layer->rasterUnitsPerPixelX() ); + mNbCellsXProvider = mInterface->xSize(); + mNbCellsYProvider = mInterface->ySize(); + + return true; +} + +QVariantMap QgsZonalHistogramAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) +{ + + std::unique_ptr< QgsFeatureSource > zones( parameterAsSource( parameters, QStringLiteral( "INPUT_VECTOR" ), context ) ); + if ( !zones ) + throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT_VECTOR" ) ) ); + + long count = zones->featureCount(); + 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.setSubsetOfAttributes( QgsAttributeList() ); + if ( zones->sourceCrs() != mCrs ) + { + request.setDestinationCrs( mCrs, context.transformContext() ); + } + QgsFeatureIterator it = zones->getFeatures( request ); + QgsFeature f; + while ( it.nextFeature( f ) ) + { + if ( feedback && feedback->isCanceled() ) + { + break; + } + feedback->setProgress( current * step ); + + if ( !f.hasGeometry() ) + { + current++; + continue; + } + + QgsGeometry featureGeometry = f.geometry(); + QgsRectangle intersectRect = featureGeometry.boundingBox().intersect( &mExtent ); + if ( intersectRect.isEmpty() ) + { + current++; + continue; + } + + int offsetX, offsetY, nCellsX, nCellsY; + // Get offset in pixels in x- and y- direction + offsetX = ( int )( ( intersectRect.xMinimum() - mExtent.xMinimum() ) / mRasterUnitsPerPixelX ); + offsetY = ( int )( ( mExtent.yMaximum() - intersectRect.yMaximum() ) / mRasterUnitsPerPixelY ); + + int maxColumn = ( int )( ( intersectRect.xMaximum() - mExtent.xMinimum() ) / mRasterUnitsPerPixelX ) + 1; + int maxRow = ( int )( ( mExtent.yMaximum() - intersectRect.yMinimum() ) / mRasterUnitsPerPixelY ) + 1; + + nCellsX = maxColumn - offsetX; + nCellsY = maxRow - offsetY; + + // Avoid access to cells outside of the raster (may occur because of rounding) + if ( ( offsetX + nCellsX ) > mNbCellsXProvider ) + { + nCellsX = mNbCellsXProvider - offsetX; + } + if ( ( offsetY + nCellsY ) > mNbCellsYProvider ) + { + nCellsY = mNbCellsYProvider - offsetY; + } + + QHash< double, qgssize > fUniqueValues; + middlePoints( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, fUniqueValues ); + + 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 + preciseIntersection( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, fUniqueValues ); + } + + 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() ); + + QString fieldPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context ); + QgsFields newFields; + for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it ) + { + newFields.append( QgsField( QStringLiteral( "%1%2" ).arg( fieldPrefix ).arg( mHasNoDataValue && *it == mNodataValue ? QStringLiteral( "NODATA" ) : QString::number( *it ) ), QVariant::LongLong, QString(), -1, 0 ) ); + } + QgsFields fields = QgsProcessingUtils::combineFields( zones->fields(), newFields ); + + QString dest; + std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, + zones->wkbType(), zones->sourceCrs() ) ); + if ( !sink ) + throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) ); + + it = zones->getFeatures( QgsFeatureRequest() ); + while ( it.nextFeature( f ) ) + { + QgsAttributes attributes = f.attributes(); + 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 ); + + sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ); + } + + QVariantMap outputs; + outputs.insert( QStringLiteral( "OUTPUT" ), dest ); + return outputs; +} + +void QgsZonalHistogramAlgorithm::middlePoints( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues ) +{ + double cellCenterX, cellCenterY; + + cellCenterY = mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - mRasterUnitsPerPixelY / 2; + + geos::unique_ptr polyGeos( QgsGeos::asGeos( poly ) ); + if ( !polyGeos ) + { + return; + } + + GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler(); + geos::prepared_unique_ptr polyGeosPrepared( GEOSPrepare_r( geosctxt, polyGeos.get() ) ); + if ( !polyGeosPrepared ) + { + return; + } + + GEOSCoordSequence *cellCenterCoords = nullptr; + geos::unique_ptr currentCellCenter; + + QgsRectangle blockRect( mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX, + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - nCellsY * mRasterUnitsPerPixelY, + mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX + nCellsX * mRasterUnitsPerPixelX, + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY ); + std::unique_ptr< QgsRasterBlock > block( mInterface->block( mBand, blockRect, nCellsX, nCellsY ) ); + for ( int i = 0; i < nCellsY; ++i ) + { + cellCenterX = mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX + mRasterUnitsPerPixelX / 2; + for ( int j = 0; j < nCellsX; ++j ) + { + if ( !std::isnan( block->value( i, j ) ) ) + { + cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 ); + GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX ); + GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY ); + currentCellCenter.reset( GEOSGeom_createPoint_r( geosctxt, cellCenterCoords ) ); + if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared.get(), currentCellCenter.get() ) ) + { + uniqueValues[block->value( i, j )]++; + } + } + cellCenterX += mRasterUnitsPerPixelX; + } + cellCenterY -= mRasterUnitsPerPixelY; + } +} + +void QgsZonalHistogramAlgorithm::preciseIntersection( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues ) +{ + double currentY = mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - mRasterUnitsPerPixelY / 2; + QgsGeometry pixelRectGeometry; + + double hCellSizeX = mRasterUnitsPerPixelX / 2.0; + double hCellSizeY = mRasterUnitsPerPixelY / 2.0; + + QgsRectangle blockRect( mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX, + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - nCellsY * mRasterUnitsPerPixelY, + mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX + nCellsX * mRasterUnitsPerPixelX, + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY ); + std::unique_ptr< QgsRasterBlock > block( mInterface->block( mBand, blockRect, nCellsX, nCellsY ) ); + for ( int i = 0; i < nCellsY; ++i ) + { + double currentX = mExtent.xMinimum() + mRasterUnitsPerPixelX / 2.0 + pixelOffsetX * mRasterUnitsPerPixelX; + for ( int j = 0; j < nCellsX; ++j ) + { + if ( !std::isnan( block->value( i, j ) ) ) + { + pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) ); + if ( !pixelRectGeometry.isNull() ) + { + //intersection + QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly ); + if ( !intersectGeometry.isNull() ) + { + double intersectionArea = intersectGeometry.area(); + if ( intersectionArea >= 0.0 ) + { + uniqueValues[block->value( i, j )]++; + } + } + pixelRectGeometry = QgsGeometry(); + } + } + currentX += mRasterUnitsPerPixelX; + } + currentY -= mRasterUnitsPerPixelY; + } +} + +///@endcond + + + diff --git a/src/analysis/processing/qgsalgorithmzonalhistogram.h b/src/analysis/processing/qgsalgorithmzonalhistogram.h new file mode 100644 index 00000000000..d6123611761 --- /dev/null +++ b/src/analysis/processing/qgsalgorithmzonalhistogram.h @@ -0,0 +1,77 @@ +/*************************************************************************** + qgsalgorithmzonalhistogram.h + --------------------- + 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. * + * * + ***************************************************************************/ + +#ifndef QGSALGORITHMZONALHISTOGRAM_H +#define QGSALGORITHMZONALHISTOGRAM_H + +#define SIP_NO_FILE + +#include "qgis.h" +#include "qgsprocessingalgorithm.h" + +///@cond PRIVATE + +/** + * Native zonal histogram algorithm. + */ +class QgsZonalHistogramAlgorithm : public QgsProcessingAlgorithm +{ + + public: + + QgsZonalHistogramAlgorithm() = default; + void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override; + QString name() const override; + QString displayName() const override; + QStringList tags() const override; + QString group() const override; + QString groupId() const override; + QString shortHelpString() const override; + QgsZonalHistogramAlgorithm *createInstance() const override SIP_FACTORY; + + protected: + + bool prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + QVariantMap processAlgorithm( const QVariantMap ¶meters, + QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override; + + //! Fetch unique values by considering the pixels where the center point is within the polygon (fast) + void middlePoints( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues ); + + //! Fetch unique values with precise pixel - polygon intersection test (slow) + void preciseIntersection( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues ); + + private: + + std::unique_ptr< QgsRasterInterface > mInterface; + int mBand; + bool mHasNoDataValue = false; + float mNodataValue = -1; + QgsRectangle mExtent; + QgsCoordinateReferenceSystem mCrs; + double mRasterUnitsPerPixelX; + double mRasterUnitsPerPixelY; + double mNbCellsXProvider; + double mNbCellsYProvider; + +}; + +///@endcond PRIVATE + +#endif // QGSALGORITHMZONALHISTOGRAM_H + + diff --git a/src/analysis/processing/qgsnativealgorithms.cpp b/src/analysis/processing/qgsnativealgorithms.cpp index 6bdb20eab6b..c7e02d8f1e6 100644 --- a/src/analysis/processing/qgsnativealgorithms.cpp +++ b/src/analysis/processing/qgsnativealgorithms.cpp @@ -79,6 +79,7 @@ #include "qgsalgorithmunion.h" #include "qgsalgorithmuniquevalueindex.h" #include "qgsalgorithmwedgebuffers.h" +#include "qgsalgorithmzonalhistogram.h" ///@cond PRIVATE @@ -186,6 +187,7 @@ void QgsNativeAlgorithms::loadAlgorithms() addAlgorithm( new QgsUnionAlgorithm() ); addAlgorithm( new QgsVariableWidthBufferByMAlgorithm() ); addAlgorithm( new QgsWedgeBuffersAlgorithm() ); + addAlgorithm( new QgsZonalHistogramAlgorithm() ); }