diff --git a/python/analysis/analysis.sip b/python/analysis/analysis.sip index 636eae955a9..8b5c2b9cecc 100644 --- a/python/analysis/analysis.sip +++ b/python/analysis/analysis.sip @@ -44,6 +44,7 @@ %Include raster/qgsderivativefilter.sip %Include raster/qgsaspectfilter.sip %Include raster/qgshillshadefilter.sip +%Include raster/qgskde.sip %Include raster/qgsninecellfilter.sip %Include raster/qgsrastercalcnode.sip %Include raster/qgsrastercalculator.sip diff --git a/python/analysis/raster/qgskde.sip b/python/analysis/raster/qgskde.sip new file mode 100644 index 00000000000..f81e4852f79 --- /dev/null +++ b/python/analysis/raster/qgskde.sip @@ -0,0 +1,104 @@ +/** + * \class QgsKernelDensityEstimation + * \ingroup analysis + * Performs Kernel Density Estimation ("heatmap") calculations on a vector layer. + * @note added in QGIS 3.0 + */ +class QgsKernelDensityEstimation +{ +%TypeHeaderCode +#include +%End + + public: + + //! Kernel shape type + enum KernelShape + { + KernelQuartic, //!< Quartic kernel + KernelTriangular, //!< Triangular kernel + KernelUniform, //!< Uniform (flat) kernel + KernelTriweight, //!< Triweight kernel + KernelEpanechnikov, //!< Epanechnikov kernel + }; + + //! Output values type + enum OutputValues + { + OutputRaw, //!< Output the raw KDE values + OutputScaled, //!< Output mathematically correct scaled values + }; + + //! Result of operation + enum Result + { + Success, //!< Operation completed successfully + DriverError, //!< Could not open the driver for the specified format + InvalidParameters, //!< Input parameters were not valid + FileCreationError, //!< Error creating output file + RasterIoError, //!< Error writing to raster + }; + + //! KDE parameters + struct Parameters + { + //! Vector point layer + QgsVectorLayer* vectorLayer; + + //! Fixed radius, in map units + double radius; + + //! Field for radius, or empty if using a fixed radius + QString radiusField; + + //! Field name for weighting field, or empty if not using weights + QString weightField; + + //! Size of pixel in output file + double pixelSize; + + //! Kernel shape + QgsKernelDensityEstimation::KernelShape shape; + + //! Decay ratio (Triangular kernels only) + double decayRatio; + + //! Type of output value + QgsKernelDensityEstimation::OutputValues outputValues; + }; + + /** + * Constructor for QgsKernelDensityEstimation. Requires a Parameters object specifying the options to use + * to generate the surface. The output path and file format are also required. + */ + QgsKernelDensityEstimation( const Parameters& parameters, const QString& outputFile, const QString& outputFormat ); + + /** + * Runs the KDE calculation across the whole layer at once. Either call this method, or manually + * call run(), addFeature() and finalise() separately. + */ + Result run(); + + /** + * Prepares the output file for writing and setups up the surface calculation. This must be called + * before adding features via addFeature(). + * @see addFeature() + * @see finalise() + */ + Result prepare(); + + /** + * Adds a single feature to the KDE surface. prepare() must be called before adding features. + * @see prepare() + * @see finalise() + */ + Result addFeature( const QgsFeature& feature ); + + /** + * Finalises the output file. Must be called after adding all features via addFeature(). + * @see prepare() + * @see addFeature() + */ + Result finalise(); + +}; diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index ff77465170a..062cdb6aa15 100644 --- a/src/analysis/CMakeLists.txt +++ b/src/analysis/CMakeLists.txt @@ -27,6 +27,7 @@ SET(QGIS_ANALYSIS_SRCS raster/qgsruggednessfilter.cpp raster/qgsderivativefilter.cpp raster/qgshillshadefilter.cpp + raster/qgskde.cpp raster/qgsslopefilter.cpp raster/qgsaspectfilter.cpp raster/qgstotalcurvaturefilter.cpp @@ -106,6 +107,7 @@ SET(QGIS_ANALYSIS_HDRS raster/qgsaspectfilter.h raster/qgsderivativefilter.h raster/qgshillshadefilter.h + raster/qgskde.h raster/qgsninecellfilter.h raster/qgsrastercalculator.h raster/qgsrelief.h diff --git a/src/analysis/raster/qgskde.cpp b/src/analysis/raster/qgskde.cpp new file mode 100644 index 00000000000..97d36b456a3 --- /dev/null +++ b/src/analysis/raster/qgskde.cpp @@ -0,0 +1,448 @@ +/*************************************************************************** + qgskde.cpp + ---------- + Date : October 2016 + Copyright : (C) 2016 by Nyall Dawson + Email : nyall dot dawson at gmail dot com + *************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#include "qgskde.h" +#include "qgsvectorlayer.h" +#include "qgsgeometry.h" + +#define NO_DATA -9999 + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800 +#define TO8F(x) (x).toUtf8().constData() +#else +#define TO8F(x) QFile::encodeName( x ).constData() +#endif + +QgsKernelDensityEstimation::QgsKernelDensityEstimation( const QgsKernelDensityEstimation::Parameters& parameters, const QString& outputFile, const QString& outputFormat ) + : mInputLayer( parameters.vectorLayer ) + , mOutputFile( outputFile ) + , mOutputFormat( outputFormat ) + , mRadiusField( -1 ) + , mWeightField( -1 ) + , mRadius( parameters.radius ) + , mPixelSize( parameters.pixelSize ) + , mShape( parameters.shape ) + , mDecay( parameters.decayRatio ) + , mOutputValues( parameters.outputValues ) + , mBufferSize( -1 ) + , mDatasetH( nullptr ) + , mRasterBandH( nullptr ) +{ + if ( !parameters.radiusField.isEmpty() ) + mRadiusField = mInputLayer->fields().lookupField( parameters.radiusField ); + if ( !parameters.weightField.isEmpty() ) + mWeightField = mInputLayer->fields().lookupField( parameters.weightField ); +} + +QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::run() +{ + Result result = prepare(); + if ( result != Success ) + return result; + + QgsAttributeList requiredAttributes; + + if ( mRadiusField >= 0 ) + requiredAttributes << mRadiusField; + + if ( mWeightField >= 0 ) + requiredAttributes << mWeightField; + + QgsFeatureIterator fit = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes ) ); + + QgsFeature f; + while ( fit.nextFeature( f ) ) + { + addFeature( f ); + } + + return finalise(); +} + +QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::prepare() +{ + GDALAllRegister(); + + GDALDriverH driver = GDALGetDriverByName( mOutputFormat.toUtf8() ); + if ( !driver ) + { + return DriverError; + } + + if ( !mInputLayer ) + return InvalidParameters; + + mBounds = calculateBounds(); + if ( mBounds.isNull() ) + return InvalidParameters; + + int rows = qMax( qRound( mBounds.height() / mPixelSize ), 1 ); + int cols = qMax( qRound( mBounds.width() / mPixelSize ), 1 ); + + if ( !createEmptyLayer( driver, mBounds, rows, cols ) ) + return FileCreationError; + + // open the raster in GA_Update mode + mDatasetH = GDALOpen( TO8F( mOutputFile ), GA_Update ); + if ( !mDatasetH ) + return FileCreationError; + mRasterBandH = GDALGetRasterBand( mDatasetH, 1 ); + if ( !mRasterBandH ) + return FileCreationError; + + mBufferSize = -1; + if ( mRadiusField < 0 ) + mBufferSize = radiusSizeInPixels( mRadius ); + + return Success; +} + +QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::addFeature( const QgsFeature& feature ) +{ + QgsGeometry featureGeometry = feature.geometry(); + if ( featureGeometry.isEmpty() ) + { + return Success; + } + + // convert the geometry to multipoint + QgsMultiPoint multiPoints; + if ( !featureGeometry.isMultipart() ) + { + QgsPoint p = featureGeometry.asPoint(); + // avoiding any empty points or out of extent points + if ( !mBounds.contains( p ) ) + return Success; + + multiPoints << p; + } + else + { + multiPoints = featureGeometry.asMultiPoint(); + } + + // if radius is variable then fetch it and calculate new pixel buffer size + double radius = mRadius; + int buffer = mBufferSize; + if ( mRadiusField >= 0 ) + { + radius = feature.attribute( mRadiusField ).toDouble(); + buffer = radiusSizeInPixels( radius ); + } + int blockSize = 2 * buffer + 1; //Block SIDE would be more appropriate + + // calculate weight + double weight = 1.0; + if ( mWeightField >= 0 ) + { + weight = feature.attribute( mWeightField ).toDouble(); + } + + Result result = Success; + + //loop through all points in multipoint + for ( QgsMultiPoint::const_iterator pointIt = multiPoints.constBegin(); pointIt != multiPoints.constEnd(); ++pointIt ) + { + // avoiding any empty points or out of extent points + if ( !mBounds.contains( *pointIt ) ) + { + continue; + } + + // calculate the pixel position + unsigned int xPosition = ((( *pointIt ).x() - mBounds.xMinimum() ) / mPixelSize ) - buffer; + unsigned int yPosition = ((( *pointIt ).y() - mBounds.yMinimum() ) / mPixelSize ) - buffer; + + // get the data + float *dataBuffer = ( float * ) CPLMalloc( sizeof( float ) * blockSize * blockSize ); + if ( GDALRasterIO( mRasterBandH, GF_Read, xPosition, yPosition, blockSize, blockSize, + dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None ) + { + result = RasterIoError; + } + + for ( int xp = 0; xp <= buffer; xp++ ) + { + for ( int yp = 0; yp <= buffer; yp++ ) + { + double distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) ); + + // is pixel outside search bandwidth of feature? + if ( distance > buffer ) + { + continue; + } + + double pixelValue = weight * calculateKernelValue( distance, buffer, mShape, mOutputValues ); + + // clearing anamolies along the axes + if ( xp == 0 && yp == 0 ) + { + pixelValue /= 4; + } + else if ( xp == 0 || yp == 0 ) + { + pixelValue /= 2; + } + + int pos[4]; + pos[0] = ( buffer + xp ) * blockSize + ( buffer + yp ); + pos[1] = ( buffer + xp ) * blockSize + ( buffer - yp ); + pos[2] = ( buffer - xp ) * blockSize + ( buffer + yp ); + pos[3] = ( buffer - xp ) * blockSize + ( buffer - yp ); + for ( int p = 0; p < 4; p++ ) + { + if ( dataBuffer[ pos[p] ] == NO_DATA ) + { + dataBuffer[ pos[p] ] = 0; + } + dataBuffer[ pos[p] ] += pixelValue; + } + } + } + if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPosition, blockSize, blockSize, + dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None ) + { + result = RasterIoError; + } + CPLFree( dataBuffer ); + } + + return result; +} + +QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::finalise() +{ + GDALClose(( GDALDatasetH ) mDatasetH ); + mDatasetH = nullptr; + mRasterBandH = nullptr; + return Success; +} + +int QgsKernelDensityEstimation::radiusSizeInPixels( double radius ) const +{ + int buffer = radius / mPixelSize; + if ( radius - ( mPixelSize * buffer ) > 0.5 ) + { + ++buffer; + } + return buffer; +} + +bool QgsKernelDensityEstimation::createEmptyLayer( GDALDriverH driver, const QgsRectangle& bounds, int rows, int columns ) const +{ + double geoTransform[6] = { bounds.xMinimum(), mPixelSize, 0, bounds.yMinimum(), 0, mPixelSize }; + GDALDatasetH emptyDataset = GDALCreate( driver, mOutputFile.toUtf8(), columns, rows, 1, GDT_Float32, nullptr ); + if ( !emptyDataset ) + return false; + + if ( GDALSetGeoTransform( emptyDataset, geoTransform ) != CE_None ) + return false; + + // Set the projection on the raster destination to match the input layer + if ( GDALSetProjection( emptyDataset, mInputLayer->crs().toWkt().toLocal8Bit().data() ) != CE_None ) + return false; + + GDALRasterBandH poBand = GDALGetRasterBand( emptyDataset, 1 ); + if ( !poBand ) + return false; + + if ( GDALSetRasterNoDataValue( poBand, NO_DATA ) != CE_None ) + return false; + + float* line = static_cast< float * >( CPLMalloc( sizeof( float ) * columns ) ); + for ( int i = 0; i < columns; i++ ) + { + line[i] = NO_DATA; + } + // Write the empty raster + for ( int i = 0; i < rows ; i++ ) + { + if ( GDALRasterIO( poBand, GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 ) != CE_None ) + { + return false; + } + } + + CPLFree( line ); + //close the dataset + GDALClose( emptyDataset ); + return true; +} + +double QgsKernelDensityEstimation::calculateKernelValue( const double distance, const int bandwidth, const QgsKernelDensityEstimation::KernelShape shape, const QgsKernelDensityEstimation::OutputValues outputType ) const +{ + switch ( shape ) + { + case KernelTriangular: + return triangularKernel( distance, bandwidth, outputType ); + + case KernelUniform: + return uniformKernel( distance, bandwidth, outputType ); + + case KernelQuartic: + return quarticKernel( distance, bandwidth, outputType ); + + case KernelTriweight: + return triweightKernel( distance, bandwidth, outputType ); + + case KernelEpanechnikov: + return epanechnikovKernel( distance, bandwidth, outputType ); + } + return 0; //no warnings +} + +/* The kernel functions below are taken from "Kernel Smoothing" by Wand and Jones (1995), p. 175 + * + * Each kernel is multiplied by a normalizing constant "k", which normalizes the kernel area + * to 1 for a given bandwidth size. + * + * k is calculated by polar double integration of the kernel function + * between a radius of 0 to the specified bandwidth and equating the area to 1. */ + +double QgsKernelDensityEstimation::uniformKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const +{ + Q_UNUSED( distance ); + switch ( outputType ) + { + case OutputScaled: + { + // Normalizing constant + double k = 2. / ( M_PI * ( double )bandwidth ); + + // Derived from Wand and Jones (1995), p. 175 + return k * ( 0.5 / ( double )bandwidth ); + } + case OutputRaw: + return 1.0; + } + return 0.0; // NO warnings!!!!! +} + +double QgsKernelDensityEstimation::quarticKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const +{ + switch ( outputType ) + { + case OutputScaled: + { + // Normalizing constant + double k = 116. / ( 5. * M_PI * pow(( double )bandwidth, 2 ) ); + + // Derived from Wand and Jones (1995), p. 175 + return k * ( 15. / 16. ) * pow( 1. - pow( distance / ( double )bandwidth, 2 ), 2 ); + } + case OutputRaw: + return pow( 1. - pow( distance / ( double )bandwidth, 2 ), 2 ); + } + return 0.0; //no, seriously, I told you NO WARNINGS! +} + +double QgsKernelDensityEstimation::triweightKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const +{ + switch ( outputType ) + { + case OutputScaled: + { + // Normalizing constant + double k = 128. / ( 35. * M_PI * pow(( double )bandwidth, 2 ) ); + + // Derived from Wand and Jones (1995), p. 175 + return k * ( 35. / 32. ) * pow( 1. - pow( distance / ( double )bandwidth, 2 ), 3 ); + } + case OutputRaw: + return pow( 1. - pow( distance / ( double )bandwidth, 2 ), 3 ); + } + return 0.0; // this is getting ridiculous... don't you ever listen to a word I say? +} + +double QgsKernelDensityEstimation::epanechnikovKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const +{ + switch ( outputType ) + { + case OutputScaled: + { + // Normalizing constant + double k = 8. / ( 3. * M_PI * pow(( double )bandwidth, 2 ) ); + + // Derived from Wand and Jones (1995), p. 175 + return k * ( 3. / 4. ) * ( 1. - pow( distance / ( double )bandwidth, 2 ) ); + } + case OutputRaw: + return ( 1. - pow( distance / ( double )bandwidth, 2 ) ); + } + + return 0.0; // la la la i'm not listening +} + +double QgsKernelDensityEstimation::triangularKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const +{ + switch ( outputType ) + { + case OutputScaled: + { + // Normalizing constant. In this case it's calculated a little different + // due to the inclusion of the non-standard "decay" parameter + + if ( mDecay >= 0 ) + { + double k = 3. / (( 1. + 2. * mDecay ) * M_PI * pow(( double )bandwidth, 2 ) ); + + // Derived from Wand and Jones (1995), p. 175 (with addition of decay parameter) + return k * ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) ); + } + else + { + // Non-standard or mathematically valid negative decay ("coolmap") + return ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) ); + } + } + case OutputRaw: + return ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) ); + } + return 0.0; // .... +} + +QgsRectangle QgsKernelDensityEstimation::calculateBounds() const +{ + if ( !mInputLayer ) + return QgsRectangle(); + + QgsRectangle bbox = mInputLayer->extent(); + + double radius = 0; + if ( mRadiusField >= 0 ) + { + // if radius is using a field, find the max value + radius = mInputLayer->maximumValue( mRadiusField ).toDouble(); + } + else + { + radius = mRadius; + } + // expand the bounds by the maximum search radius + bbox.setXMinimum( bbox.xMinimum() - radius ); + bbox.setYMinimum( bbox.yMinimum() - radius ); + bbox.setXMaximum( bbox.xMaximum() + radius ); + bbox.setYMaximum( bbox.yMaximum() + radius ); + return bbox; +} + + + + diff --git a/src/analysis/raster/qgskde.h b/src/analysis/raster/qgskde.h new file mode 100644 index 00000000000..76605c66583 --- /dev/null +++ b/src/analysis/raster/qgskde.h @@ -0,0 +1,174 @@ +/*************************************************************************** + qgskde.h + -------- + Date : October 2016 + Copyright : (C) 2016 by Nyall Dawson + Email : nyall dot dawson at gmail dot com + *************************************************************************** + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + ***************************************************************************/ + +#ifndef QGSKDE_H +#define QGSKDE_H + +#include "qgsrectangle.h" +#include + +// GDAL includes +#include +#include +#include + +class QgsVectorLayer; +class QProgressDialog; +class QgsFeature; + + +/** + * \class QgsKernelDensityEstimation + * \ingroup analysis + * Performs Kernel Density Estimation ("heatmap") calculations on a vector layer. + * @note added in QGIS 3.0 + */ +class ANALYSIS_EXPORT QgsKernelDensityEstimation +{ + public: + + //! Kernel shape type + enum KernelShape + { + KernelQuartic = 0, //!< Quartic kernel + KernelTriangular, //!< Triangular kernel + KernelUniform, //!< Uniform (flat) kernel + KernelTriweight, //!< Triweight kernel + KernelEpanechnikov, //!< Epanechnikov kernel + }; + + //! Output values type + enum OutputValues + { + OutputRaw = 0, //!< Output the raw KDE values + OutputScaled, //!< Output mathematically correct scaled values + }; + + //! Result of operation + enum Result + { + Success, //!< Operation completed successfully + DriverError, //!< Could not open the driver for the specified format + InvalidParameters, //!< Input parameters were not valid + FileCreationError, //!< Error creating output file + RasterIoError, //!< Error writing to raster + }; + + //! KDE parameters + struct Parameters + { + //! Vector point layer + QgsVectorLayer* vectorLayer; + + //! Fixed radius, in map units + double radius; + + //! Field for radius, or empty if using a fixed radius + QString radiusField; + + //! Field name for weighting field, or empty if not using weights + QString weightField; + + //! Size of pixel in output file + double pixelSize; + + //! Kernel shape + KernelShape shape; + + //! Decay ratio (Triangular kernels only) + double decayRatio; + + //! Type of output value + OutputValues outputValues; + }; + + /** + * Constructor for QgsKernelDensityEstimation. Requires a Parameters object specifying the options to use + * to generate the surface. The output path and file format are also required. + */ + QgsKernelDensityEstimation( const Parameters& parameters, const QString& outputFile, const QString& outputFormat ); + + /** + * Runs the KDE calculation across the whole layer at once. Either call this method, or manually + * call run(), addFeature() and finalise() separately. + */ + Result run(); + + /** + * Prepares the output file for writing and setups up the surface calculation. This must be called + * before adding features via addFeature(). + * @see addFeature() + * @see finalise() + */ + Result prepare(); + + /** + * Adds a single feature to the KDE surface. prepare() must be called before adding features. + * @see prepare() + * @see finalise() + */ + Result addFeature( const QgsFeature& feature ); + + /** + * Finalises the output file. Must be called after adding all features via addFeature(). + * @see prepare() + * @see addFeature() + */ + Result finalise(); + + private: + + //! Calculate the value given to a point width a given distance for a specified kernel shape + double calculateKernelValue( const double distance, const int bandwidth, const KernelShape shape, const OutputValues outputType ) const; + //! Uniform kernel function + double uniformKernel( const double distance, const int bandwidth, const OutputValues outputType ) const; + //! Quartic kernel function + double quarticKernel( const double distance, const int bandwidth, const OutputValues outputType ) const; + //! Triweight kernel function + double triweightKernel( const double distance, const int bandwidth, const OutputValues outputType ) const; + //! Epanechnikov kernel function + double epanechnikovKernel( const double distance, const int bandwidth, const OutputValues outputType ) const; + //! Triangular kernel function + double triangularKernel( const double distance, const int bandwidth, const OutputValues outputType ) const; + + QgsRectangle calculateBounds() const; + + QgsVectorLayer* mInputLayer; + + QString mOutputFile; + QString mOutputFormat; + + int mRadiusField; + int mWeightField; + double mRadius; + double mPixelSize; + QgsRectangle mBounds; + + KernelShape mShape; + double mDecay; + OutputValues mOutputValues; + + int mBufferSize; + + GDALDatasetH mDatasetH; + GDALRasterBandH mRasterBandH; + + //! Creates a new raster layer and initialises it to the no data value + bool createEmptyLayer( GDALDriverH driver, const QgsRectangle& bounds, int rows, int columns ) const; + int radiusSizeInPixels( double radius ) const; +}; + + +#endif // QGSKDE_H