|
| 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 | + |
0 commit comments