Skip to content

Commit

Permalink
merge pull request #104:
Browse files Browse the repository at this point in the history 
- redesigned the dialog, added a options for - weighted buffer - raster size
- moved the advanced stuff under advancedGroupBox
- finished the GUI code
- finally heatmap hits v0.2 Users can now
 - specify the rows and columns for the raster
 - specify the cell size of the raster
 - Choose a field as weight for each point
 - Choose a field to act as the buffer radius for each point
 - Buffer radius can be specified in real-world units like meters & mapunits
- indentation updated
  • Loading branch information
@jef-n
Arunmozhi authored and jef-n committed Mar 24, 2012
1 parent aa2d0b8 commit fc8b0c4
Show file tree
Hide file tree
Showing 5 changed files with 893 additions and 287 deletions.
341 changes: 193 additions & 148 deletions src/plugins/heatmap/heatmap.cpp
View file
Expand Up @@ -30,6 +30,9 @@
#include "qgsgeometry.h"
#include "qgsvectorlayer.h"
#include "qgsvectordataprovider.h"
#include "qgsdistancearea.h"
#include "qgscoordinatereferencesystem.h"
#include "qgslogger.h"

// Qt4 Related Includes
#include <QAction>
Expand All @@ -44,7 +47,7 @@
static const QString sName = QObject::tr( "Heatmap" );
static const QString sDescription = QObject::tr( "Creates a Heatmap raster for the input point vector" );
static const QString sCategory = QObject::tr( "Raster" );
static const QString sPluginVersion = QObject::tr( "Version 0.1" );
static const QString sPluginVersion = QObject::tr( "Version 0.2" );
static const QgisPlugin::PLUGINTYPE sPluginType = QgisPlugin::UI;
static const QString sPluginIcon = ":/heatmap/heatmap.png";

Expand Down Expand Up @@ -94,186 +97,228 @@ void Heatmap::help()
// not be enough
void Heatmap::run()
{
HeatmapGui *myPluginGui = new HeatmapGui( mQGisIface->mainWindow(), QgisGui::ModalDialogFlags );
myPluginGui->setAttribute( Qt::WA_DeleteOnClose );
HeatmapGui d( mQGisIface->mainWindow(), QgisGui::ModalDialogFlags );

// Connect the createRaster signal to createRaster Slot
connect( myPluginGui, SIGNAL( createRaster( QgsVectorLayer*, int, float, QString, QString ) ),
this, SLOT( createRaster( QgsVectorLayer*, int, float, QString, QString ) ) );

myPluginGui->show();
}

// Unload the plugin by cleaning up the GUI
void Heatmap::unload()
{
// remove the GUI
mQGisIface->removePluginRasterMenu( tr( "&Heatmap" ), mQActionPointer );
mQGisIface->removeRasterToolBarIcon( mQActionPointer );
delete mQActionPointer;
}

