Revert usage of QGIS native raster API in KDE as it causes issues

qgis · Apr 13, 2020 · 6d9fc8a · 6d9fc8a
1 parent 2fd0ff5
commit 6d9fc8a
Show file tree

Hide file tree

Showing 2 changed files with 80 additions and 37 deletions.
diff --git a/src/analysis/raster/qgskde.cpp b/src/analysis/raster/qgskde.cpp
@@ -13,13 +13,10 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <QByteArray>
-
 #include "qgskde.h"
 #include "qgsfeaturesource.h"
 #include "qgsfeatureiterator.h"
 #include "qgsgeometry.h"
-#include "qgsrasterfilewriter.h"
 
 #define NO_DATA -9999
 
@@ -35,6 +32,8 @@ 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 );
@@ -69,6 +68,14 @@ QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::run()
 
 QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::prepare()
 {
+  GDALAllRegister();
+
+  GDALDriverH driver = GDALGetDriverByName( mOutputFormat.toUtf8() );
+  if ( !driver )
+  {
+    return DriverError;
+  }
+
   if ( !mSource )
     return InvalidParameters;
 
@@ -79,28 +86,16 @@ 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 );
 
-  // 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 )
-  {
+  if ( !createEmptyLayer( driver, mBounds, rows, cols ) )
     return FileCreationError;
-  }
 
-  if ( !mProvider->setNoDataValue( 1, NO_DATA ) )
-  {
+  // 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 )
     return FileCreationError;
-  }
-
-  QgsRasterBlock block( Qgis::Float32, cols, 1 );
-  block.setNoDataValue( NO_DATA );
-  block.setIsNoData();
-  for ( int i = 0; i < rows ; i++ )
-  {
-    mProvider->writeBlock( &block, 1, 0, i );
-  }
 
   mBufferSize = -1;
   if ( mRadiusField < 0 )
@@ -167,13 +162,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
-    std::unique_ptr< QgsRasterBlock > block( mProvider->block( 1, extent, blockSize, blockSize ) );
-    QByteArray blockData = block->data();
-    float *dataBuffer = ( float * ) blockData.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;
+    }
 
     for ( int xp = 0; xp < blockSize; xp++ )
     {
@@ -199,22 +194,21 @@ QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::addFeature( const
         dataBuffer[ pos ] += pixelValue;
       }
     }
-
-    block->setData( blockData );
-    if ( !mProvider->writeBlock( block.get(), 1, xPosition, yPositionIO ) )
+    if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPositionIO, blockSize, blockSize,
+                       dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
     {
       result = RasterIoError;
     }
+    CPLFree( dataBuffer );
   }
 
   return result;
 }
 
 QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::finalise()
 {
-  mProvider->setEditable( false );
-  mProvider = nullptr;
-
+  mDatasetH.reset();
+  mRasterBandH = nullptr;
   return Success;
 }
 
@@ -228,6 +222,45 @@ 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 )
@@ -384,3 +417,7 @@ QgsRectangle QgsKernelDensityEstimation::calculateBounds() const
   bbox.setYMaximum( bbox.yMaximum() + radius );
   return bbox;
 }
+
+
+
+
diff --git a/src/analysis/raster/qgskde.h b/src/analysis/raster/qgskde.h
@@ -16,12 +16,15 @@
 #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;
@@ -165,8 +168,11 @@ class ANALYSIS_EXPORT QgsKernelDensityEstimation
 
     int mBufferSize;
 
-    QgsRasterDataProvider *mProvider = nullptr;
+    gdal::dataset_unique_ptr mDatasetH;
+    GDALRasterBandH mRasterBandH;
 
+    //! 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