use QGIS API to write raster files in QgsGridFileWriter (fix #56653)

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,67 @@ 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 )
+  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 );
+
+    std::vector<float> float32Row( mNumColumns );
+    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/test_analysis_interpolation.py

@@ -74,7 +74,7 @@ def test_idw_interpolator(self):
         interpolator = QgsIDWInterpolator([data])
         interpolator.setDistanceCoefficient(2.0)
 
-        output_file = os.path.join(tempfile.gettempdir(), "idw_interpolation.asc")
+        output_file = os.path.join(tempfile.gettempdir(), "idw_interpolation.tif")
 
         writer = QgsGridFileWriter(interpolator, output_file, extent, cols, rows)
         writer.writeFile()
@@ -84,7 +84,7 @@ def test_idw_interpolator(self):
             "gdal",
             output_file,
             "gdal",
-            os.path.join(TEST_DATA_DIR, "analysis", "idw_interpolation.asc"),
+            os.path.join(TEST_DATA_DIR, "analysis", "idw_interpolation.tif"),
         )
         self.report += checker.report()
 
@@ -118,7 +118,7 @@ def test_tin_interpolator(self):
             [data], QgsTinInterpolator.TinInterpolation.Linear
         )
 
-        output_file = os.path.join(tempfile.gettempdir(), "tin_interpolation.asc")
+        output_file = os.path.join(tempfile.gettempdir(), "tin_interpolation.tif")
 
         writer = QgsGridFileWriter(interpolator, output_file, extent, cols, rows)
         writer.writeFile()
@@ -128,7 +128,7 @@ def test_tin_interpolator(self):
             "gdal",
             output_file,
             "gdal",
-            os.path.join(TEST_DATA_DIR, "analysis", "tin_interpolation.asc"),
+            os.path.join(TEST_DATA_DIR, "analysis", "tin_interpolation.tif"),
         )
         self.report += checker.report()