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