QGIS API Documentation  3.27.0-Master (0e23467727)
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 
23 QgsKernelDensityEstimation::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 
215 int 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 
225 bool 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 
264 double 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 
280  case KernelEpanechnikov:
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 
294 double 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 
313 double 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 
331 double 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 
349 double 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 
368 double 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 
396 QgsRectangle 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(WktVariant variant=WKT1_GDAL, 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)
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:56
QgsGeometry geometry
Definition: qgsfeature.h:67
QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Definition: qgsfeature.cpp:320
int lookupField(const QString &fieldName) const
Looks up field's index from the field name.
Definition: qgsfields.cpp:349
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:125
QgsMultiPointXY asMultiPoint() const
Returns the contents of the geometry as a multi-point.
Q_GADGET bool isNull
Definition: qgsgeometry.h:127
bool isMultipart() const SIP_HOLDGIL
Returns true if WKB of the geometry is of WKBMulti* type.
QgsPointXY asPoint() const
Returns the contents of the geometry as a 2-dimensional point.
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:62
@ DriverError
Could not open the driver for the specified format.
Definition: qgskde.h:64
@ FileCreationError
Error creating output file.
Definition: qgskde.h:66
@ RasterIoError
Error writing to raster.
Definition: qgskde.h:67
@ Success
Operation completed successfully.
Definition: qgskde.h:63
@ InvalidParameters
Input parameters were not valid.
Definition: qgskde.h:65
OutputValues
Output values type.
Definition: qgskde.h:55
@ OutputRaw
Output the raw KDE values.
Definition: qgskde.h:56
@ OutputScaled
Output mathematically correct scaled values.
Definition: qgskde.h:57
KernelShape
Kernel shape type.
Definition: qgskde.h:45
@ KernelTriweight
Triweight kernel.
Definition: qgskde.h:49
@ KernelUniform
Uniform (flat) kernel.
Definition: qgskde.h:48
@ KernelEpanechnikov
Epanechnikov kernel.
Definition: qgskde.h:50
@ KernelTriangular
Triangular kernel.
Definition: qgskde.h:47
@ KernelQuartic
Quartic kernel.
Definition: qgskde.h:46
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:59
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:183
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
void setYMinimum(double y) SIP_HOLDGIL
Set the minimum y value.
Definition: qgsrectangle.h:161
bool isNull() const
Test if the rectangle is null (all coordinates zero or after call to setMinimal()).
Definition: qgsrectangle.h:479
void setXMaximum(double x) SIP_HOLDGIL
Set the maximum x value.
Definition: qgsrectangle.h:156
void setXMinimum(double x) SIP_HOLDGIL
Set the minimum x value.
Definition: qgsrectangle.h:151
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
void setYMaximum(double y) SIP_HOLDGIL
Set the maximum y value.
Definition: qgsrectangle.h:166
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
bool contains(const QgsRectangle &rect) const SIP_HOLDGIL
Returns true when rectangle contains other rectangle.
Definition: qgsrectangle.h:363
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:139
QList< int > QgsAttributeList
Definition: qgsfield.h:26
QVector< QgsPointXY > QgsMultiPointXY
A collection of QgsPoints that share a common collection of attributes.
Definition: qgsgeometry.h:82
#define NO_DATA
Definition: qgskde.cpp:21
QString radiusField
Field for radius, or empty if using a fixed radius.
Definition: qgskde.h:80
QString weightField
Field name for weighting field, or empty if not using weights.
Definition: qgskde.h:83