Skip to content

Commit fdaa57a

Browse files
committedMay 3, 2018
[FEATURE][processing] Zonal histogram algorithm
1 parent 429374e commit fdaa57a

File tree

10 files changed

+549
-0
lines changed

10 files changed

+549
-0
lines changed
 
Binary file not shown.
Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
<?xml version="1.0" encoding="utf-8" ?>
2+
<ogr:FeatureCollection
3+
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
4+
xsi:schemaLocation="http://ogr.maptools.org/ zones_histogram.xsd"
5+
xmlns:ogr="http://ogr.maptools.org/"
6+
xmlns:gml="http://www.opengis.net/gml">
7+
<gml:boundedBy>
8+
<gml:Box>
9+
<gml:coord><gml:X>270743.947032806</gml:X><gml:Y>4458959.32367457</gml:Y></gml:coord>
10+
<gml:coord><gml:X>270848.243623655</gml:X><gml:Y>4459031.472508</gml:Y></gml:coord>
11+
</gml:Box>
12+
</gml:boundedBy>
13+
14+
<gml:featureMember>
15+
<ogr:zones_histogram fid="zones.0">
16+
<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>
17+
<ogr:HISTO_826>4</ogr:HISTO_826>
18+
<ogr:HISTO_837>6</ogr:HISTO_837>
19+
<ogr:HISTO_843>6</ogr:HISTO_843>
20+
<ogr:HISTO_851>9</ogr:HISTO_851>
21+
<ogr:HISTO_880>0</ogr:HISTO_880>
22+
</ogr:zones_histogram>
23+
</gml:featureMember>
24+
<gml:featureMember>
25+
<ogr:zones_histogram fid="zones.1">
26+
<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>
27+
<ogr:HISTO_826>0</ogr:HISTO_826>
28+
<ogr:HISTO_837>0</ogr:HISTO_837>
29+
<ogr:HISTO_843>0</ogr:HISTO_843>
30+
<ogr:HISTO_851>0</ogr:HISTO_851>
31+
<ogr:HISTO_880>6</ogr:HISTO_880>
32+
</ogr:zones_histogram>
33+
</gml:featureMember>
34+
</ogr:FeatureCollection>
Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
<?xml version="1.0" encoding="UTF-8"?>
2+
<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">
3+
<xs:import namespace="http://www.opengis.net/gml" schemaLocation="http://schemas.opengis.net/gml/2.1.2/feature.xsd"/>
4+
<xs:element name="FeatureCollection" type="ogr:FeatureCollectionType" substitutionGroup="gml:_FeatureCollection"/>
5+
<xs:complexType name="FeatureCollectionType">
6+
<xs:complexContent>
7+
<xs:extension base="gml:AbstractFeatureCollectionType">
8+
<xs:attribute name="lockId" type="xs:string" use="optional"/>
9+
<xs:attribute name="scope" type="xs:string" use="optional"/>
10+
</xs:extension>
11+
</xs:complexContent>
12+
</xs:complexType>
13+
<xs:element name="zones_histogram" type="ogr:zones_histogram_Type" substitutionGroup="gml:_Feature"/>
14+
<xs:complexType name="zones_histogram_Type">
15+
<xs:complexContent>
16+
<xs:extension base="gml:AbstractFeatureType">
17+
<xs:sequence>
18+
<xs:element name="geometryProperty" type="gml:MultiPolygonPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
19+
<xs:element name="HISTO_826" nillable="true" minOccurs="0" maxOccurs="1">
20+
<xs:simpleType>
21+
<xs:restriction base="xs:long">
22+
<xs:totalDigits value="20"/>
23+
</xs:restriction>
24+
</xs:simpleType>
25+
</xs:element>
26+
<xs:element name="HISTO_837" nillable="true" minOccurs="0" maxOccurs="1">
27+
<xs:simpleType>
28+
<xs:restriction base="xs:long">
29+
<xs:totalDigits value="20"/>
30+
</xs:restriction>
31+
</xs:simpleType>
32+
</xs:element>
33+
<xs:element name="HISTO_843" nillable="true" minOccurs="0" maxOccurs="1">
34+
<xs:simpleType>
35+
<xs:restriction base="xs:long">
36+
<xs:totalDigits value="20"/>
37+
</xs:restriction>
38+
</xs:simpleType>
39+
</xs:element>
40+
<xs:element name="HISTO_851" nillable="true" minOccurs="0" maxOccurs="1">
41+
<xs:simpleType>
42+
<xs:restriction base="xs:long">
43+
<xs:totalDigits value="20"/>
44+
</xs:restriction>
45+
</xs:simpleType>
46+
</xs:element>
47+
<xs:element name="HISTO_880" nillable="true" minOccurs="0" maxOccurs="1">
48+
<xs:simpleType>
49+
<xs:restriction base="xs:long">
50+
<xs:totalDigits value="20"/>
51+
</xs:restriction>
52+
</xs:simpleType>
53+
</xs:element>
54+
</xs:sequence>
55+
</xs:extension>
56+
</xs:complexContent>
57+
</xs:complexType>
58+
</xs:schema>

