Skip to content

Commit f212746

Browse files
committedDec 6, 2016
Move heatmap generation code to analysis lib
And clean it up a lot
1 parent c558d51 commit f212746

File tree

5 files changed

+729
-0
lines changed

5 files changed

+729
-0
lines changed
 

‎python/analysis/analysis.sip

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@
4444
%Include raster/qgsderivativefilter.sip
4545
%Include raster/qgsaspectfilter.sip
4646
%Include raster/qgshillshadefilter.sip
47+
%Include raster/qgskde.sip
4748
%Include raster/qgsninecellfilter.sip
4849
%Include raster/qgsrastercalcnode.sip
4950
%Include raster/qgsrastercalculator.sip

‎python/analysis/raster/qgskde.sip

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
/**
2+
* \class QgsKernelDensityEstimation
3+
* \ingroup analysis
4+
* Performs Kernel Density Estimation ("heatmap") calculations on a vector layer.
5+
* @note added in QGIS 3.0
6+
*/
7+
class QgsKernelDensityEstimation
8+
{
9+
%TypeHeaderCode
10+
#include <qgskde.h>
11+
%End
12+
13+
public:
14+
15+
//! Kernel shape type
16+
enum KernelShape
17+
{
18+
KernelQuartic, //!< Quartic kernel
19+
KernelTriangular, //!< Triangular kernel
20+
KernelUniform, //!< Uniform (flat) kernel
21+
KernelTriweight, //!< Triweight kernel
22+
KernelEpanechnikov, //!< Epanechnikov kernel
23+
};
24+
25+
//! Output values type
26+
enum OutputValues
27+
{
28+
OutputRaw, //!< Output the raw KDE values
29+
OutputScaled, //!< Output mathematically correct scaled values
30+
};
31+
32+
//! Result of operation
33+
enum Result
34+
{
35+
Success, //!< Operation completed successfully
36+
DriverError, //!< Could not open the driver for the specified format
37+
InvalidParameters, //!< Input parameters were not valid
38+
FileCreationError, //!< Error creating output file
39+
RasterIoError, //!< Error writing to raster
40+
};
41+
42+
//! KDE parameters
43+
struct Parameters
44+
{
45+
//! Vector point layer
46+
QgsVectorLayer* vectorLayer;
47+
48+
//! Fixed radius, in map units
49+
double radius;
50+
51+
//! Field for radius, or empty if using a fixed radius
52+
QString radiusField;
53+
54+
//! Field name for weighting field, or empty if not using weights
55+
QString weightField;
56+
57+
//! Size of pixel in output file
58+
double pixelSize;
59+
60+
//! Kernel shape
61+
QgsKernelDensityEstimation::KernelShape shape;
62+
63+
//! Decay ratio (Triangular kernels only)
64+
double decayRatio;
65+
66+
//! Type of output value
67+
QgsKernelDensityEstimation::OutputValues outputValues;
68+
};
69+
70+
/**
71+
* Constructor for QgsKernelDensityEstimation. Requires a Parameters object specifying the options to use
72+
* to generate the surface. The output path and file format are also required.
73+
*/
74+
QgsKernelDensityEstimation( const Parameters& parameters, const QString& outputFile, const QString& outputFormat );
75+
76+
/**
77+
* Runs the KDE calculation across the whole layer at once. Either call this method, or manually
78+
* call run(), addFeature() and finalise() separately.
79+
*/
80+
Result run();
81+
82+
/**
83+
* Prepares the output file for writing and setups up the surface calculation. This must be called
84+
* before adding features via addFeature().
85+
* @see addFeature()
86+
* @see finalise()
87+
*/
88+
Result prepare();
89+
90+
/**
91+
* Adds a single feature to the KDE surface. prepare() must be called before adding features.
92+
* @see prepare()
93+
* @see finalise()
94+
*/
95+
Result addFeature( const QgsFeature& feature );
96+
97+
/**
98+
* Finalises the output file. Must be called after adding all features via addFeature().
99+
* @see prepare()
100+
* @see addFeature()
101+
*/
102+
Result finalise();
103+
104+
};

