Skip to content

Commit

Permalink
use QGIS API instead of GDAL API to create heatmaps
Browse files Browse the repository at this point in the history
  • Loading branch information
@alexbruy @nyalldawson
alexbruy authored and nyalldawson committed Nov 25, 2019
1 parent 42d6c39 commit 8255932
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 80 deletions.
109 changes: 38 additions & 71 deletions src/analysis/raster/qgskde.cpp
View file
Expand Up @@ -13,10 +13,13 @@
* *
***************************************************************************/

#include <QByteArray>

#include "qgskde.h"
#include "qgsfeaturesource.h"
#include "qgsfeatureiterator.h"
#include "qgsgeometry.h"
#include "qgsrasterfilewriter.h"

#define NO_DATA -9999

Expand All @@ -32,8 +35,6 @@ QgsKernelDensityEstimation::QgsKernelDensityEstimation( const QgsKernelDensityEs
, mDecay( parameters.decayRatio )
, mOutputValues( parameters.outputValues )
, mBufferSize( -1 )
, mDatasetH( nullptr )
, mRasterBandH( nullptr )
{
if ( !parameters.radiusField.isEmpty() )
mRadiusField = mSource->fields().lookupField( parameters.radiusField );
Expand Down Expand Up @@ -68,14 +69,6 @@ QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::run()

QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::prepare()
{
GDALAllRegister();

GDALDriverH driver = GDALGetDriverByName( mOutputFormat.toUtf8() );
if ( !driver )
{
return DriverError;
}

if ( !mSource )
return InvalidParameters;

Expand All @@ -86,16 +79,30 @@ QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::prepare()
int rows = std::max( std::ceil( mBounds.height() / mPixelSize ) + 1, 1.0 );
int cols = std::max( std::ceil( mBounds.width() / mPixelSize ) + 1, 1.0 );

if ( !createEmptyLayer( driver, mBounds, rows, cols ) )
// create empty raster and fill it with nodata values
QgsRasterFileWriter writer( mOutputFile );
writer.setOutputProviderKey( QStringLiteral( "gdal" ) );
writer.setOutputFormat( mOutputFormat );
mProvider = writer.createOneBandRaster( Qgis::Float32, cols, rows, mBounds, mSource->sourceCrs() );
if ( !mProvider )
{
return FileCreationError;
}

// open the raster in GA_Update mode
mDatasetH.reset( GDALOpen( mOutputFile.toUtf8().constData(), GA_Update ) );
if ( !mDatasetH )
return FileCreationError;
mRasterBandH = GDALGetRasterBand( mDatasetH.get(), 1 );
if ( !mRasterBandH )
if ( !mProvider->setNoDataValue( 1, NO_DATA ) )
{
return FileCreationError;
}

int dataTypeSize = QgsRasterBlock::typeSize( Qgis::Float32 );
std::vector<float> line( cols );
std::fill( line.begin(), line.end(), NO_DATA );
QgsRasterBlock block( Qgis::Float32, cols, 1 );
block.setData( QByteArray::fromRawData( ( char * )&line[0], dataTypeSize * cols ) );
for ( int i = 0; i < rows ; i++ )
{
mProvider->writeBlock( &block, 1, 0, i );
}

mBufferSize = -1;
if ( mRadiusField < 0 )
Expand Down Expand Up @@ -162,13 +169,13 @@ QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::addFeature( const
unsigned int yPositionIO = ( ( mBounds.yMaximum() - ( *pointIt ).y() ) / mPixelSize ) - buffer;


// calculate buffer around pixel
QgsRectangle extent( ( *pointIt ).x() - radius, ( *pointIt ).y() - radius, ( *pointIt ).x() + radius, ( *pointIt ).y() + radius );

// get the data
float *dataBuffer = ( float * ) CPLMalloc( sizeof( float ) * blockSize * blockSize );
if ( GDALRasterIO( mRasterBandH, GF_Read, xPosition, yPositionIO, blockSize, blockSize,
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
{
result = RasterIoError;
}
QgsRasterBlock *block = mProvider->block( 1, extent, blockSize, blockSize );
QByteArray blockData = block->data();
float *dataBuffer = ( float * ) blockData.data();

for ( int xp = 0; xp < blockSize; xp++ )
{
Expand All @@ -194,21 +201,24 @@ QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::addFeature( const
dataBuffer[ pos ] += pixelValue;
}
}
if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPositionIO, blockSize, blockSize,
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )

block->setData( blockData );
if ( !mProvider->writeBlock( block, 1, xPosition, yPositionIO ) )
{
result = RasterIoError;
}
CPLFree( dataBuffer );

delete block;
}

return result;
}

QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::finalise()
{
mDatasetH.reset();
mRasterBandH = nullptr;
mProvider->setEditable( false );
mProvider = nullptr;

return Success;
}

Expand All @@ -222,45 +232,6 @@ int QgsKernelDensityEstimation::radiusSizeInPixels( double radius ) const
return buffer;
}

bool QgsKernelDensityEstimation::createEmptyLayer( GDALDriverH driver, const QgsRectangle &bounds, int rows, int columns ) const
{
double geoTransform[6] = { bounds.xMinimum(), mPixelSize, 0, bounds.yMaximum(), 0, -mPixelSize };
gdal::dataset_unique_ptr emptyDataset( GDALCreate( driver, mOutputFile.toUtf8(), columns, rows, 1, GDT_Float32, nullptr ) );
if ( !emptyDataset )
return false;

if ( GDALSetGeoTransform( emptyDataset.get(), geoTransform ) != CE_None )
return false;

// Set the projection on the raster destination to match the input layer
if ( GDALSetProjection( emptyDataset.get(), mSource->sourceCrs().toWkt().toLocal8Bit().data() ) != CE_None )
return false;

GDALRasterBandH poBand = GDALGetRasterBand( emptyDataset.get(), 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 );
return true;
}

double QgsKernelDensityEstimation::calculateKernelValue( const double distance, const double bandwidth, const QgsKernelDensityEstimation::KernelShape shape, const QgsKernelDensityEstimation::OutputValues outputType ) const
{
switch ( shape )
Expand Down Expand Up @@ -417,7 +388,3 @@ QgsRectangle QgsKernelDensityEstimation::calculateBounds() const
bbox.setYMaximum( bbox.yMaximum() + radius );
return bbox;
}




12 changes: 3 additions & 9 deletions src/analysis/raster/qgskde.h
View file
Expand Up @@ -16,15 +16,12 @@
#ifndef QGSKDE_H
#define QGSKDE_H

#include "qgsrectangle.h"
#include "qgsogrutils.h"
#include <QString>

// GDAL includes
#include <gdal.h>
#include <cpl_string.h>
#include <cpl_conv.h>
#include "qgis_analysis.h"
#include "qgsrectangle.h"
#include "qgsrasterdataprovider.h"

class QgsFeatureSource;
class QgsFeature;
Expand Down Expand Up @@ -168,11 +165,8 @@ class ANALYSIS_EXPORT QgsKernelDensityEstimation

int mBufferSize;

gdal::dataset_unique_ptr mDatasetH;
GDALRasterBandH mRasterBandH;
QgsRasterDataProvider *mProvider = nullptr;

//! Creates a new raster layer and initializes it to the no data value
bool createEmptyLayer( GDALDriverH driver, const QgsRectangle &bounds, int rows, int columns ) const;
int radiusSizeInPixels( double radius ) const;

#ifdef SIP_RUN
Expand Down

0 comments on commit 8255932

Please sign in to comment.