‎python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3460,6 +3460,22 @@ tests:
34603460
geometry:
34613461
precision: 5
34623462

3463+
- algorithm: qgis:zonalhistogram
3464+
name: zonal histogram
3465+
params:
3466+
COLUMN_PREFIX: HISTO_
3467+
INPUT_RASTER:
3468+
name: raster.tif
3469+
type: raster
3470+
INPUT_VECTOR:
3471+
name: zones.gml
3472+
type: vector
3473+
RASTER_BAND: 1
3474+
results:
3475+
OUTPUT:
3476+
name: expected/zones_histogram.gml
3477+
type: vector
3478+
34633479
- algorithm: native:fixgeometries
34643480
name: Fix geometries
34653481
params:
Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
<?xml version="1.0" encoding="utf-8" ?>
2+
<ogr:FeatureCollection
3+
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
4+
xsi:schemaLocation="http://ogr.maptools.org/ zones.xsd"
5+
xmlns:ogr="http://ogr.maptools.org/"
6+
xmlns:gml="http://www.opengis.net/gml">
7+
<gml:boundedBy>
8+
<gml:Box>
9+
<gml:coord><gml:X>270743.9470328057</gml:X><gml:Y>4458959.323674565</gml:Y></gml:coord>
10+
<gml:coord><gml:X>270848.2436236552</gml:X><gml:Y>4459031.472508</gml:Y></gml:coord>
11+
</gml:Box>
12+
</gml:boundedBy>
13+
14+
<gml:featureMember>
15+
<ogr:zones fid="zones.0">
16+
<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>
17+
</ogr:zones>
18+
</gml:featureMember>
19+
<gml:featureMember>
20+
<ogr:zones fid="zones.1">
21+
<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>
22+
</ogr:zones>
23+
</gml:featureMember>
24+
</ogr:FeatureCollection>
Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
<?xml version="1.0" encoding="UTF-8"?>
2+
<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">
3+
<xs:import namespace="http://www.opengis.net/gml" schemaLocation="http://schemas.opengis.net/gml/2.1.2/feature.xsd"/>
4+
<xs:element name="FeatureCollection" type="ogr:FeatureCollectionType" substitutionGroup="gml:_FeatureCollection"/>
5+
<xs:complexType name="FeatureCollectionType">
6+
<xs:complexContent>
7+
<xs:extension base="gml:AbstractFeatureCollectionType">
8+
<xs:attribute name="lockId" type="xs:string" use="optional"/>
9+
<xs:attribute name="scope" type="xs:string" use="optional"/>
10+
</xs:extension>
11+
</xs:complexContent>
12+
</xs:complexType>
13+
<xs:element name="zones" type="ogr:zones_Type" substitutionGroup="gml:_Feature"/>
14+
<xs:complexType name="zones_Type">
15+
<xs:complexContent>
16+
<xs:extension base="gml:AbstractFeatureType">
17+
<xs:sequence>
18+
<xs:element name="geometryProperty" type="gml:MultiPolygonPropertyType" nillable="true" minOccurs="0" maxOccurs="1"/>
19+
</xs:sequence>
20+
</xs:extension>
21+
</xs:complexContent>
22+
</xs:complexType>
23+
</xs:schema>

