QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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
Looks up field&#39;s index from the field name.
Definition: qgsfields.cpp:324
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:342
virtual QgsRectangle sourceExtent() const
Returns the extent of all geometries from the source.
A rectangle specified with double values.
Definition: qgsrectangle.h:41
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.
void setXMaximum(double x)
Set the maximum x value.
Definition: qgsrectangle.h:135
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:111
QVector< QgsPointXY > QgsMultiPointXY
A collection of QgsPoints that share a common collection of attributes.
Definition: qgsgeometry.h:79
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:55
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:202
void setYMinimum(double y)
Set the minimum y value.
Definition: qgsrectangle.h:140
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
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:177
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:162
QgsPointXY asPoint() const
Returns the contents of the geometry as a 2-dimensional point.
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:145
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:436
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:167
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:172
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 the contents of the geometry as a multi-point.
QgsGeometry geometry
Definition: qgsfeature.h:67
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:262
Error creating output file.
Definition: qgskde.h:66
void setXMinimum(double x)
Set the minimum x value.
Definition: qgsrectangle.h:130
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:209