[Backport release-3_40] Use native raster API in QgsGridFileWriter

python/plugins/processing/tests/testdata/qgis_algorithm_tests2.yaml

@@ -14,9 +14,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - bfb3616a73065c0cb41eb1c5c1e9e8812fe1c63019a6177adb1dfe3e
-          - 19bbd79d15b0ba7dbde05665d559e9806fdb1cd579df9222f9cc4d92
-          - 3dfa1e041f33a72abf72d2d706c03e11fa5274f152dec1db1661b8e0
+          - e277a6774c489ee965421fe102780275163229b9407d5e33d3e268e5
         type: rasterhash
 
   - algorithm: qgis:tininterpolation
@@ -31,9 +29,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - bfb3616a73065c0cb41eb1c5c1e9e8812fe1c63019a6177adb1dfe3e
-          - 19bbd79d15b0ba7dbde05665d559e9806fdb1cd579df9222f9cc4d92
-          - 3dfa1e041f33a72abf72d2d706c03e11fa5274f152dec1db1661b8e0
+          - e277a6774c489ee965421fe102780275163229b9407d5e33d3e268e5
         type: rasterhash
 
   - algorithm: qgis:tininterpolation
@@ -48,9 +44,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - b8e49813c507b73cb39a5ad904b2d302ccb25c591a2e577b6405f982
-          - 19bbd79d15b0ba7dbde05665d559e9806fdb1cd579df9222f9cc4d92
-          - 3dfa1e041f33a72abf72d2d706c03e11fa5274f152dec1db1661b8e0
+          - e277a6774c489ee965421fe102780275163229b9407d5e33d3e268e5
         type: rasterhash
 
   - algorithm: qgis:tininterpolation
@@ -66,9 +60,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - bfb3616a73065c0cb41eb1c5c1e9e8812fe1c63019a6177adb1dfe3e
-          - 19bbd79d15b0ba7dbde05665d559e9806fdb1cd579df9222f9cc4d92
-          - 3dfa1e041f33a72abf72d2d706c03e11fa5274f152dec1db1661b8e0
+          - e277a6774c489ee965421fe102780275163229b9407d5e33d3e268e5
         type: rasterhash
 
   - algorithm: qgis:tininterpolation
@@ -84,9 +76,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - bfb3616a73065c0cb41eb1c5c1e9e8812fe1c63019a6177adb1dfe3e
-          - 19bbd79d15b0ba7dbde05665d559e9806fdb1cd579df9222f9cc4d92
-          - 3dfa1e041f33a72abf72d2d706c03e11fa5274f152dec1db1661b8e0
+          - e277a6774c489ee965421fe102780275163229b9407d5e33d3e268e5
         type: rasterhash
 
   - algorithm: qgis:tininterpolation
@@ -102,9 +92,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - b8e49813c507b73cb39a5ad904b2d302ccb25c591a2e577b6405f982
-          - 19bbd79d15b0ba7dbde05665d559e9806fdb1cd579df9222f9cc4d92
-          - 3dfa1e041f33a72abf72d2d706c03e11fa5274f152dec1db1661b8e0
+          - e277a6774c489ee965421fe102780275163229b9407d5e33d3e268e5
         type: rasterhash
 
   - algorithm: qgis:idwinterpolation
@@ -119,8 +107,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - 55dfd0e0e5ea300d2314785ab6aa1f4d941a4506e63d0d9e2b84e959
-          - b5fc06549f1b0ce3261f9fb8868bb6d082edcbfe409c89a82e2b7e7a
+          - 7836f01a89afbfdd865aaae9b57c602c96d629184adf1b6ea66ca127
         type: rasterhash
 
   - algorithm: qgis:idwinterpolation
@@ -135,8 +122,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - 55dfd0e0e5ea300d2314785ab6aa1f4d941a4506e63d0d9e2b84e959
-          - b5fc06549f1b0ce3261f9fb8868bb6d082edcbfe409c89a82e2b7e7a
+          - 7836f01a89afbfdd865aaae9b57c602c96d629184adf1b6ea66ca127
         type: rasterhash
 
   - algorithm: qgis:idwinterpolation
@@ -151,8 +137,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - d53cd7327d53bd4b8fb2de8b206e5aa0a1366963f2b63f58c2b7926e
-          - b5fc06549f1b0ce3261f9fb8868bb6d082edcbfe409c89a82e2b7e7a
+          - 7836f01a89afbfdd865aaae9b57c602c96d629184adf1b6ea66ca127
         type: rasterhash
 
   - algorithm: qgis:idwinterpolation