‎src/analysis/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,7 @@ SET(QGIS_ANALYSIS_SRCS
8282
processing/qgsalgorithmunion.cpp
8383
processing/qgsalgorithmuniquevalueindex.cpp
8484
processing/qgsalgorithmwedgebuffers.cpp
85+
processing/qgsalgorithmzonalhistogram.cpp
8586

8687
processing/qgsnativealgorithms.cpp
8788
processing/qgsoverlayutils.cpp
Lines changed: 314 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,314 @@
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 &parameters, 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 &parameters, 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+
Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
/***************************************************************************
2+
qgsalgorithmzonalhistogram.h
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+
#ifndef QGSALGORITHMZONALHISTOGRAM_H
19+
#define QGSALGORITHMZONALHISTOGRAM_H
20+
21+
#define SIP_NO_FILE
22+
23+
#include "qgis.h"
24+
#include "qgsprocessingalgorithm.h"
25+
26+
///@cond PRIVATE
27+
28+
/**
29+
* Native zonal histogram algorithm.
30+
*/
31+
class QgsZonalHistogramAlgorithm : public QgsProcessingAlgorithm
32+
{
33+
34+
public:
35+
36+
QgsZonalHistogramAlgorithm() = default;
37+
void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
38+
QString name() const override;
39+
QString displayName() const override;
40+
QStringList tags() const override;
41+
QString group() const override;
42+
QString groupId() const override;
43+
QString shortHelpString() const override;
44+
QgsZonalHistogramAlgorithm *createInstance() const override SIP_FACTORY;
45+
46+
protected:
47+
48+
bool prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
49+
QVariantMap processAlgorithm( const QVariantMap &parameters,
50+
QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
51+
52+
//! Fetch unique values by considering the pixels where the center point is within the polygon (fast)
53+
void middlePoints( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues );
54+
55+
//! Fetch unique values with precise pixel - polygon intersection test (slow)
56+
void preciseIntersection( const QgsGeometry &poly, int pixelOffsetX, int pixelOffsetY, int nCellsX, int nCellsY, QHash< double, qgssize > &uniqueValues );
57+
58+
private:
59+
60+
std::unique_ptr< QgsRasterInterface > mInterface;
61+
int mBand;
62+
bool mHasNoDataValue = false;
63+
float mNodataValue = -1;
64+
QgsRectangle mExtent;
65+
QgsCoordinateReferenceSystem mCrs;
66+
double mRasterUnitsPerPixelX;
67+
double mRasterUnitsPerPixelY;
68+
double mNbCellsXProvider;
69+
double mNbCellsYProvider;
70+
71+
};
72+
73+
///@endcond PRIVATE
74+
75+
#endif // QGSALGORITHMZONALHISTOGRAM_H
76+
77+

‎src/analysis/processing/qgsnativealgorithms.cpp

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,6 +79,7 @@
7979
#include "qgsalgorithmunion.h"
8080
#include "qgsalgorithmuniquevalueindex.h"
8181
#include "qgsalgorithmwedgebuffers.h"
82+
#include "qgsalgorithmzonalhistogram.h"
8283

8384

8485
///@cond PRIVATE
@@ -186,6 +187,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
186187
addAlgorithm( new QgsUnionAlgorithm() );
187188
addAlgorithm( new QgsVariableWidthBufferByMAlgorithm() );
188189
addAlgorithm( new QgsWedgeBuffersAlgorithm() );
190+
addAlgorithm( new QgsZonalHistogramAlgorithm() );
189191
}
190192

191193

0 commit comments

Comments
 (0)
Please sign in to comment.