QGIS API Documentation 3.41.0-Master (cea29feecf2)
Loading...
Searching...
No Matches
qgskde.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgskde.cpp
3 ----------
4 Date : October 2016
5 Copyright : (C) 2016 by Nyall Dawson
6 Email : nyall dot dawson at gmail dot com
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15
16#include "qgskde.h"
17#include "qgsfeaturesource.h"
18#include "qgsfeatureiterator.h"
19#include "qgsgeometry.h"
20
21#define NO_DATA -9999
22
23QgsKernelDensityEstimation::QgsKernelDensityEstimation( const QgsKernelDensityEstimation::Parameters &parameters, const QString &outputFile, const QString &outputFormat )
24 : mSource( parameters.source )
25 , mOutputFile( outputFile )
26 , mOutputFormat( outputFormat )
27 , mRadiusField( -1 )
28 , mWeightField( -1 )
29 , mRadius( parameters.radius )
30 , mPixelSize( parameters.pixelSize )
31 , mShape( parameters.shape )
32 , mDecay( parameters.decayRatio )
33 , mOutputValues( parameters.outputValues )
34 , mBufferSize( -1 )
35 , mDatasetH( nullptr )
36 , mRasterBandH( nullptr )
37{
38 if ( !parameters.radiusField.isEmpty() )
39 mRadiusField = mSource->fields().lookupField( parameters.radiusField );
40 if ( !parameters.weightField.isEmpty() )
41 mWeightField = mSource->fields().lookupField( parameters.weightField );
42}
43
45{
46 const Result result = prepare();
47 if ( result != Success )
48 return result;
49
50 QgsAttributeList requiredAttributes;
51
52 if ( mRadiusField >= 0 )
53 requiredAttributes << mRadiusField;
54
55 if ( mWeightField >= 0 )
56 requiredAttributes << mWeightField;
57
58 QgsFeatureIterator fit = mSource->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes ) );
59
60 QgsFeature f;
61 while ( fit.nextFeature( f ) )
62 {
63 addFeature( f );
64 }
65
66 return finalise();
67}
68
70{
71 GDALAllRegister();
72
73 GDALDriverH driver = GDALGetDriverByName( mOutputFormat.toUtf8() );
74 if ( !driver )
75 {
76 return DriverError;
77 }
78
79 if ( !mSource )
80 return InvalidParameters;
81
82 mBounds = calculateBounds();
83 if ( mBounds.isNull() )
84 return InvalidParameters;
85
86 const int rows = std::max( std::ceil( mBounds.height() / mPixelSize ) + 1, 1.0 );
87 const int cols = std::max( std::ceil( mBounds.width() / mPixelSize ) + 1, 1.0 );
88
89 if ( !createEmptyLayer( driver, mBounds, rows, cols ) )
90 return FileCreationError;
91
92 // open the raster in GA_Update mode
93 mDatasetH.reset( GDALOpen( mOutputFile.toUtf8().constData(), GA_Update ) );
94 if ( !mDatasetH )
95 return FileCreationError;
96 mRasterBandH = GDALGetRasterBand( mDatasetH.get(), 1 );
97 if ( !mRasterBandH )
98 return FileCreationError;
99
100 mBufferSize = -1;
101 if ( mRadiusField < 0 )
102 mBufferSize = radiusSizeInPixels( mRadius );
103
104 return Success;
105}
106
108{
109 const QgsGeometry featureGeometry = feature.geometry();
110 if ( featureGeometry.isNull() )
111 {
112 return Success;
113 }
114
115 // convert the geometry to multipoint
116 QgsMultiPointXY multiPoints;
117 if ( !featureGeometry.isMultipart() )
118 {
119 const QgsPointXY p = featureGeometry.asPoint();
120 // avoiding any empty points or out of extent points
121 if ( !mBounds.contains( p ) )
122 return Success;
123
124 multiPoints << p;
125 }
126 else
127 {
128 multiPoints = featureGeometry.asMultiPoint();
129 }
130
131 // if radius is variable then fetch it and calculate new pixel buffer size
132 double radius = mRadius;
133 int buffer = mBufferSize;
134 if ( mRadiusField >= 0 )
135 {
136 radius = feature.attribute( mRadiusField ).toDouble();
137 buffer = radiusSizeInPixels( radius );
138 }
139 const int blockSize = 2 * buffer + 1; //Block SIDE would be more appropriate
140
141 // calculate weight
142 double weight = 1.0;
143 if ( mWeightField >= 0 )
144 {
145 weight = feature.attribute( mWeightField ).toDouble();
146 }
147
148 Result result = Success;
149
150 //loop through all points in multipoint
151 for ( QgsMultiPointXY::const_iterator pointIt = multiPoints.constBegin(); pointIt != multiPoints.constEnd(); ++pointIt )
152 {
153 // avoiding any empty points or out of extent points
154 if ( !mBounds.contains( *pointIt ) )
155 {
156 continue;
157 }
158
159 // calculate the pixel position
160 const unsigned int xPosition = ( ( ( *pointIt ).x() - mBounds.xMinimum() ) / mPixelSize ) - buffer;
161 const unsigned int yPosition = ( ( ( *pointIt ).y() - mBounds.yMinimum() ) / mPixelSize ) - buffer;
162 const unsigned int yPositionIO = ( ( mBounds.yMaximum() - ( *pointIt ).y() ) / mPixelSize ) - buffer;
163
164
165 // get the data
166 float *dataBuffer = ( float * ) CPLMalloc( sizeof( float ) * blockSize * blockSize );
167 if ( GDALRasterIO( mRasterBandH, GF_Read, xPosition, yPositionIO, blockSize, blockSize, dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
168 {
169 result = RasterIoError;
170 }
171
172 for ( int xp = 0; xp < blockSize; xp++ )
173 {
174 for ( int yp = 0; yp < blockSize; yp++ )
175 {
176 const double pixelCentroidX = ( xPosition + xp + 0.5 ) * mPixelSize + mBounds.xMinimum();
177 const double pixelCentroidY = ( yPosition + yp + 0.5 ) * mPixelSize + mBounds.yMinimum();
178
179 const double distance = std::sqrt( std::pow( pixelCentroidX - ( *pointIt ).x(), 2.0 ) + std::pow( pixelCentroidY - ( *pointIt ).y(), 2.0 ) );
180
181 // is pixel outside search bandwidth of feature?
182 if ( distance > radius )
183 {
184 continue;
185 }
186
187 const double pixelValue = weight * calculateKernelValue( distance, radius, mShape, mOutputValues );
188 const int pos = xp + blockSize * yp;
189 if ( dataBuffer[pos] == NO_DATA )
190 {
191 dataBuffer[pos] = 0;
192 }
193 dataBuffer[pos] += pixelValue;
194 }
195 }
196 if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPositionIO, blockSize, blockSize, dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
197 {
198 result = RasterIoError;
199 }
200 CPLFree( dataBuffer );
201 }
202
203 return result;
204}
205
207{
208 mDatasetH.reset();
209 mRasterBandH = nullptr;
210 return Success;
211}
212
213int QgsKernelDensityEstimation::radiusSizeInPixels( double radius ) const
214{
215 int buffer = radius / mPixelSize;
216 if ( radius - ( mPixelSize * buffer ) > 0.5 )
217 {
218 ++buffer;
219 }
220 return buffer;
221}
222
223bool QgsKernelDensityEstimation::createEmptyLayer( GDALDriverH driver, const QgsRectangle &bounds, int rows, int columns ) const
224{
225 double geoTransform[6] = { bounds.xMinimum(), mPixelSize, 0, bounds.yMaximum(), 0, -mPixelSize };
226 const gdal::dataset_unique_ptr emptyDataset( GDALCreate( driver, mOutputFile.toUtf8(), columns, rows, 1, GDT_Float32, nullptr ) );
227 if ( !emptyDataset )
228 return false;
229
230 if ( GDALSetGeoTransform( emptyDataset.get(), geoTransform ) != CE_None )
231 return false;
232
233 // Set the projection on the raster destination to match the input layer
234 if ( GDALSetProjection( emptyDataset.get(), mSource->sourceCrs().toWkt().toLocal8Bit().data() ) != CE_None )
235 return false;
236
237 GDALRasterBandH poBand = GDALGetRasterBand( emptyDataset.get(), 1 );
238 if ( !poBand )
239 return false;
240
241 if ( GDALSetRasterNoDataValue( poBand, NO_DATA ) != CE_None )
242 return false;
243
244 float *line = static_cast<float *>( CPLMalloc( sizeof( float ) * columns ) );
245 for ( int i = 0; i < columns; i++ )
246 {
247 line[i] = NO_DATA;
248 }
249 // Write the empty raster
250 for ( int i = 0; i < rows; i++ )
251 {
252 if ( GDALRasterIO( poBand, GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 ) != CE_None )
253 {
254 return false;
255 }
256 }
257
258 CPLFree( line );
259 return true;
260}
261
262double QgsKernelDensityEstimation::calculateKernelValue( const double distance, const double bandwidth, const QgsKernelDensityEstimation::KernelShape shape, const QgsKernelDensityEstimation::OutputValues outputType ) const
263{
264 switch ( shape )
265 {
266 case KernelTriangular:
267 return triangularKernel( distance, bandwidth, outputType );
268
269 case KernelUniform:
270 return uniformKernel( distance, bandwidth, outputType );
271
272 case KernelQuartic:
273 return quarticKernel( distance, bandwidth, outputType );
274
275 case KernelTriweight:
276 return triweightKernel( distance, bandwidth, outputType );
277
279 return epanechnikovKernel( distance, bandwidth, outputType );
280 }
281 return 0; //no warnings
282}
283
284/* The kernel functions below are taken from "Kernel Smoothing" by Wand and Jones (1995), p. 175
285 *
286 * Each kernel is multiplied by a normalizing constant "k", which normalizes the kernel area
287 * to 1 for a given bandwidth size.
288 *
289 * k is calculated by polar double integration of the kernel function
290 * between a radius of 0 to the specified bandwidth and equating the area to 1. */
291
292double QgsKernelDensityEstimation::uniformKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
293{
294 Q_UNUSED( distance )
295 switch ( outputType )
296 {
297 case OutputScaled:
298 {
299 // Normalizing constant
300 const double k = 2. / ( M_PI * bandwidth );
301
302 // Derived from Wand and Jones (1995), p. 175
303 return k * ( 0.5 / bandwidth );
304 }
305 case OutputRaw:
306 return 1.0;
307 }
308 return 0.0; // NO warnings!!!!!
309}
310
311double QgsKernelDensityEstimation::quarticKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
312{
313 switch ( outputType )
314 {
315 case OutputScaled:
316 {
317 // Normalizing constant
318 const double k = 116. / ( 5. * M_PI * std::pow( bandwidth, 2 ) );
319
320 // Derived from Wand and Jones (1995), p. 175
321 return k * ( 15. / 16. ) * std::pow( 1. - std::pow( distance / bandwidth, 2 ), 2 );
322 }
323 case OutputRaw:
324 return std::pow( 1. - std::pow( distance / bandwidth, 2 ), 2 );
325 }
326 return 0.0; //no, seriously, I told you NO WARNINGS!
327}
328
329double QgsKernelDensityEstimation::triweightKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
330{
331 switch ( outputType )
332 {
333 case OutputScaled:
334 {
335 // Normalizing constant
336 const double k = 128. / ( 35. * M_PI * std::pow( bandwidth, 2 ) );
337
338 // Derived from Wand and Jones (1995), p. 175
339 return k * ( 35. / 32. ) * std::pow( 1. - std::pow( distance / bandwidth, 2 ), 3 );
340 }
341 case OutputRaw:
342 return std::pow( 1. - std::pow( distance / bandwidth, 2 ), 3 );
343 }
344 return 0.0; // this is getting ridiculous... don't you ever listen to a word I say?
345}
346
347double QgsKernelDensityEstimation::epanechnikovKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
348{
349 switch ( outputType )
350 {
351 case OutputScaled:
352 {
353 // Normalizing constant
354 const double k = 8. / ( 3. * M_PI * std::pow( bandwidth, 2 ) );
355
356 // Derived from Wand and Jones (1995), p. 175
357 return k * ( 3. / 4. ) * ( 1. - std::pow( distance / bandwidth, 2 ) );
358 }
359 case OutputRaw:
360 return ( 1. - std::pow( distance / bandwidth, 2 ) );
361 }
362
363 return 0.0; // la la la i'm not listening
364}
365
366double QgsKernelDensityEstimation::triangularKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
367{
368 switch ( outputType )
369 {
370 case OutputScaled:
371 {
372 // Normalizing constant. In this case it's calculated a little different
373 // due to the inclusion of the non-standard "decay" parameter
374
375 if ( mDecay >= 0 )
376 {
377 const double k = 3. / ( ( 1. + 2. * mDecay ) * M_PI * std::pow( bandwidth, 2 ) );
378
379 // Derived from Wand and Jones (1995), p. 175 (with addition of decay parameter)
380 return k * ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
381 }
382 else
383 {
384 // Non-standard or mathematically valid negative decay ("coolmap")
385 return ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
386 }
387 }
388 case OutputRaw:
389 return ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
390 }
391 return 0.0; // ....
392}
393
394QgsRectangle QgsKernelDensityEstimation::calculateBounds() const
395{
396 if ( !mSource )
397 return QgsRectangle();
398
399 QgsRectangle bbox = mSource->sourceExtent();
400
401 double radius = 0;
402 if ( mRadiusField >= 0 )
403 {
404 // if radius is using a field, find the max value
405 radius = mSource->maximumValue( mRadiusField ).toDouble();
406 }
407 else
408 {
409 radius = mRadius;
410 }
411 // expand the bounds by the maximum search radius
412 bbox.setXMinimum( bbox.xMinimum() - radius );
413 bbox.setYMinimum( bbox.yMinimum() - radius );
414 bbox.setXMaximum( bbox.xMaximum() + radius );
415 bbox.setYMaximum( bbox.yMaximum() + radius );
416 return bbox;
417}
QString toWkt(Qgis::CrsWktVariant variant=Qgis::CrsWktVariant::Wkt1Gdal, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
virtual QgsFields fields() const =0
Returns the fields associated with features in the source.
virtual QgsCoordinateReferenceSystem sourceCrs() const =0
Returns the coordinate reference system for features in the source.
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
virtual QVariant maximumValue(int fieldIndex) const
Returns the maximum value for an attribute column or an invalid variant in case of error.
virtual QgsRectangle sourceExtent() const
Returns the extent of all geometries from the source.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsGeometry geometry
Definition qgsfeature.h:69
Q_INVOKABLE QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Q_INVOKABLE int lookupField(const QString &fieldName) const
Looks up field's index from the field name.
A geometry is the spatial representation of a feature.
QgsMultiPointXY asMultiPoint() const
Returns the contents of the geometry as a multi-point.
QgsPointXY asPoint() const
Returns the contents of the geometry as a 2-dimensional point.
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
QgsKernelDensityEstimation(const Parameters &parameters, const QString &outputFile, const QString &outputFormat)
Constructor for QgsKernelDensityEstimation.
Definition qgskde.cpp:23
Result run()
Runs the KDE calculation across the whole layer at once.
Definition qgskde.cpp:44
Result addFeature(const QgsFeature &feature)
Adds a single feature to the KDE surface.
Definition qgskde.cpp:107
Result finalise()
Finalises the output file.
Definition qgskde.cpp:206
Result
Result of operation.
Definition qgskde.h:60
@ DriverError
Could not open the driver for the specified format.
Definition qgskde.h:62
@ FileCreationError
Error creating output file.
Definition qgskde.h:64
@ RasterIoError
Error writing to raster.
Definition qgskde.h:65
@ Success
Operation completed successfully.
Definition qgskde.h:61
@ InvalidParameters
Input parameters were not valid.
Definition qgskde.h:63
OutputValues
Output values type.
Definition qgskde.h:53
@ OutputRaw
Output the raw KDE values.
Definition qgskde.h:54
@ OutputScaled
Output mathematically correct scaled values.
Definition qgskde.h:55
KernelShape
Kernel shape type.
Definition qgskde.h:43
@ KernelTriweight
Triweight kernel.
Definition qgskde.h:47
@ KernelUniform
Uniform (flat) kernel.
Definition qgskde.h:46
@ KernelEpanechnikov
Epanechnikov kernel.
Definition qgskde.h:48
@ KernelTriangular
Triangular kernel.
Definition qgskde.h:45
@ KernelQuartic
Quartic kernel.
Definition qgskde.h:44
Result prepare()
Prepares the output file for writing and setups up the surface calculation.
Definition qgskde.cpp:69
A class to represent a 2D point.
Definition qgspointxy.h:60
A rectangle specified with double values.
bool contains(const QgsRectangle &rect) const
Returns true when rectangle contains other rectangle.
double xMinimum
double yMinimum
double xMaximum
void setYMinimum(double y)
Set the minimum y value.
void setXMinimum(double x)
Set the minimum x value.
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
double yMaximum
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
QList< int > QgsAttributeList
Definition qgsfield.h:27
QVector< QgsPointXY > QgsMultiPointXY
A collection of QgsPoints that share a common collection of attributes.
Definition qgsgeometry.h:80
#define NO_DATA
Definition qgskde.cpp:21
QString radiusField
Field for radius, or empty if using a fixed radius.
Definition qgskde.h:78
QString weightField
Field name for weighting field, or empty if not using weights.
Definition qgskde.h:81