24 : mSource( parameters.source )
25 , mOutputFile( outputFile )
26 , mOutputFormat( outputFormat )
29 , mRadius( parameters.radius )
30 , mPixelSize( parameters.pixelSize )
31 , mShape( parameters.shape )
32 , mDecay( parameters.decayRatio )
33 , mOutputValues( parameters.outputValues )
35 , mDatasetH( nullptr )
36 , mRasterBandH( nullptr )
52 if ( mRadiusField >= 0 )
53 requiredAttributes << mRadiusField;
55 if ( mWeightField >= 0 )
56 requiredAttributes << mWeightField;
73 GDALDriverH driver = GDALGetDriverByName( mOutputFormat.toUtf8() );
82 mBounds = calculateBounds();
86 const int rows = std::max( std::ceil( mBounds.
height() / mPixelSize ) + 1, 1.0 );
87 const int cols = std::max( std::ceil( mBounds.
width() / mPixelSize ) + 1, 1.0 );
89 if ( !createEmptyLayer( driver, mBounds, rows, cols ) )
93 mDatasetH.reset( GDALOpen( mOutputFile.toUtf8().constData(), GA_Update ) );
96 mRasterBandH = GDALGetRasterBand( mDatasetH.get(), 1 );
101 if ( mRadiusField < 0 )
102 mBufferSize = radiusSizeInPixels( mRadius );
110 if ( featureGeometry.
isNull() )
132 double radius = mRadius;
133 int buffer = mBufferSize;
134 if ( mRadiusField >= 0 )
136 radius = feature.
attribute( mRadiusField ).toDouble();
137 buffer = radiusSizeInPixels( radius );
139 const int blockSize = 2 * buffer + 1;
143 if ( mWeightField >= 0 )
145 weight = feature.
attribute( mWeightField ).toDouble();
151 for ( QgsMultiPointXY::const_iterator pointIt = multiPoints.constBegin(); pointIt != multiPoints.constEnd(); ++pointIt )
154 if ( !mBounds.
contains( *pointIt ) )
160 const unsigned int xPosition = ( ( ( *pointIt ).x() - mBounds.
xMinimum() ) / mPixelSize ) - buffer;
161 const unsigned int yPosition = ( ( ( *pointIt ).y() - mBounds.
yMinimum() ) / mPixelSize ) - buffer;
162 const unsigned int yPositionIO = ( ( mBounds.
yMaximum() - ( *pointIt ).y() ) / mPixelSize ) - buffer;
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 )
173 for (
int xp = 0; xp < blockSize; xp++ )
175 for (
int yp = 0; yp < blockSize; yp++ )
177 const double pixelCentroidX = ( xPosition + xp + 0.5 ) * mPixelSize + mBounds.
xMinimum();
178 const double pixelCentroidY = ( yPosition + yp + 0.5 ) * mPixelSize + mBounds.
yMinimum();
180 const double distance = std::sqrt( std::pow( pixelCentroidX - ( *pointIt ).x(), 2.0 ) + std::pow( pixelCentroidY - ( *pointIt ).y(), 2.0 ) );
183 if ( distance > radius )
188 const double pixelValue = weight * calculateKernelValue( distance, radius, mShape, mOutputValues );
189 const int pos = xp + blockSize * yp;
190 if ( dataBuffer[ pos ] ==
NO_DATA )
192 dataBuffer[ pos ] = 0;
194 dataBuffer[ pos ] += pixelValue;
197 if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPositionIO, blockSize, blockSize,
198 dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
202 CPLFree( dataBuffer );
211 mRasterBandH =
nullptr;
215int QgsKernelDensityEstimation::radiusSizeInPixels(
double radius )
const
217 int buffer = radius / mPixelSize;
218 if ( radius - ( mPixelSize * buffer ) > 0.5 )
225bool QgsKernelDensityEstimation::createEmptyLayer( GDALDriverH driver,
const QgsRectangle &bounds,
int rows,
int columns )
const
227 double geoTransform[6] = { bounds.
xMinimum(), mPixelSize, 0, bounds.
yMaximum(), 0, -mPixelSize };
228 const gdal::dataset_unique_ptr emptyDataset( GDALCreate( driver, mOutputFile.toUtf8(), columns, rows, 1, GDT_Float32,
nullptr ) );
232 if ( GDALSetGeoTransform( emptyDataset.get(), geoTransform ) != CE_None )
236 if ( GDALSetProjection( emptyDataset.get(), mSource->
sourceCrs().
toWkt().toLocal8Bit().data() ) != CE_None )
239 GDALRasterBandH poBand = GDALGetRasterBand( emptyDataset.get(), 1 );
243 if ( GDALSetRasterNoDataValue( poBand,
NO_DATA ) != CE_None )
246 float *line =
static_cast< float *
>( CPLMalloc(
sizeof(
float ) * columns ) );
247 for (
int i = 0; i < columns; i++ )
252 for (
int i = 0; i < rows ; i++ )
254 if ( GDALRasterIO( poBand, GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 ) != CE_None )
269 return triangularKernel( distance, bandwidth, outputType );
272 return uniformKernel( distance, bandwidth, outputType );
275 return quarticKernel( distance, bandwidth, outputType );
278 return triweightKernel( distance, bandwidth, outputType );
281 return epanechnikovKernel( distance, bandwidth, outputType );
297 switch ( outputType )
302 const double k = 2. / ( M_PI * bandwidth );
305 return k * ( 0.5 / bandwidth );
315 switch ( outputType )
320 const double k = 116. / ( 5. * M_PI * std::pow( bandwidth, 2 ) );
323 return k * ( 15. / 16. ) * std::pow( 1. - std::pow( distance / bandwidth, 2 ), 2 );
326 return std::pow( 1. - std::pow( distance / bandwidth, 2 ), 2 );
333 switch ( outputType )
338 const double k = 128. / ( 35. * M_PI * std::pow( bandwidth, 2 ) );
341 return k * ( 35. / 32. ) * std::pow( 1. - std::pow( distance / bandwidth, 2 ), 3 );
344 return std::pow( 1. - std::pow( distance / bandwidth, 2 ), 3 );
351 switch ( outputType )
356 const double k = 8. / ( 3. * M_PI * std::pow( bandwidth, 2 ) );
359 return k * ( 3. / 4. ) * ( 1. - std::pow( distance / bandwidth, 2 ) );
362 return ( 1. - std::pow( distance / bandwidth, 2 ) );
370 switch ( outputType )
379 const double k = 3. / ( ( 1. + 2. * mDecay ) * M_PI * std::pow( bandwidth, 2 ) );
382 return k * ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
387 return ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
391 return ( 1. - ( 1. - mDecay ) * ( distance / bandwidth ) );
396QgsRectangle QgsKernelDensityEstimation::calculateBounds()
const
404 if ( mRadiusField >= 0 )
407 radius = mSource->
maximumValue( mRadiusField ).toDouble();
QString toWkt(Qgis::CrsWktVariant variant=Qgis::CrsWktVariant::Wkt1Gdal, bool multiline=false, int indentationWidth=4) const
Returns a WKT representation of this CRS.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
virtual QgsFields fields() const =0
Returns the fields associated with features in the source.
virtual QgsCoordinateReferenceSystem sourceCrs() const =0
Returns the coordinate reference system for features in the source.
virtual QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const =0
Returns an iterator for the features in the source.
virtual QVariant maximumValue(int fieldIndex) const
Returns the maximum value for an attribute column or an invalid variant in case of error.
virtual QgsRectangle sourceExtent() const
Returns the extent of all geometries from the source.
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Q_INVOKABLE QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Q_INVOKABLE int lookupField(const QString &fieldName) const
Looks up field's index from the field name.
A geometry is the spatial representation of a feature.
QgsMultiPointXY asMultiPoint() const
Returns the contents of the geometry as a multi-point.
QgsPointXY asPoint() const
Returns the contents of the geometry as a 2-dimensional point.
bool isMultipart() const
Returns true if WKB of the geometry is of WKBMulti* type.
QgsKernelDensityEstimation(const Parameters ¶meters, const QString &outputFile, const QString &outputFormat)
Constructor for QgsKernelDensityEstimation.
Result run()
Runs the KDE calculation across the whole layer at once.
Result addFeature(const QgsFeature &feature)
Adds a single feature to the KDE surface.
Result finalise()
Finalises the output file.
Result
Result of operation.
@ DriverError
Could not open the driver for the specified format.
@ FileCreationError
Error creating output file.
@ RasterIoError
Error writing to raster.
@ Success
Operation completed successfully.
@ InvalidParameters
Input parameters were not valid.
OutputValues
Output values type.
@ OutputRaw
Output the raw KDE values.
@ OutputScaled
Output mathematically correct scaled values.
KernelShape
Kernel shape type.
@ KernelTriweight
Triweight kernel.
@ KernelUniform
Uniform (flat) kernel.
@ KernelEpanechnikov
Epanechnikov kernel.
@ KernelTriangular
Triangular kernel.
@ KernelQuartic
Quartic kernel.
Result prepare()
Prepares the output file for writing and setups up the surface calculation.
A class to represent a 2D point.
A rectangle specified with double values.
bool contains(const QgsRectangle &rect) const
Returns true when rectangle contains other rectangle.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
void setYMinimum(double y)
Set the minimum y value.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
void setXMinimum(double x)
Set the minimum x value.
double width() const
Returns the width of the rectangle.
double xMaximum() const
Returns the x maximum value (right side of rectangle).
bool isNull() const
Test if the rectangle is null (holding no spatial information).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
double height() const
Returns the height of the rectangle.
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
QList< int > QgsAttributeList
QVector< QgsPointXY > QgsMultiPointXY
A collection of QgsPoints that share a common collection of attributes.
QString radiusField
Field for radius, or empty if using a fixed radius.
QString weightField
Field name for weighting field, or empty if not using weights.