[GDAL provider] in create(), use newly create dataset handle to insta…

…nciate the provider (fix #17103) Previously we created an empty file, close it and re-opened it immediately in open mode. However for GeoPackage, if you create for example a 1 band or 3 band dataset, upon reopening you get a 4 band RGBA dataset, due to the fact that GeoPackage doesn't explicitly store the number of bands. Thus the 4th band (alpha) was never written by the QgsRasterFileWriter logic.
qgis · Nov 9, 2017 · c3e15f5 · c3e15f5
1 parent 75c41f5
commit c3e15f5
Show file tree

Hide file tree

Showing 3 changed files with 68 additions and 26 deletions.
diff --git a/src/providers/gdal/qgsgdalprovider.cpp b/src/providers/gdal/qgsgdalprovider.cpp
@@ -117,7 +117,7 @@ QgsGdalProvider::QgsGdalProvider( const QString &uri, const QgsError &error )
   setError( error );
 }
 
-QgsGdalProvider::QgsGdalProvider( const QString &uri, bool update )
+QgsGdalProvider::QgsGdalProvider( const QString &uri, bool update, GDALDatasetH dataset )
   : QgsRasterDataProvider( uri )
   , mUpdate( update )
 {
@@ -156,29 +156,35 @@ QgsGdalProvider::QgsGdalProvider( const QString &uri, bool update )
   }
 
   mGdalDataset = nullptr;
-
-  // Try to open using VSIFileHandler (see qgsogrprovider.cpp)
-  QString vsiPrefix = QgsZipItem::vsiPrefix( uri );
-  if ( !vsiPrefix.isEmpty() )
+  if ( dataset )
   {
-    if ( !uri.startsWith( vsiPrefix ) )
-      setDataSourceUri( vsiPrefix + uri );
-    QgsDebugMsg( QString( "Trying %1 syntax, uri= %2" ).arg( vsiPrefix, dataSourceUri() ) );
+    mGdalBaseDataset = dataset;
   }
+  else
+  {
+    // Try to open using VSIFileHandler (see qgsogrprovider.cpp)
+    QString vsiPrefix = QgsZipItem::vsiPrefix( uri );
+    if ( !vsiPrefix.isEmpty() )
+    {
+      if ( !uri.startsWith( vsiPrefix ) )
+        setDataSourceUri( vsiPrefix + uri );
+      QgsDebugMsg( QString( "Trying %1 syntax, uri= %2" ).arg( vsiPrefix, dataSourceUri() ) );
+    }
 
-  QString gdalUri = dataSourceUri();
+    QString gdalUri = dataSourceUri();
 
-  CPLErrorReset();
-  mGdalBaseDataset = gdalOpen( gdalUri.toUtf8().constData(), mUpdate ? GA_Update : GA_ReadOnly );
+    CPLErrorReset();
+    mGdalBaseDataset = gdalOpen( gdalUri.toUtf8().constData(), mUpdate ? GA_Update : GA_ReadOnly );
 
-  if ( !mGdalBaseDataset )
-  {
-    QString msg = QStringLiteral( "Cannot open GDAL dataset %1:\n%2" ).arg( dataSourceUri(), QString::fromUtf8( CPLGetLastErrorMsg() ) );
-    appendError( ERRMSG( msg ) );
-    return;
-  }
+    if ( !mGdalBaseDataset )
+    {
+      QString msg = QStringLiteral( "Cannot open GDAL dataset %1:\n%2" ).arg( dataSourceUri(), QString::fromUtf8( CPLGetLastErrorMsg() ) );
+      appendError( ERRMSG( msg ) );
+      return;
+    }
 
-  QgsDebugMsg( "GdalDataset opened" );
+    QgsDebugMsg( "GdalDataset opened" );
+  }
   initBaseDataset();
 }
 
