sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-04-14 00:07:35 -04:00
Code Issues Packages Projects Releases Wiki Activity

[FEATURE][processing] Zonal histogram algorithm

Browse Source
This commit is contained in:
nirvn 2018-05-01 11:42:31 +07:00 committed by Mathieu Pellerin
parent 429374ecd1
commit fdaa57a273
10 changed files with 549 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
python/plugins/processing/tests/testdata/custom/circular_strings.gpkg vendored
View File

Binary file not shown.

34
python/plugins/processing/tests/testdata/expected/zones_histogram.gml vendored Normal file
View File

@ -0,0 +1,34 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://ogr.maptools.org/ zones_histogram.xsd"
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>270743.947032806</gml:X><gml:Y>4458959.32367457</gml:Y></gml:coord>
<gml:coord><gml:X>270848.243623655</gml:X><gml:Y>4459031.472508</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>
<gml:featureMember>
<ogr:zones_histogram fid="zones.0">
<ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:23030"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>270744.399022365,4459030.56852888 270797.959785163,4459031.472508 270797.959785163,4458978.58972954 270743.947032806,4458979.0417191 270744.399022365,4459030.56852888</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
<ogr:HISTO_826>4</ogr:HISTO_826>
<ogr:HISTO_837>6</ogr:HISTO_837>
<ogr:HISTO_843>6</ogr:HISTO_843>
<ogr:HISTO_851>9</ogr:HISTO_851>
<ogr:HISTO_880>0</ogr:HISTO_880>
</ogr:zones_histogram>
</gml:featureMember>
<gml:featureMember>
<ogr:zones_histogram fid="zones.1">
<ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:23030"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>270818.41231273,4458978.87222301 270818.18631795,4458959.54966935 270848.243623655,4458959.32367457 270848.017628875,4458978.87222301 270818.41231273,4458978.87222301</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
<ogr:HISTO_826>0</ogr:HISTO_826>
<ogr:HISTO_837>0</ogr:HISTO_837>
<ogr:HISTO_843>0</ogr:HISTO_843>
<ogr:HISTO_851>0</ogr:HISTO_851>
<ogr:HISTO_880>6</ogr:HISTO_880>
</ogr:zones_histogram>
</gml:featureMember>
</ogr:FeatureCollection>

58
python/plugins/processing/tests/testdata/expected/zones_histogram.xsd vendored Normal file
View File

@ -0,0 +1,58 @@
<?xml version="1.0" encoding="UTF-8"?>
<xs:schema targetNamespace="http://ogr.maptools.org/" xmlns:ogr="http://ogr.maptools.org/" xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:gml="http://www.opengis.net/gml" elementFormDefault="qualified" version="1.0">
<xs:import namespace="http://www.opengis.net/gml" schemaLocation="http://schemas.opengis.net/gml/2.1.2/feature.xsd"/>
<xs:element name="FeatureCollection" type="ogr:FeatureCollectionType" substitutionGroup="gml:_FeatureCollection"/>
<xs:complexType name="FeatureCollectionType">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureCollectionType">
<xs:attribute name="lockId" type="xs:string" use="optional"/>
<xs:attribute name="scope" type="xs:string" use="optional"/>
</xs:extension>
</xs:complexContent>
</xs:complexType>
<xs:element name="zones_histogram" type="ogr:zones_histogram_Type" substitutionGroup="gml:_Feature"/>
<xs:complexType name="zones_histogram_Type">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureType">
<xs:sequence>
<xs:element name="geometryProperty" type="gml:MultiPolygonPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
<xs:element name="HISTO_826" nillable="true" minOccurs="0" maxOccurs="1">
<xs:simpleType>
<xs:restriction base="xs:long">
<xs:totalDigits value="20"/>
</xs:restriction>
</xs:simpleType>
</xs:element>
<xs:element name="HISTO_837" nillable="true" minOccurs="0" maxOccurs="1">
<xs:simpleType>
<xs:restriction base="xs:long">
<xs:totalDigits value="20"/>
</xs:restriction>
</xs:simpleType>
</xs:element>
<xs:element name="HISTO_843" nillable="true" minOccurs="0" maxOccurs="1">
<xs:simpleType>
<xs:restriction base="xs:long">
<xs:totalDigits value="20"/>
</xs:restriction>
</xs:simpleType>
</xs:element>
<xs:element name="HISTO_851" nillable="true" minOccurs="0" maxOccurs="1">
<xs:simpleType>
<xs:restriction base="xs:long">
<xs:totalDigits value="20"/>
</xs:restriction>
</xs:simpleType>
</xs:element>
<xs:element name="HISTO_880" nillable="true" minOccurs="0" maxOccurs="1">
<xs:simpleType>
<xs:restriction base="xs:long">
<xs:totalDigits value="20"/>
</xs:restriction>
</xs:simpleType>
</xs:element>
</xs:sequence>
</xs:extension>
</xs:complexContent>
</xs:complexType>
</xs:schema>

