QGIS API Documentation  3.2.0-Bonn (bc43194)
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  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  int rows = std::max( std::ceil( mBounds.height() / mPixelSize ) + 1, 1.0 );
87  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  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  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  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  unsigned int xPosition = ( ( ( *pointIt ).x() - mBounds.xMinimum() ) / mPixelSize ) - buffer;
161  unsigned int yPosition = ( ( ( *pointIt ).y() - mBounds.yMinimum() ) / mPixelSize ) - buffer;
162  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  double pixelCentroidX = ( xPosition + xp + 0.5 ) * mPixelSize + mBounds.xMinimum();
178  double pixelCentroidY = ( yPosition + yp + 0.5 ) * mPixelSize + mBounds.yMinimum();
179 
180  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  double pixelValue = weight * calculateKernelValue( distance, radius, mShape, mOutputValues );
189  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  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  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  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  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  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  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 
int lookupField(const QString &fieldName) const
Look up field&#39;s index from the field name.
Definition: qgsfields.cpp:299
QgsKernelDensityEstimation(const Parameters &parameters, const QString &outputFile, const QString &outputFormat)
Constructor for QgsKernelDensityEstimation.
Definition: qgskde.cpp:23
Wrapper for iterator of features from vector data provider or vector layer.
bool contains(const QgsRectangle &rect) const
Returns true when rectangle contains other rectangle.
Definition: qgsrectangle.h:335
virtual QgsRectangle sourceExtent() const
Returns the extent of all geometries from the source.
A rectangle specified with double values.
Definition: qgsrectangle.h:40
QString weightField
Field name for weighting field, or empty if not using weights.
Definition: qgskde.h:83
virtual QgsFields fields() const =0
Returns the fields associated with features in the source.
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
void setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:134
Uniform (flat) kernel.
Definition: qgskde.h:48
#define NO_DATA
Definition: qgskde.cpp:21
A class to represent a 2D point.
Definition: qgspointxy.h:43
Error writing to raster.
Definition: qgskde.h:67
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:104
QVector< QgsPointXY > QgsMultiPointXY
A collection of QgsPoints that share a common collection of attributes.
Definition: qgsgeometry.h:72
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:62
KernelShape
Kernel shape type.
Definition: qgskde.h:44
OutputValues
Output values type.
Definition: qgskde.h:54
Input parameters were not valid.
Definition: qgskde.h:65
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:201
void setYMinimum(double y)
Set the minimum y value.
Definition: qgsrectangle.h:139
Result finalise()
Finalises the output file.
Definition: qgskde.cpp:208
This class wraps a request for features to a vector layer (or directly its vector data provider)...
Output the raw KDE values.
Definition: qgskde.h:56
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:101
Could not open the driver for the specified format.
Definition: qgskde.h:64
Output mathematically correct scaled values.
Definition: qgskde.h:57
virtual QgsCoordinateReferenceSystem sourceCrs() const =0
Returns the coordinate reference system for features in the source.
Operation completed successfully.
Definition: qgskde.h:63
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:176
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:161
QgsPointXY asPoint() const
Returns contents of the geometry as a point if wkbType is WKBPoint, otherwise returns [0...
QString radiusField
Field for radius, or empty if using a fixed radius.
Definition: qgskde.h:80
void setYMaximum(double y)
Set the maximum y value.
Definition: qgsrectangle.h:144
Result addFeature(const QgsFeature &feature)
Adds a single feature to the KDE surface.
Definition: qgskde.cpp:107
QString toWkt() const
Returns a WKT representation of this CRS.
bool isNull() const
Test if the rectangle is null (all coordinates zero or after call to setMinimal()).
Definition: qgsrectangle.h:429
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:166
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:171
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:134
Result run()
Runs the KDE calculation across the whole layer at once.
Definition: qgskde.cpp:44
QgsMultiPointXY asMultiPoint() const
Returns contents of the geometry as a multi point if wkbType is WKBMultiPoint, otherwise an empty lis...
virtual QVariant maximumValue(int fieldIndex) const
Returns the maximum value for an attribute column or an invalid variant in case of error...
QList< int > QgsAttributeList
Definition: qgsfield.h:27
bool nextFeature(QgsFeature &f)
Result
Result of operation.
Definition: qgskde.h:61
Result prepare()
Prepares the output file for writing and setups up the surface calculation.
Definition: qgskde.cpp:69
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:255
Error creating output file.
Definition: qgskde.h:66
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:129
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:208