@@ -168,8 +153,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - 55dfd0e0e5ea300d2314785ab6aa1f4d941a4506e63d0d9e2b84e959
-          - b5fc06549f1b0ce3261f9fb8868bb6d082edcbfe409c89a82e2b7e7a
+          - 7836f01a89afbfdd865aaae9b57c602c96d629184adf1b6ea66ca127
         type: rasterhash
 
   - algorithm: qgis:idwinterpolation
@@ -185,8 +169,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - 55dfd0e0e5ea300d2314785ab6aa1f4d941a4506e63d0d9e2b84e959
-          - b5fc06549f1b0ce3261f9fb8868bb6d082edcbfe409c89a82e2b7e7a
+          - 7836f01a89afbfdd865aaae9b57c602c96d629184adf1b6ea66ca127
         type: rasterhash
 
   - algorithm: qgis:idwinterpolation
@@ -202,8 +185,7 @@ tests:
     results:
       OUTPUT:
         hash:
-          - d53cd7327d53bd4b8fb2de8b206e5aa0a1366963f2b63f58c2b7926e
-          - b5fc06549f1b0ce3261f9fb8868bb6d082edcbfe409c89a82e2b7e7a
+          - 7836f01a89afbfdd865aaae9b57c602c96d629184adf1b6ea66ca127
         type: rasterhash
 
   - algorithm: native:removenullgeometries

src/analysis/interpolation/qgsgridfilewriter.cpp

@@ -19,7 +19,9 @@
 #include "qgsinterpolator.h"
 #include "qgsvectorlayer.h"
 #include "qgsfeedback.h"
-#include <QFile>
+#include "qgsrasterfilewriter.h"
+#include "qgsrasterdataprovider.h"
+#include "qgsrasterblock.h"
 #include <QFileInfo>
 
 QgsGridFileWriter::QgsGridFileWriter( QgsInterpolator *i, const QString &outputPath, const QgsRectangle &extent, int nCols, int nRows )
