QGIS API Documentation 3.41.0-Master (181b2f43d8e)
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,
168 dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
169 {
170 result = RasterIoError;
171 }
172
173 for ( int xp = 0; xp < blockSize; xp++ )
174 {
175 for ( int yp = 0; yp < blockSize; yp++ )
176 {
177 const double pixelCentroidX = ( xPosition + xp + 0.5 ) * mPixelSize + mBounds.xMinimum();
178 const double pixelCentroidY = ( yPosition + yp + 0.5 ) * mPixelSize + mBounds.yMinimum();
179
180 const double distance = std::sqrt( std::pow( pixelCentroidX - ( *pointIt ).x(), 2.0 ) + std::pow( pixelCentroidY - ( *pointIt ).y(), 2.0 ) );
181
182 // is pixel outside search bandwidth of feature?
183 if ( distance > radius )
184 {
185 continue;
186 }
187
188 const double pixelValue = weight * calculateKernelValue( distance, radius, mShape, mOutputValues );
189 const int pos = xp + blockSize * yp;
190 if ( dataBuffer[ pos ] == NO_DATA )
191 {
192 dataBuffer[ pos ] = 0;
193 }
194 dataBuffer[ pos ] += pixelValue;
195 }
196 }
197 if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPositionIO, blockSize, blockSize,
198 dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
199 {
200 result = RasterIoError;
201 }
202 CPLFree( dataBuffer );
203 }
204
205 return result;
206}
207
209{
210 mDatasetH.reset();
211 mRasterBandH = nullptr;
212 return Success;
213}
214
215int QgsKernelDensityEstimation::radiusSizeInPixels( double radius ) const
216{
217 int buffer = radius / mPixelSize;
218 if ( radius - ( mPixelSize * buffer ) > 0.5 )
219 {
220 ++buffer;
221 }
222 return buffer;
223}
224
225bool QgsKernelDensityEstimation::createEmptyLayer( GDALDriverH driver, const QgsRectangle &bounds, int rows, int columns ) const
226{
227 double geoTransform[6] = { bounds.xMinimum(), mPixelSize, 0, bounds.yMaximum(), 0, -mPixelSize };
228 const gdal::dataset_unique_ptr emptyDataset( GDALCreate( driver, mOutputFile.toUtf8(), columns, rows, 1, GDT_Float32, nullptr ) );
229 if ( !emptyDataset )
230 return false;
231
232 if ( GDALSetGeoTransform( emptyDataset.get(), geoTransform ) != CE_None )
233 return false;
234
235 // Set the projection on the raster destination to match the input layer
236 if ( GDALSetProjection( emptyDataset.get(), mSource->sourceCrs().toWkt().toLocal8Bit().data() ) != CE_None )
237 return false;
238
239 GDALRasterBandH poBand = GDALGetRasterBand( emptyDataset.get(), 1 );
240 if ( !poBand )
241 return false;
242
243 if ( GDALSetRasterNoDataValue( poBand, NO_DATA ) != CE_None )
244 return false;
245
246 float *line = static_cast< float * >( CPLMalloc( sizeof( float ) * columns ) );
247 for ( int i = 0; i < columns; i++ )
248 {
249 line[i] = NO_DATA;
250 }
251 // Write the empty raster
252 for ( int i = 0; i < rows ; i++ )
253 {
254 if ( GDALRasterIO( poBand, GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 ) != CE_None )
255 {
256 return false;
257 }
258 }
259
260 CPLFree( line );
261 return true;
262}
263
264double QgsKernelDensityEstimation::calculateKernelValue( const double distance, const double bandwidth, const QgsKernelDensityEstimation::KernelShape shape, const QgsKernelDensityEstimation::OutputValues outputType ) const
265{
266 switch ( shape )
267 {
268 case KernelTriangular:
269 return triangularKernel( distance, bandwidth, outputType );
270
271 case KernelUniform:
272 return uniformKernel( distance, bandwidth, outputType );
273
274 case KernelQuartic:
275 return quarticKernel( distance, bandwidth, outputType );
276
277 case KernelTriweight:
278 return triweightKernel( distance, bandwidth, outputType );
279
281 return epanechnikovKernel( distance, bandwidth, outputType );
282 }
283 return 0; //no warnings
284}
285
286/* The kernel functions below are taken from "Kernel Smoothing" by Wand and Jones (1995), p. 175
287 *
288 * Each kernel is multiplied by a normalizing constant "k", which normalizes the kernel area
289 * to 1 for a given bandwidth size.
290 *
291 * k is calculated by polar double integration of the kernel function
292 * between a radius of 0 to the specified bandwidth and equating the area to 1. */
293
294double QgsKernelDensityEstimation::uniformKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
295{
296 Q_UNUSED( distance )
297 switch ( outputType )
298 {
299 case OutputScaled:
300 {
301 // Normalizing constant
302 const double k = 2. / ( M_PI * bandwidth );
303
304 // Derived from Wand and Jones (1995), p. 175
305 return k * ( 0.5 / bandwidth );
306 }
307 case OutputRaw:
308 return 1.0;
309 }
310 return 0.0; // NO warnings!!!!!
311}
312
313double QgsKernelDensityEstimation::quarticKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
314{
315 switch ( outputType )
316 {
317 case OutputScaled:
318 {
319 // Normalizing constant
320 const double k = 116. / ( 5. * M_PI * std::pow( bandwidth, 2 ) );
321
322 // Derived from Wand and Jones (1995), p. 175
323 return k * ( 15. / 16. ) * std::pow( 1. - std::pow( distance / bandwidth, 2 ), 2 );
324 }
325 case OutputRaw:
326 return std::pow( 1. - std::pow( distance / bandwidth, 2 ), 2 );
327 }
328 return 0.0; //no, seriously, I told you NO WARNINGS!
329}
330
331double QgsKernelDensityEstimation::triweightKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
332{
333 switch ( outputType )
334 {
335 case OutputScaled:
336 {
337 // Normalizing constant
338 const double k = 128. / ( 35. * M_PI * std::pow( bandwidth, 2 ) );
339
340 // Derived from Wand and Jones (1995), p. 175
341 return k * ( 35. / 32. ) * std::pow( 1. - std::pow( distance / bandwidth, 2 ), 3 );
342 }
343 case OutputRaw:
344 return std::pow( 1. - std::pow( distance / bandwidth, 2 ), 3 );
345 }
346 return 0.0; // this is getting ridiculous... don't you ever listen to a word I say?
347}
348
349double QgsKernelDensityEstimation::epanechnikovKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
350{
351 switch ( outputType )
352 {
353 case OutputScaled:
354 {
355 // Normalizing constant
356 const double k = 8. / ( 3. * M_PI * std::pow( bandwidth, 2 ) );
357
358 // Derived from Wand and Jones (1995), p. 175
359 return k * ( 3. / 4. ) * ( 1. - std::pow( distance / bandwidth, 2 ) );
360 }
361 case OutputRaw:
362 return ( 1. - std::pow( distance / bandwidth, 2 ) );
363 }
364
365 return 0.0; // la la la i'm not listening
366}
367
368double QgsKernelDensityEstimation::triangularKernel( const double distance, const double bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
369{
370 switch ( outputType )
371 {
372 case OutputScaled:
373 {
374 // Normalizing constant. In this case it's calculated a little different
375 // due to the inclusion of the non-standard "decay" parameter
376
377 if ( mDecay >= 0 )
378 {
379 const double k = 3. / ( ( 1. + 2. * mDecay ) * M_PI * std::pow( bandwidth, 2 ) );
380
381 // Derived from Wand and Jones (1995), p. 175 (with addition of decay parameter)
382 return k * ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
383 }
384 else
385 {
386 // Non-standard or mathematically valid negative decay ("coolmap")
387 return ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
388 }
389 }
390 case OutputRaw:
391 return ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
392 }
393 return 0.0; // ....
394}
395
396QgsRectangle QgsKernelDensityEstimation::calculateBounds() const
397{
398 if ( !mSource )
399 return QgsRectangle();
400
401 QgsRectangle bbox = mSource->sourceExtent();
402
403 double radius = 0;
404 if ( mRadiusField >= 0 )
405 {
406 // if radius is using a field, find the max value
407 radius = mSource->maximumValue( mRadiusField ).toDouble();
408 }
409 else
410 {
411 radius = mRadius;
412 }
413 // expand the bounds by the maximum search radius
414 bbox.setXMinimum( bbox.xMinimum() - radius );
415 bbox.setYMinimum( bbox.yMinimum() - radius );
416 bbox.setXMaximum( bbox.xMaximum() + radius );
417 bbox.setYMaximum( bbox.yMaximum() + radius );
418 return bbox;
419}
420
421
422
423
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:208
Result
Result of operation.
Definition qgskde.h:61
@ DriverError
Could not open the driver for the specified format.
Definition qgskde.h:63
@ FileCreationError
Error creating output file.
Definition qgskde.h:65
@ RasterIoError
Error writing to raster.
Definition qgskde.h:66
@ Success
Operation completed successfully.
Definition qgskde.h:62
@ InvalidParameters
Input parameters were not valid.
Definition qgskde.h:64
OutputValues
Output values type.
Definition qgskde.h:54
@ OutputRaw
Output the raw KDE values.
Definition qgskde.h:55
@ OutputScaled
Output mathematically correct scaled values.
Definition qgskde.h:56
KernelShape
Kernel shape type.
Definition qgskde.h:44
@ KernelTriweight
Triweight kernel.
Definition qgskde.h:48
@ KernelUniform
Uniform (flat) kernel.
Definition qgskde.h:47
@ KernelEpanechnikov
Epanechnikov kernel.
Definition qgskde.h:49
@ KernelTriangular
Triangular kernel.
Definition qgskde.h:46
@ KernelQuartic
Quartic kernel.
Definition qgskde.h:45
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:79
QString weightField
Field name for weighting field, or empty if not using weights.
Definition qgskde.h:82