Skip to content
Permalink
Browse files

Fix heatmap failing for multipoint tables, fix incorrect calculation

of heatmap input layer.
  • Loading branch information
nyalldawson committed Sep 24, 2014
1 parent ace704a commit 19125b82e04046a73b3b95deaf27cd1ad8bbfcc1
Showing with 79 additions and 51 deletions.
  1. +79 −51 src/plugins/heatmap/heatmap.cpp
@@ -112,9 +112,8 @@ void Heatmap::run()
{
HeatmapGui d( mQGisIface->mainWindow(), QgisGui::ModalDialogFlags, &mSessionSettings );

// Start working on the input vector
QgsVectorLayer* inputLayer = d.inputVectorLayer();
if ( !inputLayer )
//check that dialog found a suitable vector layer
if ( !d.inputVectorLayer() )
{
mQGisIface->messageBar()->pushMessage( tr( "Layer not found" ), tr( "The heatmap plugin requires at least one point vector layer" ), QgsMessageBar::INFO, mQGisIface->messageTimeout() );
return;
@@ -123,6 +122,7 @@ void Heatmap::run()
if ( d.exec() == QDialog::Accepted )
{
// everything runs here
QgsVectorLayer* inputLayer = d.inputVectorLayer();

// Get the required data from the dialog
QgsRectangle myBBox = d.bbox();
@@ -132,6 +132,9 @@ void Heatmap::run()
mDecay = d.decayRatio();
int kernelShape = d.kernelShape();

//is input layer multipoint?
bool isMultiPoint = inputLayer->wkbType() == QGis::WKBMultiPoint || inputLayer->wkbType() == QGis::WKBMultiPoint25D;

// Getting the rasterdataset in place
GDALAllRegister();

@@ -235,18 +238,30 @@ void Heatmap::run()
break;
}

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() )
|| ( myPoint.x() > myBBox.xMaximum() ) || ( myPoint.y() > myBBox.yMaximum() ) )
QgsGeometry* featureGeometry = myFeature.geometry();
if ( !featureGeometry )
{
continue;
}

// convert the geometry to multipoint
QgsMultiPoint multiPoints;
if ( !isMultiPoint )
{
QgsPoint myPoint = featureGeometry->asPoint();
// avoiding any empty points or out of extent points
if (( myPoint.x() < myBBox.xMinimum() ) || ( myPoint.y() < myBBox.yMinimum() )
|| ( myPoint.x() > myBBox.xMaximum() ) || ( myPoint.y() > myBBox.yMaximum() ) )
{
continue;
}
multiPoints << myPoint;
}
else
{
multiPoints = featureGeometry->asMultiPoint();
}

// If radius is variable then fetch it and calculate new pixel buffer size
if ( d.variableRadius() )
{
@@ -255,66 +270,79 @@ void Heatmap::run()
}

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 );

double weight = 1.0;
if ( d.weighted() )
{
weight = myFeature.attribute( wField ).toDouble();
}

for ( int xp = 0; xp <= myBuffer; xp++ )
//loop through all points in multipoint
for ( QgsMultiPoint::const_iterator pointIt = multiPoints.constBegin(); pointIt != multiPoints.constEnd(); ++pointIt )
{
for ( int yp = 0; yp <= myBuffer; yp++ )
// avoiding any empty points or out of extent points
if ((( *pointIt ).x() < myBBox.xMinimum() ) || (( *pointIt ).y() < myBBox.yMinimum() )
|| (( *pointIt ).x() > myBBox.xMaximum() ) || (( *pointIt ).y() > myBBox.yMaximum() ) )
{
double distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
continue;
}

// is pixel outside search bandwidth of feature?
if ( distance > myBuffer )
{
continue;
}
// calculate the pixel position
unsigned int xPosition, yPosition;
xPosition = ((( *pointIt ).x() - myBBox.xMinimum() ) / cellsize ) - myBuffer;
yPosition = ((( *pointIt ).y() - myBBox.yMinimum() ) / cellsize ) - myBuffer;

double pixelValue = weight * calculateKernelValue( distance, myBuffer, kernelShape );
// 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 );

// clearing anamolies along the axes
if ( xp == 0 && yp == 0 )
{
pixelValue /= 4;
}
else if ( xp == 0 || yp == 0 )
for ( int xp = 0; xp <= myBuffer; xp++ )
{
for ( int yp = 0; yp <= myBuffer; yp++ )
{
pixelValue /= 2;
}
double distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );

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++ )
{
if ( dataBuffer[ pos[p] ] == NO_DATA )
// is pixel outside search bandwidth of feature?
if ( distance > myBuffer )
{
continue;
}

double pixelValue = weight * calculateKernelValue( distance, myBuffer, kernelShape );

// 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] = ( 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 );
}

poBand->RasterIO( GF_Write, xPosition, yPosition, blockSize, blockSize,
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );
CPLFree( dataBuffer );
}


// Finally close the dataset
GDALClose(( GDALDatasetH ) heatmapDS );

0 comments on commit 19125b8

Please sign in to comment.
You can’t perform that action at this time.