// The worker
void Heatmap::createRaster( QgsVectorLayer* theVectorLayer, int theBuffer, float theDecay, QString theOutputFilename, QString theOutputFormat )
{
// generic variables
int xSize, ySize;
double xResolution, yResolution;
double rasterX, rasterY;

// Getting the rasterdataset in place
GDALAllRegister();
if ( d.exec() == QDialog::Accepted )
{
// everything runs here

GDALDataset *emptyDataset;
GDALDriver *myDriver;
// Get the required data from the dialog
QgsRectangle myBBox = d.bbox();
int columns = d.columns();
int rows = d.rows();
float cellsize = d.cellSizeX(); // or d.cellSizeY(); both have the same value
float myDecay = d.decayRatio();

myDriver = GetGDALDriverManager()->GetDriverByName( theOutputFormat.toUtf8() );
if ( myDriver == NULL )
{
QMessageBox::information( 0, tr( "GDAL driver error" ), tr( "Cannot open the driver for the specified format" ) );
return;
}
// Getting the rasterdataset in place
GDALAllRegister();

// bounding box info
QgsRectangle myBBox = theVectorLayer->extent();
// fixing a base width of 500 px/cells
xSize = 500;
xResolution = myBBox.width() / xSize;
yResolution = xResolution;
ySize = myBBox.height() / yResolution;
// add extra extend to cover the corner points' heat region
xSize = xSize + ( theBuffer * 2 ) + 10 ;
ySize = ySize + ( theBuffer * 2 ) + 10 ;
// Define the new lat,lon for the buffered raster area
rasterX = myBBox.xMinimum() - ( theBuffer + 5 ) * xResolution;
rasterY = myBBox.yMinimum() - ( theBuffer + 5 ) * yResolution;

double geoTransform[6] = { rasterX, xResolution, 0, rasterY, 0, yResolution };

emptyDataset = myDriver->Create( theOutputFilename.toUtf8(), xSize, ySize, 1, GDT_Float32, NULL );

emptyDataset->SetGeoTransform( geoTransform );

GDALRasterBand *poBand;
poBand = emptyDataset->GetRasterBand( 1 );
poBand->SetNoDataValue( NO_DATA );

float* line = ( float * ) CPLMalloc( sizeof( float ) * xSize );
for ( int i = 0; i < xSize; i++ )
line[i] = NO_DATA;
// Write the empty raster
for ( int i = 0; i < ySize ; i++ )
{
poBand->RasterIO( GF_Write, 0, 0, xSize, 1, line, xSize, 1, GDT_Float32, 0, 0 );
}
GDALDataset *emptyDataset;
GDALDriver *myDriver;

CPLFree( line );
//close the dataset
GDALClose(( GDALDatasetH ) emptyDataset );
myDriver = GetGDALDriverManager()->GetDriverByName( d.outputFormat().toUtf8() );
if ( myDriver == NULL )
{
QMessageBox::information( 0, tr( "GDAL driver error" ), tr( "Cannot open the driver for the specified format" ) );
return;
}

// open the raster in GA_Update mode
GDALDataset *heatmapDS;
heatmapDS = ( GDALDataset * ) GDALOpen( theOutputFilename.toUtf8(), GA_Update );
if ( !heatmapDS )
{
QMessageBox::information( 0, tr( "Raster update error" ), tr( "Could not open the created raster for updating. The heatmap was not generated." ) );
return;
}
poBand = heatmapDS->GetRasterBand( 1 );
// Get the data buffer ready
int blockSize = 2 * theBuffer + 1; // block SIDE would have been more appropriate
// Open the vector features
QgsVectorDataProvider* myVectorProvider = theVectorLayer->dataProvider();
if ( !myVectorProvider )
{
QMessageBox::information( 0, tr( "Point layer error" ), tr( "Could not identify the vector data provider." ) );
return;
}
QgsAttributeList dummyList;
myVectorProvider->select( dummyList );
double geoTransform[6] = { myBBox.xMinimum(), cellsize, 0, myBBox.yMinimum(), 0, cellsize };
emptyDataset = myDriver->Create( d.outputFilename().toUtf8(), columns, rows, 1, GDT_Float32, NULL );
emptyDataset->SetGeoTransform( geoTransform );

int totalFeatures = myVectorProvider->featureCount();
int counter = 0;
GDALRasterBand *poBand;
poBand = emptyDataset->GetRasterBand( 1 );
poBand->SetNoDataValue( NO_DATA );

QProgressDialog p( "Creating Heatmap ... ", "Abort", 0, totalFeatures );
p.setWindowModality( Qt::WindowModal );
float* line = ( 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++ )
{
poBand->RasterIO( GF_Write, 0, 0, columns, 1, line, columns, 1, GDT_Float32, 0, 0 );
}

QgsFeature myFeature;
CPLFree( line );
//close the dataset
GDALClose(( GDALDatasetH ) emptyDataset );

while ( myVectorProvider->nextFeature( myFeature ) )
{
counter++;
p.setValue( counter );
if ( p.wasCanceled() )
// open the raster in GA_Update mode
GDALDataset *heatmapDS;
heatmapDS = ( GDALDataset * ) GDALOpen( d.outputFilename().toUtf8(), GA_Update );
if ( !heatmapDS )
{
QMessageBox::information( 0, tr( "Raster update error" ), tr( "Could not open the created raster for updating. The heatmap was not generated." ) );
return;
}
poBand = heatmapDS->GetRasterBand( 1 );
// Start working on the input vector
QgsVectorLayer* inputLayer = d.inputVectorLayer();
QgsVectorDataProvider* myVectorProvider = inputLayer->dataProvider();
if ( !myVectorProvider )
{
QMessageBox::information( 0, tr( "Heatmap generation aborted" ), tr( "QGIS will now load the partially-computed raster." ) );
break;
QMessageBox::information( 0, tr( "Point layer error" ), tr( "Could not identify the vector data provider." ) );
return;
}

QgsGeometry* myPointGeometry;
myPointGeometry = myFeature.geometry();
// convert the geometry to point
QgsPoint myPoint;
myPoint = myPointGeometry->asPoint();
// avoiding any empty points or out of extent points
if (( myPoint.x() < rasterX ) || ( myPoint.y() < rasterY ) )
QgsAttributeList myAttrList;
int rField = 0;
int wField = 0;
if ( d.variableRadius() )
{
rField = d.radiusField();
myAttrList.append( rField );
QgsDebugMsg( tr( "Radius Field index received: %1" ).arg( rField ) );
}
if ( d.weighted() )
{
continue;
wField = d.weightField();
myAttrList.append( wField );
}
// calculate the pixel position
unsigned int xPosition, yPosition;
xPosition = (( myPoint.x() - rasterX ) / xResolution ) - theBuffer;
yPosition = (( myPoint.y() - rasterY ) / yResolution ) - theBuffer;
// This might have attributes or mightnot have attibutes at all
// based on the variableRadius() and weighted()
myVectorProvider->select( myAttrList );
int totalFeatures = myVectorProvider->featureCount();
int counter = 0;

// get the data
float *dataBuffer = ( float * ) CPLMalloc( sizeof( float ) * blockSize * blockSize );
poBand->RasterIO( GF_Read, xPosition, yPosition, blockSize, blockSize, dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );
QProgressDialog p( "Creating Heatmap ... ", "Abort", 0, totalFeatures );
p.setWindowModality( Qt::WindowModal );

for ( int xp = 0; xp <= theBuffer; xp++ )
QgsFeature myFeature;

while ( myVectorProvider->nextFeature( myFeature ) )
{
for ( int yp = 0; yp <= theBuffer; yp++ )
counter++;
p.setValue( counter );
if ( p.wasCanceled() )
{
float distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
float pixelValue = 1 - (( 1 - theDecay ) * distance / theBuffer );
QMessageBox::information( 0, tr( "Heatmap generation aborted" ), tr( "QGIS will now load the partially-computed raster." ) );
break;
}

// clearing anamolies along the axes
if ( xp == 0 && yp == 0 )
{
pixelValue /= 4;
}
else if ( xp == 0 || yp == 0 )
{
pixelValue /= 2;
}
QgsGeometry* myPointGeometry;
myPointGeometry = myFeature.geometry();
// convert the geometry to point
QgsPoint myPoint;
myPoint = myPointGeometry->asPoint();
// avoiding any empty points or out of extent points
if (( myPoint.x() < myBBox.xMinimum() ) || ( myPoint.y() < myBBox.yMinimum() ) )
{
continue;
}
float radius;
if ( d.variableRadius() )
{
QgsAttributeMap myAttrMap = myFeature.attributeMap();
radius = myAttrMap.value( rField ).toFloat();
}
else
{
radius = d.radius();
}
//convert the radius to map units if it is in meters
if ( d.radiusUnit() == HeatmapGui::Meters )
{
radius = mapUnitsOf( radius, inputLayer->crs() );
}
// convert radius in map units to pixel count
int myBuffer = radius / cellsize;
if ( radius - ( cellsize * myBuffer ) > 0.5 )
{
++myBuffer;
}
int blockSize = 2 * myBuffer + 1; //Block SIDE would be more appropriate
// calculate the pixel position
unsigned int xPosition, yPosition;
xPosition = (( myPoint.x() - myBBox.xMinimum() ) / cellsize ) - myBuffer;
yPosition = (( myPoint.y() - myBBox.yMinimum() ) / cellsize ) - myBuffer;

// get the data
float *dataBuffer = ( float * ) CPLMalloc( sizeof( float ) * blockSize * blockSize );
poBand->RasterIO( GF_Read, xPosition, yPosition, blockSize, blockSize,
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );

float weight = 1.0;
if ( d.weighted() )
{
QgsAttributeMap myAttrMap = myFeature.attributeMap();
weight = myAttrMap.value( wField ).toFloat();
}

if ( distance <= theBuffer )
for ( int xp = 0; xp <= myBuffer; xp++ )
{
for ( int yp = 0; yp <= myBuffer; yp++ )
{
int pos[4];
pos[0] = ( theBuffer + xp ) * blockSize + ( theBuffer + yp );
pos[1] = ( theBuffer + xp ) * blockSize + ( theBuffer - yp );
pos[2] = ( theBuffer - xp ) * blockSize + ( theBuffer + yp );
pos[3] = ( theBuffer - xp ) * blockSize + ( theBuffer - yp );
for ( int p = 0; p < 4; p++ )
float distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
float pixelValue = weight * ( 1 - (( 1 - myDecay ) * distance / myBuffer ) );

// clearing anamolies along the axes
if ( xp == 0 && yp == 0 )
{
pixelValue /= 4;
}
else if ( xp == 0 || yp == 0 )
{
if ( dataBuffer[ pos[p] ] == NO_DATA )
pixelValue /= 2;
}

if ( distance <= myBuffer )
{
int pos[4];
pos[0] = ( myBuffer + xp ) * blockSize + ( myBuffer + yp );
pos[1] = ( myBuffer + xp ) * blockSize + ( myBuffer - yp );
pos[2] = ( myBuffer - xp ) * blockSize + ( myBuffer + yp );
pos[3] = ( myBuffer - xp ) * blockSize + ( myBuffer - yp );
for ( int p = 0; p < 4; p++ )
{
dataBuffer[ pos[p] ] = 0;
if ( dataBuffer[ pos[p] ] == NO_DATA )
{
dataBuffer[ pos[p] ] = 0;
}
dataBuffer[ pos[p] ] += pixelValue;
}
dataBuffer[ pos[p] ] += pixelValue;
}
}
}

poBand->RasterIO( GF_Write, xPosition, yPosition, blockSize, blockSize,
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );
CPLFree( dataBuffer );
}
//Finally close the dataset
GDALClose(( GDALDatasetH ) heatmapDS );

poBand->RasterIO( GF_Write, xPosition, yPosition, blockSize, blockSize, dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );
CPLFree( dataBuffer );
// Open the file in QGIS window
mQGisIface->addRasterLayer( d.outputFilename(), QFileInfo( d.outputFilename() ).baseName() );
}
}

//Finally close the dataset
GDALClose(( GDALDatasetH ) heatmapDS );
/*
*
* Local functions
*
*/
float Heatmap::mapUnitsOf( float meters, QgsCoordinateReferenceSystem layerCrs )
{
// Worker to transform metres input to mapunits
QgsDistanceArea da;
da.setSourceCrs( layerCrs.srsid() );
da.setEllipsoid( layerCrs.ellipsoidAcronym() );
if ( da.geographic() )
{
da.setProjectionsEnabled( true );
}
return meters / da.measureLine( QgsPoint( 0.0, 0.0 ), QgsPoint( 0.0, 1.0 ) );
}

// Open the file in QGIS window
mQGisIface->addRasterLayer( theOutputFilename, QFileInfo( theOutputFilename ).baseName() );
// Unload the plugin by cleaning up the GUI
void Heatmap::unload()
{
// remove the GUI
mQGisIface->removePluginRasterMenu( tr( "&Heatmap" ), mQActionPointer );
mQGisIface->removeRasterToolBarIcon( mQActionPointer );
delete mQActionPointer;
}

//////////////////////////////////////////////////////////////////////////
Expand Down

0 comments on commit fc8b0c4

Please sign in to comment.