‎src/analysis/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ SET(QGIS_ANALYSIS_SRCS
2727
raster/qgsruggednessfilter.cpp
2828
raster/qgsderivativefilter.cpp
2929
raster/qgshillshadefilter.cpp
30+
raster/qgskde.cpp
3031
raster/qgsslopefilter.cpp
3132
raster/qgsaspectfilter.cpp
3233
raster/qgstotalcurvaturefilter.cpp
@@ -106,6 +107,7 @@ SET(QGIS_ANALYSIS_HDRS
106107
raster/qgsaspectfilter.h
107108
raster/qgsderivativefilter.h
108109
raster/qgshillshadefilter.h
110+
raster/qgskde.h
109111
raster/qgsninecellfilter.h
110112
raster/qgsrastercalculator.h
111113
raster/qgsrelief.h

‎src/analysis/raster/qgskde.cpp

Lines changed: 448 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,448 @@
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 "qgsvectorlayer.h"
18+
#include "qgsgeometry.h"
19+
20+
#define NO_DATA -9999
21+
22+
#ifndef M_PI
23+
#define M_PI 3.14159265358979323846
24+
#endif
25+
26+
#if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
27+
#define TO8F(x) (x).toUtf8().constData()
28+
#else
29+
#define TO8F(x) QFile::encodeName( x ).constData()
30+
#endif
31+
32+
QgsKernelDensityEstimation::QgsKernelDensityEstimation( const QgsKernelDensityEstimation::Parameters& parameters, const QString& outputFile, const QString& outputFormat )
33+
: mInputLayer( parameters.vectorLayer )
34+
, mOutputFile( outputFile )
35+
, mOutputFormat( outputFormat )
36+
, mRadiusField( -1 )
37+
, mWeightField( -1 )
38+
, mRadius( parameters.radius )
39+
, mPixelSize( parameters.pixelSize )
40+
, mShape( parameters.shape )
41+
, mDecay( parameters.decayRatio )
42+
, mOutputValues( parameters.outputValues )
43+
, mBufferSize( -1 )
44+
, mDatasetH( nullptr )
45+
, mRasterBandH( nullptr )
46+
{
47+
if ( !parameters.radiusField.isEmpty() )
48+
mRadiusField = mInputLayer->fields().lookupField( parameters.radiusField );
49+
if ( !parameters.weightField.isEmpty() )
50+
mWeightField = mInputLayer->fields().lookupField( parameters.weightField );
51+
}
52+
53+
QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::run()
54+
{
55+
Result result = prepare();
56+
if ( result != Success )
57+
return result;
58+
59+
QgsAttributeList requiredAttributes;
60+
61+
if ( mRadiusField >= 0 )
62+
requiredAttributes << mRadiusField;
63+
64+
if ( mWeightField >= 0 )
65+
requiredAttributes << mWeightField;
66+
67+
QgsFeatureIterator fit = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes ) );
68+
69+
QgsFeature f;
70+
while ( fit.nextFeature( f ) )
71+
{
72+
addFeature( f );
73+
}
74+
75+
return finalise();
76+
}
77+
78+
QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::prepare()
79+
{
80+
GDALAllRegister();
81+
82+
GDALDriverH driver = GDALGetDriverByName( mOutputFormat.toUtf8() );
83+
if ( !driver )
84+
{
85+
return DriverError;
86+
}
87+
88+
if ( !mInputLayer )
89+
return InvalidParameters;
90+
91+
mBounds = calculateBounds();
92+
if ( mBounds.isNull() )
93+
return InvalidParameters;
94+
95+
int rows = qMax( qRound( mBounds.height() / mPixelSize ), 1 );
96+
int cols = qMax( qRound( mBounds.width() / mPixelSize ), 1 );
97+
98+
if ( !createEmptyLayer( driver, mBounds, rows, cols ) )
99+
return FileCreationError;
100+
101+
// open the raster in GA_Update mode
102+
mDatasetH = GDALOpen( TO8F( mOutputFile ), GA_Update );
103+
if ( !mDatasetH )
104+
return FileCreationError;
105+
mRasterBandH = GDALGetRasterBand( mDatasetH, 1 );
106+
if ( !mRasterBandH )
107+
return FileCreationError;
108+
109+
mBufferSize = -1;
110+
if ( mRadiusField < 0 )
111+
mBufferSize = radiusSizeInPixels( mRadius );
112+
113+
return Success;
114+
}
115+
116+
QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::addFeature( const QgsFeature& feature )
117+
{
118+
QgsGeometry featureGeometry = feature.geometry();
119+
if ( featureGeometry.isEmpty() )
120+
{
121+
return Success;
122+
}
123+
124+
// convert the geometry to multipoint
125+
QgsMultiPoint multiPoints;
126+
if ( !featureGeometry.isMultipart() )
127+
{
128+
QgsPoint p = featureGeometry.asPoint();
129+
// avoiding any empty points or out of extent points
130+
if ( !mBounds.contains( p ) )
131+
return Success;
132+
133+
multiPoints << p;
134+
}
135+
else
136+
{
137+
multiPoints = featureGeometry.asMultiPoint();
138+
}
139+
140+
// if radius is variable then fetch it and calculate new pixel buffer size
141+
double radius = mRadius;
142+
int buffer = mBufferSize;
143+
if ( mRadiusField >= 0 )
144+
{
145+
radius = feature.attribute( mRadiusField ).toDouble();
146+
buffer = radiusSizeInPixels( radius );
147+
}
148+
int blockSize = 2 * buffer + 1; //Block SIDE would be more appropriate
149+
150+
// calculate weight
151+
double weight = 1.0;
152+
if ( mWeightField >= 0 )
153+
{
154+
weight = feature.attribute( mWeightField ).toDouble();
155+
}
156+
157+
Result result = Success;
158+
159+
//loop through all points in multipoint
160+
for ( QgsMultiPoint::const_iterator pointIt = multiPoints.constBegin(); pointIt != multiPoints.constEnd(); ++pointIt )
161+
{
162+
// avoiding any empty points or out of extent points
163+
if ( !mBounds.contains( *pointIt ) )
164+
{
165+
continue;
166+
}
167+
168+
// calculate the pixel position
169+
unsigned int xPosition = ((( *pointIt ).x() - mBounds.xMinimum() ) / mPixelSize ) - buffer;
170+
unsigned int yPosition = ((( *pointIt ).y() - mBounds.yMinimum() ) / mPixelSize ) - buffer;
171+
172+
// get the data
173+
float *dataBuffer = ( float * ) CPLMalloc( sizeof( float ) * blockSize * blockSize );
174+
if ( GDALRasterIO( mRasterBandH, GF_Read, xPosition, yPosition, blockSize, blockSize,
175+
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
176+
{
177+
result = RasterIoError;
178+
}
179+
180+
for ( int xp = 0; xp <= buffer; xp++ )
181+
{
182+
for ( int yp = 0; yp <= buffer; yp++ )
183+
{
184+
double distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
185+
186+
// is pixel outside search bandwidth of feature?
187+
if ( distance > buffer )
188+
{
189+
continue;
190+
}
191+
192+
double pixelValue = weight * calculateKernelValue( distance, buffer, mShape, mOutputValues );
193+
194+
// clearing anamolies along the axes
195+
if ( xp == 0 && yp == 0 )
196+
{
197+
pixelValue /= 4;
198+
}
199+
else if ( xp == 0 || yp == 0 )
200+
{
201+
pixelValue /= 2;
202+
}
203+
204+
int pos[4];
205+
pos[0] = ( buffer + xp ) * blockSize + ( buffer + yp );
206+
pos[1] = ( buffer + xp ) * blockSize + ( buffer - yp );
207+
pos[2] = ( buffer - xp ) * blockSize + ( buffer + yp );
208+
pos[3] = ( buffer - xp ) * blockSize + ( buffer - yp );
209+
for ( int p = 0; p < 4; p++ )
210+
{
211+
if ( dataBuffer[ pos[p] ] == NO_DATA )
212+
{
213+
dataBuffer[ pos[p] ] = 0;
214+
}
215+
dataBuffer[ pos[p] ] += pixelValue;
216+
}
217+
}
218+
}
219+
if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPosition, blockSize, blockSize,
220+
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
221+
{
222+
result = RasterIoError;
223+
}
224+
CPLFree( dataBuffer );
225+
}
226+
227+
return result;
228+
}
229+
230+
QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::finalise()
231+
{
232+
GDALClose(( GDALDatasetH ) mDatasetH );
233+
mDatasetH = nullptr;
234+
mRasterBandH = nullptr;
235+
return Success;
236+
}
237+
238+
int QgsKernelDensityEstimation::radiusSizeInPixels( double radius ) const
239+
{
240+
int buffer = radius / mPixelSize;
241+
if ( radius - ( mPixelSize * buffer ) > 0.5 )
242+
{
243+
++buffer;
244+
}
245+
return buffer;
246+
}
247+
248+
bool QgsKernelDensityEstimation::createEmptyLayer( GDALDriverH driver, const QgsRectangle& bounds, int rows, int columns ) const
249+
{
250+
double geoTransform[6] = { bounds.xMinimum(), mPixelSize, 0, bounds.yMinimum(), 0, mPixelSize };
251+
GDALDatasetH emptyDataset = GDALCreate( driver, mOutputFile.toUtf8(), columns, rows, 1, GDT_Float32, nullptr );
252+
if ( !emptyDataset )
253+
return false;
254+
255+
if ( GDALSetGeoTransform( emptyDataset, geoTransform ) != CE_None )
256+
return false;
257+
258+
// Set the projection on the raster destination to match the input layer
259+
if ( GDALSetProjection( emptyDataset, mInputLayer->crs().toWkt().toLocal8Bit().data() ) != CE_None )
260+
return false;
261+
262+
GDALRasterBandH poBand = GDALGetRasterBand( emptyDataset, 1 );
263+
if ( !poBand )
264+
return false;
265+
266+
if ( GDALSetRasterNoDataValue( poBand, NO_DATA ) != CE_None )
267+
return false;
268+
269+
float* line = static_cast< float * >( CPLMalloc( sizeof( float ) * columns ) );
270+
for ( int i = 0; i < columns; i++ )
271+
{
272+
line[i] = NO_DATA;
273+
}
274+
// Write the empty raster
275+
for ( int i = 0; i < rows ; i++ )
276+
{
277+
if ( GDALRasterIO( poBand, GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 ) != CE_None )
278+
{
279+
return false;
280+
}
281+
}
282+
283+
CPLFree( line );
284+
//close the dataset
285+
GDALClose( emptyDataset );
286+
return true;
287+
}
288+
289+
double QgsKernelDensityEstimation::calculateKernelValue( const double distance, const int bandwidth, const QgsKernelDensityEstimation::KernelShape shape, const QgsKernelDensityEstimation::OutputValues outputType ) const
290+
{
291+
switch ( shape )
292+
{
293+
case KernelTriangular:
294+
return triangularKernel( distance, bandwidth, outputType );
295+
296+
case KernelUniform:
297+
return uniformKernel( distance, bandwidth, outputType );
298+
299+
case KernelQuartic:
300+
return quarticKernel( distance, bandwidth, outputType );
301+
302+
case KernelTriweight:
303+
return triweightKernel( distance, bandwidth, outputType );
304+
305+
case KernelEpanechnikov:
306+
return epanechnikovKernel( distance, bandwidth, outputType );
307+
}
308+
return 0; //no warnings
309+
}
310+
311+
/* The kernel functions below are taken from "Kernel Smoothing" by Wand and Jones (1995), p. 175
312+
*
313+
* Each kernel is multiplied by a normalizing constant "k", which normalizes the kernel area
314+
* to 1 for a given bandwidth size.
315+
*
316+
* k is calculated by polar double integration of the kernel function
317+
* between a radius of 0 to the specified bandwidth and equating the area to 1. */
318+
319+
double QgsKernelDensityEstimation::uniformKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
320+
{
321+
Q_UNUSED( distance );
322+
switch ( outputType )
323+
{
324+
case OutputScaled:
325+
{
326+
// Normalizing constant
327+
double k = 2. / ( M_PI * ( double )bandwidth );
328+
329+
// Derived from Wand and Jones (1995), p. 175
330+
return k * ( 0.5 / ( double )bandwidth );
331+
}
332+
case OutputRaw:
333+
return 1.0;
334+
}
335+
return 0.0; // NO warnings!!!!!
336+
}
337+
338+
double QgsKernelDensityEstimation::quarticKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
339+
{
340+
switch ( outputType )
341+
{
342+
case OutputScaled:
343+
{
344+
// Normalizing constant
345+
double k = 116. / ( 5. * M_PI * pow(( double )bandwidth, 2 ) );
346+
347+
// Derived from Wand and Jones (1995), p. 175
348+
return k * ( 15. / 16. ) * pow( 1. - pow( distance / ( double )bandwidth, 2 ), 2 );
349+
}
350+
case OutputRaw:
351+
return pow( 1. - pow( distance / ( double )bandwidth, 2 ), 2 );
352+
}
353+
return 0.0; //no, seriously, I told you NO WARNINGS!
354+
}
355+
356+
double QgsKernelDensityEstimation::triweightKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
357+
{
358+
switch ( outputType )
359+
{
360+
case OutputScaled:
361+
{
362+
// Normalizing constant
363+
double k = 128. / ( 35. * M_PI * pow(( double )bandwidth, 2 ) );
364+
365+
// Derived from Wand and Jones (1995), p. 175
366+
return k * ( 35. / 32. ) * pow( 1. - pow( distance / ( double )bandwidth, 2 ), 3 );
367+
}
368+
case OutputRaw:
369+
return pow( 1. - pow( distance / ( double )bandwidth, 2 ), 3 );
370+
}
371+
return 0.0; // this is getting ridiculous... don't you ever listen to a word I say?
372+
}
373+
374+
double QgsKernelDensityEstimation::epanechnikovKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
375+
{
376+
switch ( outputType )
377+
{
378+
case OutputScaled:
379+
{
380+
// Normalizing constant
381+
double k = 8. / ( 3. * M_PI * pow(( double )bandwidth, 2 ) );
382+
383+
// Derived from Wand and Jones (1995), p. 175
384+
return k * ( 3. / 4. ) * ( 1. - pow( distance / ( double )bandwidth, 2 ) );
385+
}
386+
case OutputRaw:
387+
return ( 1. - pow( distance / ( double )bandwidth, 2 ) );
388+
}
389+
390+
return 0.0; // la la la i'm not listening
391+
}
392+
393+
double QgsKernelDensityEstimation::triangularKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
394+
{
395+
switch ( outputType )
396+
{
397+
case OutputScaled:
398+
{
399+
// Normalizing constant. In this case it's calculated a little different
400+
// due to the inclusion of the non-standard "decay" parameter
401+
402+
if ( mDecay >= 0 )
403+
{
404+
double k = 3. / (( 1. + 2. * mDecay ) * M_PI * pow(( double )bandwidth, 2 ) );
405+
406+
// Derived from Wand and Jones (1995), p. 175 (with addition of decay parameter)
407+
return k * ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) );
408+
}
409+
else
410+
{
411+
// Non-standard or mathematically valid negative decay ("coolmap")
412+
return ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) );
413+
}
414+
}
415+
case OutputRaw:
416+
return ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) );
417+
}
418+
return 0.0; // ....
419+
}
420+
421+
QgsRectangle QgsKernelDensityEstimation::calculateBounds() const
422+
{
423+
if ( !mInputLayer )
424+
return QgsRectangle();
425+
426+
QgsRectangle bbox = mInputLayer->extent();
427+
428+
double radius = 0;
429+
if ( mRadiusField >= 0 )
430+
{
431+
// if radius is using a field, find the max value
432+
radius = mInputLayer->maximumValue( mRadiusField ).toDouble();
433+
}
434+
else
435+
{
436+
radius = mRadius;
437+
}
438+
// expand the bounds by the maximum search radius
439+
bbox.setXMinimum( bbox.xMinimum() - radius );
440+
bbox.setYMinimum( bbox.yMinimum() - radius );
441+
bbox.setXMaximum( bbox.xMaximum() + radius );
442+
bbox.setYMaximum( bbox.yMaximum() + radius );
443+
return bbox;
444+
}
445+
446+
447+
448+