@@ -34,95 +36,66 @@ QgsGridFileWriter::QgsGridFileWriter( QgsInterpolator *i, const QString &outputP
 
 int QgsGridFileWriter::writeFile( QgsFeedback *feedback )
 {
-  QFile outputFile( mOutputFilePath );
+  const QFileInfo fi( mOutputFilePath );
+  const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
+
+  QgsInterpolator::LayerData ld = mInterpolator->layerData().at( 0 );
+  const QgsCoordinateReferenceSystem crs = ld.source->sourceCrs();
 
-  if ( !outputFile.open( QFile::WriteOnly | QIODevice::Truncate ) )
+  std::unique_ptr<QgsRasterFileWriter> writer = std::make_unique<QgsRasterFileWriter>( mOutputFilePath );
+  writer->setOutputProviderKey( QStringLiteral( "gdal" ) );
+  writer->setOutputFormat( outputFormat );
+
+  std::unique_ptr<QgsRasterDataProvider> provider( writer->createOneBandRaster( Qgis::DataType::Float32, mNumColumns, mNumRows, mInterpolationExtent, crs ) );
+  if ( !provider )
   {
+    QgsDebugMsgLevel( QStringLiteral( "Could not create raster output: %1" ).arg( mOutputFilePath ), 2 );
     return 1;
   }
-
-  if ( !mInterpolator )
+  if ( !provider->isValid() )
   {
-    outputFile.remove();
+    QgsDebugMsgLevel( QStringLiteral( "Could not create raster output: %1: %2" ).arg( mOutputFilePath, provider->error().message( QgsErrorMessage::Text ) ), 2 );
     return 2;
   }
 
-  QTextStream outStream( &outputFile );
-#if QT_VERSION < QT_VERSION_CHECK( 6, 0, 0 )
-  outStream.setCodec( "UTF-8" );
-#endif
-  outStream.setRealNumberPrecision( 8 );
-  writeHeader( outStream );
+  provider->setNoDataValue( 1, -9999 );
 
   double currentYValue = mInterpolationExtent.yMaximum() - mCellSizeY / 2.0; //calculate value in the center of the cell
   double currentXValue;
   double interpolatedValue;
 
-  for ( int i = 0; i < mNumRows; ++i )
+  std::vector<float> float32Row( mNumColumns );
+  const double step = mNumRows > 0 ? 100.0 / mNumRows : 1;
+  for ( int row = 0; row < mNumRows; row++ )
   {
+    if ( feedback && feedback->isCanceled() )
+    {
+      break;
+    }
+
     currentXValue = mInterpolationExtent.xMinimum() + mCellSizeX / 2.0; //calculate value in the center of the cell
-    for ( int j = 0; j < mNumColumns; ++j )
+    QgsRasterBlock block( Qgis::DataType::Float32, mNumColumns, 1 );
+
+    for ( int col = 0; col < mNumColumns; col++ )
     {
       if ( mInterpolator->interpolatePoint( currentXValue, currentYValue, interpolatedValue, feedback ) == 0 )
       {
-        outStream << interpolatedValue << ' ';
+        float32Row[col] = interpolatedValue;
       }
       else
       {
-        outStream << "-9999 ";
+        float32Row[col] = -9999;
       }
       currentXValue += mCellSizeX;
     }
-
-    outStream << Qt::endl;
+    block.setData( QByteArray( reinterpret_cast<const char *>( float32Row.data() ), QgsRasterBlock::typeSize( Qgis::DataType::Float32 ) * mNumColumns ) );
+    provider->writeBlock( &block, 1, 0, row );
     currentYValue -= mCellSizeY;
-
     if ( feedback )
     {
-      if ( feedback->isCanceled() )
-      {
-        outputFile.remove();
-        return 3;
-      }
-      feedback->setProgress( 100.0 * i / static_cast<double>( mNumRows ) );
+      feedback->setProgress( row * step );
     }
   }
 
-  // create prj file
-  QgsInterpolator::LayerData ld;
-  ld = mInterpolator->layerData().at( 0 );
-  QgsFeatureSource *source = ld.source;
-  const QString crs = source->sourceCrs().toWkt();
-  const QFileInfo fi( mOutputFilePath );
-  const QString fileName = fi.absolutePath() + '/' + fi.completeBaseName() + ".prj";
-  QFile prjFile( fileName );
-  if ( !prjFile.open( QFile::WriteOnly | QIODevice::Truncate ) )
-  {
-    return 1;
-  }
-  QTextStream prjStream( &prjFile );
-  prjStream << crs;
-  prjStream << Qt::endl;
-  prjFile.close();
-
-  return 0;
-}
-
-int QgsGridFileWriter::writeHeader( QTextStream &outStream )
-{
-  outStream << "NCOLS " << mNumColumns << Qt::endl;
-  outStream << "NROWS " << mNumRows << Qt::endl;
-  outStream << "XLLCORNER " << mInterpolationExtent.xMinimum() << Qt::endl;
-  outStream << "YLLCORNER " << mInterpolationExtent.yMinimum() << Qt::endl;
-  if ( mCellSizeX == mCellSizeY ) //standard way
-  {
-    outStream << "CELLSIZE " << mCellSizeX << Qt::endl;
-  }
-  else //this is supported by GDAL but probably not by other products
-  {
-    outStream << "DX " << mCellSizeX << Qt::endl;
-    outStream << "DY " << mCellSizeY << Qt::endl;
-  }
-  outStream << "NODATA_VALUE -9999" << Qt::endl;
   return 0;
 }

src/analysis/interpolation/qgsgridfilewriter.h

@@ -56,8 +56,6 @@ class ANALYSIS_EXPORT QgsGridFileWriter
   private:
     QgsGridFileWriter() = delete;
 
-    int writeHeader( QTextStream &outStream );
-
     QgsInterpolator *mInterpolator = nullptr;
     QString mOutputFilePath;
     QgsRectangle mInterpolationExtent;

tests/src/python/CMakeLists.txt

@@ -403,6 +403,7 @@ ADD_PYTHON_TEST(PyQgsAttributeEditorAction test_qgsattributeeditoraction.py)
 ADD_PYTHON_TEST(PyQgsVectorTile test_qgsvectortile.py)
 ADD_PYTHON_TEST(PyQgsVtpk test_qgsvtpk.py)
 ADD_PYTHON_TEST(PyQgsProcessingAlgsGdalGdalUtils test_processing_algs_gdal_gdalutils.py)
+ADD_PYTHON_TEST(PyQgsInterpolation test_analysis_interpolation.py)
 
 if (NOT WIN32)
   ADD_PYTHON_TEST(PyQgsLogger test_qgslogger.py)