|
| 1 | +/*************************************************************************** |
| 2 | + qgsalgorithmzonalhistogram.cpp |
| 3 | + --------------------- |
| 4 | + begin : May, 2018 |
| 5 | + copyright : (C) 2018 by Mathieu Pellerin |
| 6 | + email : nirvn dot asia at gmail dot com |
| 7 | + ***************************************************************************/ |
| 8 | + |
| 9 | +/*************************************************************************** |
| 10 | + * * |
| 11 | + * This program is free software; you can redistribute it and/or modify * |
| 12 | + * it under the terms of the GNU General Public License as published by * |
| 13 | + * the Free Software Foundation; either version 2 of the License, or * |
| 14 | + * (at your option) any later version. * |
| 15 | + * * |
| 16 | + ***************************************************************************/ |
| 17 | + |
| 18 | +#include "qgsalgorithmzonalhistogram.h" |
| 19 | +#include "qgsgeos.h" |
| 20 | +#include "qgslogger.h" |
| 21 | + |
| 22 | +///@cond PRIVATE |
| 23 | + |
| 24 | +QString QgsZonalHistogramAlgorithm::name() const |
| 25 | +{ |
| 26 | + return QStringLiteral( "zonalhistogram" ); |
| 27 | +} |
| 28 | + |
| 29 | +QString QgsZonalHistogramAlgorithm::displayName() const |
| 30 | +{ |
| 31 | + return QObject::tr( "Zonal histogram" ); |
| 32 | +} |
| 33 | + |
| 34 | +QStringList QgsZonalHistogramAlgorithm::tags() const |
| 35 | +{ |
| 36 | + return QObject::tr( "raster,unique,values,count,area,statistics" ).split( ',' ); |
| 37 | +} |
| 38 | + |
| 39 | +QString QgsZonalHistogramAlgorithm::group() const |
| 40 | +{ |
| 41 | + return QObject::tr( "Raster analysis" ); |
| 42 | +} |
| 43 | + |
| 44 | +QString QgsZonalHistogramAlgorithm::groupId() const |
| 45 | +{ |
| 46 | + return QStringLiteral( "rasteranalysis" ); |
| 47 | +} |
| 48 | + |
| 49 | +void QgsZonalHistogramAlgorithm::initAlgorithm( const QVariantMap & ) |
| 50 | +{ |
| 51 | + addParameter( new QgsProcessingParameterRasterLayer( QStringLiteral( "INPUT_RASTER" ), |
| 52 | + QObject::tr( "Raster layer" ) ) ); |
| 53 | + addParameter( new QgsProcessingParameterBand( QStringLiteral( "RASTER_BAND" ), |
| 54 | + QObject::tr( "Band number" ), 1, QStringLiteral( "INPUT_RASTER" ) ) ); |
| 55 | + |
| 56 | + addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT_VECTOR" ), |
| 57 | + QObject::tr( "Vector layer containing zones" ), QList< int >() << QgsProcessing::TypeVectorPolygon ) ); |
| 58 | + |
| 59 | + addParameter( new QgsProcessingParameterString( QStringLiteral( "COLUMN_PREFIX" ), QObject::tr( "Output column prefix" ), QStringLiteral( "HISTO_" ), false, true ) ); |
| 60 | + |
| 61 | + addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Output zones" ), QgsProcessing::TypeVectorPolygon ) ); |
| 62 | +} |
| 63 | + |
| 64 | +QString QgsZonalHistogramAlgorithm::shortHelpString() const |
| 65 | +{ |
| 66 | + return QObject::tr( "This algorithm appends fields representing counts of each unique value from a raster layer contained within zones defined as polygons." ); |
| 67 | +} |
| 68 | + |
| 69 | +QgsZonalHistogramAlgorithm *QgsZonalHistogramAlgorithm::createInstance() const |
| 70 | +{ |
| 71 | + return new QgsZonalHistogramAlgorithm(); |
| 72 | +} |
| 73 | + |
| 74 | +bool QgsZonalHistogramAlgorithm::prepareAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback * ) |
| 75 | +{ |
| 76 | + QgsRasterLayer *layer = parameterAsRasterLayer( parameters, QStringLiteral( "INPUT_RASTER" ), context ); |
| 77 | + mBand = parameterAsInt( parameters, QStringLiteral( "RASTER_BAND" ), context ); |
| 78 | + |
| 79 | + if ( !layer ) |
| 80 | + throw QgsProcessingException( QObject::tr( "Could not load raster layer" ) ); |
| 81 | + |
| 82 | + mInterface.reset( layer->dataProvider()->clone() ); |
| 83 | + mHasNoDataValue = layer->dataProvider()->sourceHasNoDataValue( mBand ); |
| 84 | + mNodataValue = layer->dataProvider()->sourceNoDataValue( mBand ); |
| 85 | + mExtent = layer->extent(); |
| 86 | + mCrs = layer->crs(); |
| 87 | + mRasterUnitsPerPixelX = std::abs( layer->rasterUnitsPerPixelX() ); |
| 88 | + mRasterUnitsPerPixelY = std::abs( layer->rasterUnitsPerPixelX() ); |
| 89 | + mNbCellsXProvider = mInterface->xSize(); |
| 90 | + mNbCellsYProvider = mInterface->ySize(); |
| 91 | + |
| 92 | + return true; |
| 93 | +} |
| 94 | + |
| 95 | +QVariantMap QgsZonalHistogramAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) |
| 96 | +{ |
| 97 | + |
| 98 | + std::unique_ptr< QgsFeatureSource > zones( parameterAsSource( parameters, QStringLiteral( "INPUT_VECTOR" ), context ) ); |
| 99 | + if ( !zones ) |
| 100 | + throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT_VECTOR" ) ) ); |
| 101 | + |
| 102 | + long count = zones->featureCount(); |
| 103 | + double step = count > 0 ? 100.0 / count : 1; |
| 104 | + long current = 0; |
| 105 | + |
| 106 | + QList< double > uniqueValues; |
| 107 | + QMap< QgsFeatureId, QHash< double, qgssize > > featuresUniqueValues; |
| 108 | + |
| 109 | + // First loop through the zones to build up a list of unique values across all zones to determine sink fields list |
| 110 | + QgsFeatureRequest request; |
| 111 | + request.setSubsetOfAttributes( QgsAttributeList() ); |
| 112 | + if ( zones->sourceCrs() != mCrs ) |
| 113 | + { |
| 114 | + request.setDestinationCrs( mCrs, context.transformContext() ); |
| 115 | + } |
| 116 | + QgsFeatureIterator it = zones->getFeatures( request ); |
| 117 | + QgsFeature f; |
| 118 | + while ( it.nextFeature( f ) ) |
| 119 | + { |
| 120 | + if ( feedback && feedback->isCanceled() ) |
| 121 | + { |
| 122 | + break; |
| 123 | + } |
| 124 | + feedback->setProgress( current * step ); |
| 125 | + |
| 126 | + if ( !f.hasGeometry() ) |
| 127 | + { |
| 128 | + current++; |
| 129 | + continue; |
| 130 | + } |
| 131 | + |
| 132 | + QgsGeometry featureGeometry = f.geometry(); |
| 133 | + QgsRectangle intersectRect = featureGeometry.boundingBox().intersect( &mExtent ); |
| 134 | + if ( intersectRect.isEmpty() ) |
| 135 | + { |
| 136 | + current++; |
| 137 | + continue; |
| 138 | + } |
| 139 | + |
| 140 | + int offsetX, offsetY, nCellsX, nCellsY; |
| 141 | + // Get offset in pixels in x- and y- direction |
| 142 | + offsetX = ( int )( ( intersectRect.xMinimum() - mExtent.xMinimum() ) / mRasterUnitsPerPixelX ); |
| 143 | + offsetY = ( int )( ( mExtent.yMaximum() - intersectRect.yMaximum() ) / mRasterUnitsPerPixelY ); |
| 144 | + |
| 145 | + int maxColumn = ( int )( ( intersectRect.xMaximum() - mExtent.xMinimum() ) / mRasterUnitsPerPixelX ) + 1; |
| 146 | + int maxRow = ( int )( ( mExtent.yMaximum() - intersectRect.yMinimum() ) / mRasterUnitsPerPixelY ) + 1; |
| 147 | + |
| 148 | + nCellsX = maxColumn - offsetX; |
| 149 | + nCellsY = maxRow - offsetY; |
| 150 | + |
| 151 | + // Avoid access to cells outside of the raster (may occur because of rounding) |
| 152 | + if ( ( offsetX + nCellsX ) > mNbCellsXProvider ) |
| 153 | + { |
| 154 | + nCellsX = mNbCellsXProvider - offsetX; |
| 155 | + } |
| 156 | + if ( ( offsetY + nCellsY ) > mNbCellsYProvider ) |
| 157 | + { |
| 158 | + nCellsY = mNbCellsYProvider - offsetY; |
| 159 | + } |
| 160 | + |
| 161 | + QHash< double, qgssize > fUniqueValues; |
| 162 | + middlePoints( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, fUniqueValues ); |
| 163 | + |
| 164 | + if ( fUniqueValues.count() < 1 ) |
| 165 | + { |
| 166 | + // The cell resolution is probably larger than the polygon area. We switch to slower precise pixel - polygon intersection in this case |
| 167 | + preciseIntersection( featureGeometry, offsetX, offsetY, nCellsX, nCellsY, fUniqueValues ); |
| 168 | + } |
| 169 | + |
| 170 | + for ( auto it = fUniqueValues.constBegin(); it != fUniqueValues.constEnd(); ++it ) |
| 171 | + { |
| 172 | + if ( uniqueValues.indexOf( it.key() ) == -1 ) |
| 173 | + { |
| 174 | + uniqueValues << it.key(); |
| 175 | + } |
| 176 | + featuresUniqueValues[f.id()][it.key()] += it.value(); |
| 177 | + } |
| 178 | + |
| 179 | + current++; |
| 180 | + } |
| 181 | + |
| 182 | + std::sort( uniqueValues.begin(), uniqueValues.end() ); |
| 183 | + |
| 184 | + QString fieldPrefix = parameterAsString( parameters, QStringLiteral( "COLUMN_PREFIX" ), context ); |
| 185 | + QgsFields newFields; |
| 186 | + for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it ) |
| 187 | + { |
| 188 | + newFields.append( QgsField( QStringLiteral( "%1%2" ).arg( fieldPrefix ).arg( mHasNoDataValue && *it == mNodataValue ? QStringLiteral( "NODATA" ) : QString::number( *it ) ), QVariant::LongLong, QString(), -1, 0 ) ); |
| 189 | + } |
| 190 | + QgsFields fields = QgsProcessingUtils::combineFields( zones->fields(), newFields ); |
| 191 | + |
| 192 | + QString dest; |
| 193 | + std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, |
| 194 | + zones->wkbType(), zones->sourceCrs() ) ); |
| 195 | + if ( !sink ) |
| 196 | + throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) ); |
| 197 | + |
| 198 | + it = zones->getFeatures( QgsFeatureRequest() ); |
| 199 | + while ( it.nextFeature( f ) ) |
| 200 | + { |
| 201 | + QgsAttributes attributes = f.attributes(); |
| 202 | + QHash< double, qgssize > fUniqueValues = featuresUniqueValues.value( f.id() ); |
| 203 | + for ( auto it = uniqueValues.constBegin(); it != uniqueValues.constEnd(); ++it ) |
| 204 | + { |
| 205 | + attributes += fUniqueValues.value( *it, 0 ); |
| 206 | + } |
| 207 | + |
| 208 | + QgsFeature outputFeature; |
| 209 | + outputFeature.setGeometry( f.geometry() ); |
| 210 | + outputFeature.setAttributes( attributes ); |
| 211 | + |
| 212 | + sink->addFeature( outputFeature, QgsFeatureSink::FastInsert ); |
| 213 | + } |
| 214 | + |
| 215 | + QVariantMap outputs; |
| 216 | + outputs.insert( QStringLiteral( "OUTPUT" ), dest ); |
| 217 | + return outputs; |
| 218 | +} |
| 219 | + |
| 220 | +void QgsZonalHistogramAlgorithm::middlePoints( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues ) |
| 221 | +{ |
| 222 | + double cellCenterX, cellCenterY; |
| 223 | + |
| 224 | + cellCenterY = mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - mRasterUnitsPerPixelY / 2; |
| 225 | + |
| 226 | + geos::unique_ptr polyGeos( QgsGeos::asGeos( poly ) ); |
| 227 | + if ( !polyGeos ) |
| 228 | + { |
| 229 | + return; |
| 230 | + } |
| 231 | + |
| 232 | + GEOSContextHandle_t geosctxt = QgsGeos::getGEOSHandler(); |
| 233 | + geos::prepared_unique_ptr polyGeosPrepared( GEOSPrepare_r( geosctxt, polyGeos.get() ) ); |
| 234 | + if ( !polyGeosPrepared ) |
| 235 | + { |
| 236 | + return; |
| 237 | + } |
| 238 | + |
| 239 | + GEOSCoordSequence *cellCenterCoords = nullptr; |
| 240 | + geos::unique_ptr currentCellCenter; |
| 241 | + |
| 242 | + QgsRectangle blockRect( mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX, |
| 243 | + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - nCellsY * mRasterUnitsPerPixelY, |
| 244 | + mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX + nCellsX * mRasterUnitsPerPixelX, |
| 245 | + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY ); |
| 246 | + std::unique_ptr< QgsRasterBlock > block( mInterface->block( mBand, blockRect, nCellsX, nCellsY ) ); |
| 247 | + for ( int i = 0; i < nCellsY; ++i ) |
| 248 | + { |
| 249 | + cellCenterX = mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX + mRasterUnitsPerPixelX / 2; |
| 250 | + for ( int j = 0; j < nCellsX; ++j ) |
| 251 | + { |
| 252 | + if ( !std::isnan( block->value( i, j ) ) ) |
| 253 | + { |
| 254 | + cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 ); |
| 255 | + GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX ); |
| 256 | + GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY ); |
| 257 | + currentCellCenter.reset( GEOSGeom_createPoint_r( geosctxt, cellCenterCoords ) ); |
| 258 | + if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared.get(), currentCellCenter.get() ) ) |
| 259 | + { |
| 260 | + uniqueValues[block->value( i, j )]++; |
| 261 | + } |
| 262 | + } |
| 263 | + cellCenterX += mRasterUnitsPerPixelX; |
| 264 | + } |
| 265 | + cellCenterY -= mRasterUnitsPerPixelY; |
| 266 | + } |
| 267 | +} |
| 268 | + |
| 269 | +void QgsZonalHistogramAlgorithm::preciseIntersection( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues ) |
| 270 | +{ |
| 271 | + double currentY = mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - mRasterUnitsPerPixelY / 2; |
| 272 | + QgsGeometry pixelRectGeometry; |
| 273 | + |
| 274 | + double hCellSizeX = mRasterUnitsPerPixelX / 2.0; |
| 275 | + double hCellSizeY = mRasterUnitsPerPixelY / 2.0; |
| 276 | + |
| 277 | + QgsRectangle blockRect( mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX, |
| 278 | + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY - nCellsY * mRasterUnitsPerPixelY, |
| 279 | + mExtent.xMinimum() + pixelOffsetX * mRasterUnitsPerPixelX + nCellsX * mRasterUnitsPerPixelX, |
| 280 | + mExtent.yMaximum() - pixelOffsetY * mRasterUnitsPerPixelY ); |
| 281 | + std::unique_ptr< QgsRasterBlock > block( mInterface->block( mBand, blockRect, nCellsX, nCellsY ) ); |
| 282 | + for ( int i = 0; i < nCellsY; ++i ) |
| 283 | + { |
| 284 | + double currentX = mExtent.xMinimum() + mRasterUnitsPerPixelX / 2.0 + pixelOffsetX * mRasterUnitsPerPixelX; |
| 285 | + for ( int j = 0; j < nCellsX; ++j ) |
| 286 | + { |
| 287 | + if ( !std::isnan( block->value( i, j ) ) ) |
| 288 | + { |
| 289 | + pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) ); |
| 290 | + if ( !pixelRectGeometry.isNull() ) |
| 291 | + { |
| 292 | + //intersection |
| 293 | + QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly ); |
| 294 | + if ( !intersectGeometry.isNull() ) |
| 295 | + { |
| 296 | + double intersectionArea = intersectGeometry.area(); |
| 297 | + if ( intersectionArea >= 0.0 ) |
| 298 | + { |
| 299 | + uniqueValues[block->value( i, j )]++; |
| 300 | + } |
| 301 | + } |
| 302 | + pixelRectGeometry = QgsGeometry(); |
| 303 | + } |
| 304 | + } |
| 305 | + currentX += mRasterUnitsPerPixelX; |
| 306 | + } |
| 307 | + currentY -= mRasterUnitsPerPixelY; |
| 308 | + } |
| 309 | +} |
| 310 | + |
| 311 | +///@endcond |
| 312 | + |
| 313 | + |
| 314 | + |
0 commit comments