16
python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml vendored
View File

@ -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:

24
python/plugins/processing/tests/testdata/zones.gml vendored Normal file
View File

@ -0,0 +1,24 @@
<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://ogr.maptools.org/ zones.xsd"
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:boundedBy>
<gml:Box>
<gml:coord><gml:X>270743.9470328057</gml:X><gml:Y>4458959.323674565</gml:Y></gml:coord>
<gml:coord><gml:X>270848.2436236552</gml:X><gml:Y>4459031.472508</gml:Y></gml:coord>
</gml:Box>
</gml:boundedBy>
<gml:featureMember>
<ogr:zones fid="zones.0">
<ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:23030"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>270744.399022365,4459030.56852888 270797.959785163,4459031.472508 270797.959785163,4458978.58972954 270743.947032806,4458979.0417191 270744.399022365,4459030.56852888</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
</ogr:zones>
</gml:featureMember>
<gml:featureMember>
<ogr:zones fid="zones.1">
<ogr:geometryProperty><gml:MultiPolygon srsName="EPSG:23030"><gml:polygonMember><gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>270818.41231273,4458978.87222301 270818.18631795,4458959.54966935 270848.243623655,4458959.32367457 270848.017628875,4458978.87222301 270818.41231273,4458978.87222301</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon></gml:polygonMember></gml:MultiPolygon></ogr:geometryProperty>
</ogr:zones>
</gml:featureMember>
</ogr:FeatureCollection>

23
python/plugins/processing/tests/testdata/zones.xsd vendored Normal file
View File

@ -0,0 +1,23 @@
<?xml version="1.0" encoding="UTF-8"?>
<xs:schema targetNamespace="http://ogr.maptools.org/" xmlns:ogr="http://ogr.maptools.org/" xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:gml="http://www.opengis.net/gml" elementFormDefault="qualified" version="1.0">
<xs:import namespace="http://www.opengis.net/gml" schemaLocation="http://schemas.opengis.net/gml/2.1.2/feature.xsd"/>
<xs:element name="FeatureCollection" type="ogr:FeatureCollectionType" substitutionGroup="gml:_FeatureCollection"/>
<xs:complexType name="FeatureCollectionType">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureCollectionType">
<xs:attribute name="lockId" type="xs:string" use="optional"/>
<xs:attribute name="scope" type="xs:string" use="optional"/>
</xs:extension>
</xs:complexContent>
</xs:complexType>
<xs:element name="zones" type="ogr:zones_Type" substitutionGroup="gml:_Feature"/>
<xs:complexType name="zones_Type">
<xs:complexContent>
<xs:extension base="gml:AbstractFeatureType">
<xs:sequence>
<xs:element name="geometryProperty" type="gml:MultiPolygonPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
</xs:sequence>
</xs:extension>
</xs:complexContent>
</xs:complexType>
</xs:schema>

1
src/analysis/CMakeLists.txt
View File

@ -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

314
src/analysis/processing/qgsalgorithmzonalhistogram.cpp Normal file
View File

@ -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 &parameters, 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 &parameters, 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

77
src/analysis/processing/qgsalgorithmzonalhistogram.h Normal file
View File

@ -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 &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
QVariantMap processAlgorithm( const QVariantMap &parameters,
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

2
src/analysis/processing/qgsnativealgorithms.cpp
View File

@ -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() );
}