‎src/analysis/raster/qgskde.h

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
/***************************************************************************
2+
qgskde.h
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+
#ifndef QGSKDE_H
17+
#define QGSKDE_H
18+
19+
#include "qgsrectangle.h"
20+
#include <QString>
21+
22+
// GDAL includes
23+
#include <gdal.h>
24+
#include <cpl_string.h>
25+
#include <cpl_conv.h>
26+
27+
class QgsVectorLayer;
28+
class QProgressDialog;
29+
class QgsFeature;
30+
31+
32+
/**
33+
* \class QgsKernelDensityEstimation
34+
* \ingroup analysis
35+
* Performs Kernel Density Estimation ("heatmap") calculations on a vector layer.
36+
* @note added in QGIS 3.0
37+
*/
38+
class ANALYSIS_EXPORT QgsKernelDensityEstimation
39+
{
40+
public:
41+
42+
//! Kernel shape type
43+
enum KernelShape
44+
{
45+
KernelQuartic = 0, //!< Quartic kernel
46+
KernelTriangular, //!< Triangular kernel
47+
KernelUniform, //!< Uniform (flat) kernel
48+
KernelTriweight, //!< Triweight kernel
49+
KernelEpanechnikov, //!< Epanechnikov kernel
50+
};
51+
52+
//! Output values type
53+
enum OutputValues
54+
{
55+
OutputRaw = 0, //!< Output the raw KDE values
56+
OutputScaled, //!< Output mathematically correct scaled values
57+
};
58+
59+
//! Result of operation
60+
enum Result
61+
{
62+
Success, //!< Operation completed successfully
63+
DriverError, //!< Could not open the driver for the specified format
64+
InvalidParameters, //!< Input parameters were not valid
65+
FileCreationError, //!< Error creating output file
66+
RasterIoError, //!< Error writing to raster
67+
};
68+
69+
//! KDE parameters
70+
struct Parameters
71+
{
72+
//! Vector point layer
73+
QgsVectorLayer* vectorLayer;
74+
75+
//! Fixed radius, in map units
76+
double radius;
77+
78+
//! Field for radius, or empty if using a fixed radius
79+
QString radiusField;
80+
81+
//! Field name for weighting field, or empty if not using weights
82+
QString weightField;
83+
84+
//! Size of pixel in output file
85+
double pixelSize;
86+
87+
//! Kernel shape
88+
KernelShape shape;
89+
90+
//! Decay ratio (Triangular kernels only)
91+
double decayRatio;
92+
93+
//! Type of output value
94+
OutputValues outputValues;
95+
};
96+
97+
/**
98+
* Constructor for QgsKernelDensityEstimation. Requires a Parameters object specifying the options to use
99+
* to generate the surface. The output path and file format are also required.
100+
*/
101+
QgsKernelDensityEstimation( const Parameters& parameters, const QString& outputFile, const QString& outputFormat );
102+
103+
/**
104+
* Runs the KDE calculation across the whole layer at once. Either call this method, or manually
105+
* call run(), addFeature() and finalise() separately.
106+
*/
107+
Result run();
108+
109+
/**
110+
* Prepares the output file for writing and setups up the surface calculation. This must be called
111+
* before adding features via addFeature().
112+
* @see addFeature()
113+
* @see finalise()
114+
*/
115+
Result prepare();
116+
117+
/**
118+
* Adds a single feature to the KDE surface. prepare() must be called before adding features.
119+
* @see prepare()
120+
* @see finalise()
121+
*/
122+
Result addFeature( const QgsFeature& feature );
123+
124+
/**
125+
* Finalises the output file. Must be called after adding all features via addFeature().
126+
* @see prepare()
127+
* @see addFeature()
128+
*/
129+
Result finalise();
130+
131+
private:
132+
133+
//! Calculate the value given to a point width a given distance for a specified kernel shape
134+
double calculateKernelValue( const double distance, const int bandwidth, const KernelShape shape, const OutputValues outputType ) const;
135+
//! Uniform kernel function
136+
double uniformKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
137+
//! Quartic kernel function
138+
double quarticKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
139+
//! Triweight kernel function
140+
double triweightKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
141+
//! Epanechnikov kernel function
142+
double epanechnikovKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
143+
//! Triangular kernel function
144+
double triangularKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
145+
146+
QgsRectangle calculateBounds() const;
147+
148+
QgsVectorLayer* mInputLayer;
149+
150+
QString mOutputFile;
151+
QString mOutputFormat;
152+
153+
int mRadiusField;
154+
int mWeightField;
155+
double mRadius;
156+
double mPixelSize;
157+
QgsRectangle mBounds;
158+
159+
KernelShape mShape;
160+
double mDecay;
161+
OutputValues mOutputValues;
162+
163+
int mBufferSize;
164+
165+
GDALDatasetH mDatasetH;
166+
GDALRasterBandH mRasterBandH;
167+
168+
//! Creates a new raster layer and initialises it to the no data value
169+
bool createEmptyLayer( GDALDriverH driver, const QgsRectangle& bounds, int rows, int columns ) const;
170+
int radiusSizeInPixels( double radius ) const;
171+
};
172+
173+
174+
#endif // QGSKDE_H

0 commit comments

Comments
 (0)
Please sign in to comment.