@@ -2726,7 +2732,7 @@ QGISEXTERN QgsGdalProvider *create(
   //create dataset
   CPLErrorReset();
   char **papszOptions = papszFromStringList( createOptions );
-  gdal::dataset_unique_ptr dataset( GDALCreate( driver, uri.toUtf8().constData(), width, height, nBands, ( GDALDataType )type, papszOptions ) );
+  GDALDatasetH dataset = GDALCreate( driver, uri.toUtf8().constData(), width, height, nBands, ( GDALDataType )type, papszOptions );
   CSLDestroy( papszOptions );
   if ( !dataset )
   {
@@ -2735,11 +2741,10 @@ QGISEXTERN QgsGdalProvider *create(
     return new QgsGdalProvider( uri, error );
   }
 
-  GDALSetGeoTransform( dataset.get(), geoTransform );
-  GDALSetProjection( dataset.get(), crs.toWkt().toLocal8Bit().data() );
-  dataset.reset();
+  GDALSetGeoTransform( dataset, geoTransform );
+  GDALSetProjection( dataset, crs.toWkt().toLocal8Bit().data() );
 
-  return new QgsGdalProvider( uri, true );
+  return new QgsGdalProvider( uri, true, dataset );
 }
 
 bool QgsGdalProvider::write( void *data, int band, int width, int height, int xOffset, int yOffset )

diff --git a/src/providers/gdal/qgsgdalprovider.h b/src/providers/gdal/qgsgdalprovider.h
@@ -63,11 +63,12 @@ class QgsGdalProvider : public QgsRasterDataProvider, QgsGdalProviderBase
     /**
      * Constructor for the provider.
      *
-     * \param   uri   HTTP URL of the Web Server.  If needed a proxy will be used
-     *                otherwise we contact the host directly.
+     * \param   uri         file name
+     * \param   update      whether to open in update mode
+     * \param   newDataset  handle of newly created dataset.
      *
      */
-    QgsGdalProvider( QString const &uri = QString(), bool update = false );
+    QgsGdalProvider( QString const &uri = QString(), bool update = false, GDALDatasetH newDataset = nullptr );
 
     //! Create invalid provider with error
     QgsGdalProvider( QString const &uri, const QgsError &error );

diff --git a/tests/src/python/test_qgsrasterfilewriter.py b/tests/src/python/test_qgsrasterfilewriter.py
@@ -16,7 +16,9 @@
 
 import os
 import glob
+import tempfile
 
+from osgeo import gdal
 from qgis.PyQt.QtCore import QTemporaryFile, QDir
 from qgis.core import (QgsRasterLayer,
                        QgsRasterChecker,
@@ -107,6 +109,40 @@ def testDriverForExtension(self):
         self.assertEqual(QgsRasterFileWriter.driverForExtension('not a format'), '')
         self.assertEqual(QgsRasterFileWriter.driverForExtension(''), '')
 
+    def testImportIntoGpkg(self):
+        # init target file
+        test_gpkg = tempfile.mktemp(suffix='.gpkg', dir=self.testDataDir)
+        gdal.GetDriverByName('GPKG').Create(test_gpkg, 1, 1, 1)
+        source = QgsRasterLayer(os.path.join(self.testDataDir, 'raster', 'band3_byte_noct_epsg4326.tif'), 'my', 'gdal')
+        self.assertTrue(source.isValid())
+        provider = source.dataProvider()
+        fw = QgsRasterFileWriter(test_gpkg)
+        fw.setOutputFormat('gpkg')
+        fw.setCreateOptions(['RASTER_TABLE=imported_table', 'APPEND_SUBDATASET=YES'])
+
+        pipe = QgsRasterPipe()
+        self.assertTrue(pipe.set(provider.clone()))
+
+        projector = QgsRasterProjector()
+        projector.setCrs(provider.crs(), provider.crs())
+        self.assertTrue(pipe.insert(2, projector))
+
+        self.assertEqual(fw.writeRaster(pipe,
+                                        provider.xSize(),
+                                        provider.ySize(),
+                                        provider.extent(),
+                                        provider.crs()), 0)
+
+        # Check that the test geopackage contains the raster layer and compare
+        rlayer = QgsRasterLayer('GPKG:%s:imported_table' % test_gpkg)
+        self.assertTrue(rlayer.isValid())
+        out_provider = rlayer.dataProvider()
+        self.assertEqual(provider.block(1, provider.extent(), source.width(), source.height()).data(),
+                         out_provider.block(1, out_provider.extent(), rlayer.width(), rlayer.height()).data())
+
+        # remove result file
+        os.unlink(test_gpkg)
+
 
 if __name__ == '__main__':
     unittest.main()