Heatmap v0.2

src/plugins/heatmap/heatmap.cpp

@@ -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>
@@ -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";
 
@@ -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 recieved: %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;
 }
 
 //////////////////////////////////////////////////////////////////////////