QGIS API Documentation  2.99.0-Master (8ec3eaf)
qgsrasterfilewriter.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrasterfilewriter.cpp
3  ---------------------
4  begin : July 2012
5  copyright : (C) 2012 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 #include <typeinfo>
16 
17 #include "qgsrasterfilewriter.h"
18 #include "qgscoordinatetransform.h"
19 #include "qgsproviderregistry.h"
20 #include "qgsrasterinterface.h"
21 #include "qgsrasteriterator.h"
22 #include "qgsrasterlayer.h"
23 #include "qgsrasterprojector.h"
24 #include "qgsrasterdataprovider.h"
25 #include "qgsrasternuller.h"
26 
27 #include <QCoreApplication>
28 #include <QProgressDialog>
29 #include <QTextStream>
30 #include <QMessageBox>
31 
32 QgsRasterFileWriter::QgsRasterFileWriter( const QString& outputUrl )
33  : mMode( Raw )
34  , mOutputUrl( outputUrl )
35  , mOutputProviderKey( QStringLiteral( "gdal" ) )
36  , mOutputFormat( QStringLiteral( "GTiff" ) )
37  , mTiledMode( false )
38  , mMaxTileWidth( 500 )
39  , mMaxTileHeight( 500 )
40  , mBuildPyramidsFlag( QgsRaster::PyramidsFlagNo )
41  , mPyramidsFormat( QgsRaster::PyramidsGTiff )
42  , mProgressDialog( nullptr )
43  , mPipe( nullptr )
44  , mInput( nullptr )
45 {
46 
47 }
48 
50  : mMode( Raw )
51  , mOutputProviderKey( QStringLiteral( "gdal" ) )
52  , mOutputFormat( QStringLiteral( "GTiff" ) )
53  , mTiledMode( false )
54  , mMaxTileWidth( 500 )
55  , mMaxTileHeight( 500 )
56  , mBuildPyramidsFlag( QgsRaster::PyramidsFlagNo )
57  , mPyramidsFormat( QgsRaster::PyramidsGTiff )
58  , mProgressDialog( nullptr )
59  , mPipe( nullptr )
60  , mInput( nullptr )
61 {
62 
63 }
64 
65 QgsRasterFileWriter::WriterError QgsRasterFileWriter::writeRaster( const QgsRasterPipe* pipe, int nCols, int nRows, const QgsRectangle& outputExtent,
66  const QgsCoordinateReferenceSystem& crs, QProgressDialog* progressDialog )
67 {
68  QgsDebugMsgLevel( "Entered", 4 );
69 
70  if ( !pipe )
71  {
72  return SourceProviderError;
73  }
74  mPipe = pipe;
75 
76  //const QgsRasterInterface* iface = iter->input();
77  const QgsRasterInterface* iface = pipe->last();
78  if ( !iface )
79  {
80  return SourceProviderError;
81  }
82  mInput = iface;
83 
84  if ( QgsRasterBlock::typeIsColor( iface->dataType( 1 ) ) )
85  {
86  mMode = Image;
87  }
88  else
89  {
90  mMode = Raw;
91  }
92 
93  QgsDebugMsgLevel( QString( "reading from %1" ).arg( typeid( *iface ).name() ), 4 );
94 
95  if ( !iface->sourceInput() )
96  {
97  QgsDebugMsg( "iface->srcInput() == 0" );
98  return SourceProviderError;
99  }
100 #ifdef QGISDEBUG
101  const QgsRasterInterface &srcInput = *iface->sourceInput();
102  QgsDebugMsgLevel( QString( "srcInput = %1" ).arg( typeid( srcInput ).name() ), 4 );
103 #endif
104 
105  mProgressDialog = progressDialog;
106 
107  QgsRasterIterator iter( pipe->last() );
108 
109  //create directory for output files
110  if ( mTiledMode )
111  {
112  QFileInfo fileInfo( mOutputUrl );
113  if ( !fileInfo.exists() )
114  {
115  QDir dir = fileInfo.dir();
116  if ( !dir.mkdir( fileInfo.fileName() ) )
117  {
118  QgsDebugMsg( "Cannot create output VRT directory " + fileInfo.fileName() + " in " + dir.absolutePath() );
119  return CreateDatasourceError;
120  }
121  }
122  }
123 
124  if ( mMode == Image )
125  {
126  WriterError e = writeImageRaster( &iter, nCols, nRows, outputExtent, crs, progressDialog );
127  mProgressDialog = nullptr;
128  return e;
129  }
130  else
131  {
132  mProgressDialog = nullptr;
133  WriterError e = writeDataRaster( pipe, &iter, nCols, nRows, outputExtent, crs, progressDialog );
134  return e;
135  }
136 }
137 
138 QgsRasterFileWriter::WriterError QgsRasterFileWriter::writeDataRaster( const QgsRasterPipe* pipe, QgsRasterIterator* iter, int nCols, int nRows, const QgsRectangle& outputExtent,
139  const QgsCoordinateReferenceSystem& crs, QProgressDialog* progressDialog )
140 {
141  QgsDebugMsgLevel( "Entered", 4 );
142  if ( !iter )
143  {
144  return SourceProviderError;
145  }
146 
147  const QgsRasterInterface* iface = pipe->last();
148  if ( !iface )
149  {
150  return SourceProviderError;
151  }
152 
153  QgsRasterDataProvider* srcProvider = const_cast<QgsRasterDataProvider*>( dynamic_cast<const QgsRasterDataProvider*>( iface->sourceInput() ) );
154  if ( !srcProvider )
155  {
156  QgsDebugMsg( "Cannot get source data provider" );
157  return SourceProviderError;
158  }
159 
160  iter->setMaximumTileWidth( mMaxTileWidth );
161  iter->setMaximumTileHeight( mMaxTileHeight );
162 
163  int nBands = iface->bandCount();
164  if ( nBands < 1 )
165  {
166  return SourceProviderError;
167  }
168 
169 
170  //check if all the bands have the same data type size, otherwise we cannot write it to the provider
171  //(at least not with the current interface)
172  int dataTypeSize = QgsRasterBlock::typeSize( srcProvider->sourceDataType( 1 ) );
173  for ( int i = 2; i <= nBands; ++i )
174  {
175  if ( QgsRasterBlock::typeSize( srcProvider->sourceDataType( 1 ) ) != dataTypeSize )
176  {
177  return DestProviderError;
178  }
179  }
180 
181  // Output data type - source data type is preferred but it may happen that we need
182  // to set 'no data' value (which was not set on source data) if output extent
183  // is larger than source extent (with or without reprojection) and there is no 'free'
184  // (not used) value available
185  QList<bool> destHasNoDataValueList;
186  QList<double> destNoDataValueList;
187  QList<Qgis::DataType> destDataTypeList;
188  destDataTypeList.reserve( nBands );
189  destHasNoDataValueList.reserve( nBands );
190  destNoDataValueList.reserve( nBands );
191 
192  for ( int bandNo = 1; bandNo <= nBands; bandNo++ )
193  {
194  QgsRasterNuller *nuller = pipe->nuller();
195 
196  bool srcHasNoDataValue = srcProvider->sourceHasNoDataValue( bandNo );
197  bool destHasNoDataValue = false;
198  double destNoDataValue = std::numeric_limits<double>::quiet_NaN();
199  Qgis::DataType destDataType = srcProvider->sourceDataType( bandNo );
200  // TODO: verify what happens/should happen if srcNoDataValue is disabled by setUseSrcNoDataValue
201  QgsDebugMsgLevel( QString( "srcHasNoDataValue = %1 srcNoDataValue = %2" ).arg( srcHasNoDataValue ).arg( srcProvider->sourceNoDataValue( bandNo ) ), 4 );
202  if ( srcHasNoDataValue )
203  {
204 
205  // If source has no data value, it is used by provider
206  destNoDataValue = srcProvider->sourceNoDataValue( bandNo );
207  destHasNoDataValue = true;
208  }
209  else if ( nuller && !nuller->noData( bandNo ).isEmpty() )
210  {
211  // Use one user defined no data value
212  destNoDataValue = nuller->noData( bandNo ).value( 0 ).min();
213  destHasNoDataValue = true;
214  }
215  else
216  {
217  // Verify if we realy need no data value, i.e.
218  QgsRectangle srcExtent = outputExtent;
219  QgsRasterProjector *projector = pipe->projector();
220  if ( projector && projector->destinationCrs() != projector->sourceCrs() )
221  {
222  QgsCoordinateTransform ct( projector->destinationCrs(), projector->sourceCrs() );
223  srcExtent = ct.transformBoundingBox( outputExtent );
224  }
225  if ( !srcProvider->extent().contains( srcExtent ) )
226  {
227  // Destination extent is larger than source extent, we need destination no data values
228  // Get src sample statistics (estimation from sample)
229  QgsRasterBandStats stats = srcProvider->bandStatistics( bandNo, QgsRasterBandStats::Min | QgsRasterBandStats::Max, srcExtent, 250000 );
230 
231  // Test if we have free (not used) values
232  double typeMinValue = QgsContrastEnhancement::maximumValuePossible( static_cast< Qgis::DataType >( srcProvider->sourceDataType( bandNo ) ) );
233  double typeMaxValue = QgsContrastEnhancement::maximumValuePossible( static_cast< Qgis::DataType >( srcProvider->sourceDataType( bandNo ) ) );
234  if ( stats.minimumValue > typeMinValue )
235  {
236  destNoDataValue = typeMinValue;
237  }
238  else if ( stats.maximumValue < typeMaxValue )
239  {
240  destNoDataValue = typeMaxValue;
241  }
242  else
243  {
244  // We have to use wider type
245  destDataType = QgsRasterBlock::typeWithNoDataValue( destDataType, &destNoDataValue );
246  }
247  destHasNoDataValue = true;
248  }
249  }
250 
251  if ( nuller && destHasNoDataValue )
252  {
253  nuller->setOutputNoDataValue( bandNo, destNoDataValue );
254  }
255 
256  QgsDebugMsgLevel( QString( "bandNo = %1 destDataType = %2 destHasNoDataValue = %3 destNoDataValue = %4" ).arg( bandNo ).arg( destDataType ).arg( destHasNoDataValue ).arg( destNoDataValue ), 4 );
257  destDataTypeList.append( destDataType );
258  destHasNoDataValueList.append( destHasNoDataValue );
259  destNoDataValueList.append( destNoDataValue );
260  }
261 
262  Qgis::DataType destDataType = destDataTypeList.value( 0 );
263  // Currently write API supports one output type for dataset only -> find the widest
264  for ( int i = 1; i < nBands; i++ )
265  {
266  if ( destDataTypeList.value( i ) > destDataType )
267  {
268  destDataType = destDataTypeList.value( i );
269  // no data value may be left per band (for future)
270  }
271  }
272 
273  //create destProvider for whole dataset here
274  QgsRasterDataProvider* destProvider = nullptr;
275  double pixelSize;
276  double geoTransform[6];
277  globalOutputParameters( outputExtent, nCols, nRows, geoTransform, pixelSize );
278 
279  // initOutput() returns 0 in tile mode!
280  destProvider = initOutput( nCols, nRows, crs, geoTransform, nBands, destDataType, destHasNoDataValueList, destNoDataValueList );
281 
282  WriterError error = writeDataRaster( pipe, iter, nCols, nRows, outputExtent, crs, destDataType, destHasNoDataValueList, destNoDataValueList, destProvider, progressDialog );
283 
284  if ( error == NoDataConflict )
285  {
286  // The value used for no data was found in source data, we must use wider data type
287  if ( destProvider ) // no tiles
288  {
289  destProvider->remove();
290  delete destProvider;
291  destProvider = nullptr;
292  }
293  else // VRT
294  {
295  // TODO: remove created VRT
296  }
297 
298  // But we don't know which band -> wider all
299  for ( int i = 0; i < nBands; i++ )
300  {
301  double destNoDataValue;
302  Qgis::DataType destDataType = QgsRasterBlock::typeWithNoDataValue( destDataTypeList.value( i ), &destNoDataValue );
303  destDataTypeList.replace( i, destDataType );
304  destNoDataValueList.replace( i, destNoDataValue );
305  }
306  destDataType = destDataTypeList.value( 0 );
307 
308  // Try again
309  destProvider = initOutput( nCols, nRows, crs, geoTransform, nBands, destDataType, destHasNoDataValueList, destNoDataValueList );
310  error = writeDataRaster( pipe, iter, nCols, nRows, outputExtent, crs, destDataType, destHasNoDataValueList, destNoDataValueList, destProvider, progressDialog );
311  }
312 
313  if ( destProvider )
314  delete destProvider;
315 
316  return error;
317 }
318 
319 QgsRasterFileWriter::WriterError QgsRasterFileWriter::writeDataRaster(
320  const QgsRasterPipe* pipe,
321  QgsRasterIterator* iter,
322  int nCols, int nRows,
323  const QgsRectangle& outputExtent,
324  const QgsCoordinateReferenceSystem& crs,
325  Qgis::DataType destDataType,
326  const QList<bool>& destHasNoDataValueList,
327  const QList<double>& destNoDataValueList,
328  QgsRasterDataProvider* destProvider,
329  QProgressDialog* progressDialog )
330 {
331  Q_UNUSED( pipe );
332  Q_UNUSED( destHasNoDataValueList );
333  QgsDebugMsgLevel( "Entered", 4 );
334 
335  const QgsRasterInterface* iface = iter->input();
336  const QgsRasterDataProvider *srcProvider = dynamic_cast<const QgsRasterDataProvider*>( iface->sourceInput() );
337  int nBands = iface->bandCount();
338  QgsDebugMsgLevel( QString( "nBands = %1" ).arg( nBands ), 4 );
339 
340  //Get output map units per pixel
341  int iterLeft = 0;
342  int iterTop = 0;
343  int iterCols = 0;
344  int iterRows = 0;
345 
346  QList<QgsRasterBlock*> blockList;
347  blockList.reserve( nBands );
348  for ( int i = 1; i <= nBands; ++i )
349  {
350  iter->startRasterRead( i, nCols, nRows, outputExtent );
351  blockList.push_back( nullptr );
352  if ( destProvider && destHasNoDataValueList.value( i - 1 ) ) // no tiles
353  {
354  destProvider->setNoDataValue( i, destNoDataValueList.value( i - 1 ) );
355  }
356  }
357 
358  int nParts = 0;
359  int fileIndex = 0;
360  if ( progressDialog )
361  {
362  int nPartsX = nCols / iter->maximumTileWidth() + 1;
363  int nPartsY = nRows / iter->maximumTileHeight() + 1;
364  nParts = nPartsX * nPartsY;
365  progressDialog->setMaximum( nParts );
366  progressDialog->show();
367  progressDialog->setLabelText( QObject::tr( "Reading raster part %1 of %2" ).arg( fileIndex + 1 ).arg( nParts ) );
368  }
369 
370  // hmm why is there a for(;;) here ..
371  // not good coding practice IMHO, it might be better to use [ for() and break ] or [ while (test) ]
372  Q_FOREVER
373  {
374  for ( int i = 1; i <= nBands; ++i )
375  {
376  if ( !iter->readNextRasterPart( i, iterCols, iterRows, &( blockList[i - 1] ), iterLeft, iterTop ) )
377  {
378  // No more parts, create VRT and return
379  if ( mTiledMode )
380  {
381  QString vrtFilePath( mOutputUrl + '/' + vrtFileName() );
382  writeVRT( vrtFilePath );
383  if ( mBuildPyramidsFlag == QgsRaster::PyramidsFlagYes )
384  {
385  buildPyramids( vrtFilePath );
386  }
387  }
388  else
389  {
390  if ( mBuildPyramidsFlag == QgsRaster::PyramidsFlagYes )
391  {
392  buildPyramids( mOutputUrl );
393  }
394  }
395 
396  QgsDebugMsgLevel( "Done", 4 );
397  return NoError; //reached last tile, bail out
398  }
399  // TODO: verify if NoDataConflict happened, to do that we need the whole pipe or nuller interface
400  }
401 
402  if ( progressDialog && fileIndex < ( nParts - 1 ) )
403  {
404  progressDialog->setValue( fileIndex + 1 );
405  progressDialog->setLabelText( QObject::tr( "Reading raster part %1 of %2" ).arg( fileIndex + 2 ).arg( nParts ) );
406  QCoreApplication::processEvents( QEventLoop::AllEvents, 1000 );
407  if ( progressDialog->wasCanceled() )
408  {
409  for ( int i = 0; i < nBands; ++i )
410  {
411  delete blockList[i];
412  }
413  break;
414  }
415  }
416 
417  // It may happen that internal data type (dataType) is wider than destDataType
418  QList<QgsRasterBlock*> destBlockList;
419  for ( int i = 1; i <= nBands; ++i )
420  {
421  if ( srcProvider && srcProvider->dataType( i ) == destDataType )
422  {
423  destBlockList.push_back( blockList[i-1] );
424  }
425  else
426  {
427  // TODO: this conversion should go to QgsRasterDataProvider::write with additional input data type param
428  blockList[i-1]->convert( destDataType );
429  destBlockList.push_back( blockList[i-1] );
430  }
431  blockList[i-1] = nullptr;
432  }
433 
434  if ( mTiledMode ) //write to file
435  {
436  QgsRasterDataProvider* partDestProvider = createPartProvider( outputExtent,
437  nCols, iterCols, iterRows,
438  iterLeft, iterTop, mOutputUrl,
439  fileIndex, nBands, destDataType, crs );
440 
441  if ( partDestProvider )
442  {
443  //write data to output file. todo: loop over the data list
444  for ( int i = 1; i <= nBands; ++i )
445  {
446  if ( destHasNoDataValueList.value( i - 1 ) )
447  {
448  partDestProvider->setNoDataValue( i, destNoDataValueList.value( i - 1 ) );
449  }
450  partDestProvider->write( destBlockList[i - 1]->bits( 0 ), i, iterCols, iterRows, 0, 0 );
451  delete destBlockList[i - 1];
452  addToVRT( partFileName( fileIndex ), i, iterCols, iterRows, iterLeft, iterTop );
453  }
454  delete partDestProvider;
455  }
456  }
457  else if ( destProvider )
458  {
459  //loop over data
460  for ( int i = 1; i <= nBands; ++i )
461  {
462  destProvider->write( destBlockList[i - 1]->bits( 0 ), i, iterCols, iterRows, iterLeft, iterTop );
463  delete destBlockList[i - 1];
464  }
465  }
466  ++fileIndex;
467  }
468 
469  QgsDebugMsgLevel( "Done", 4 );
470  return NoError;
471 }
472 
473 QgsRasterFileWriter::WriterError QgsRasterFileWriter::writeImageRaster( QgsRasterIterator* iter, int nCols, int nRows, const QgsRectangle& outputExtent,
474  const QgsCoordinateReferenceSystem& crs, QProgressDialog* progressDialog )
475 {
476  QgsDebugMsgLevel( "Entered", 4 );
477  if ( !iter )
478  {
479  return SourceProviderError;
480  }
481 
482  const QgsRasterInterface* iface = iter->input();
483  if ( !iface )
484  return SourceProviderError;
485 
486  Qgis::DataType inputDataType = iface->dataType( 1 );
487  if ( inputDataType != Qgis::ARGB32 && inputDataType != Qgis::ARGB32_Premultiplied )
488  {
489  return SourceProviderError;
490  }
491 
492  iter->setMaximumTileWidth( mMaxTileWidth );
493  iter->setMaximumTileHeight( mMaxTileHeight );
494 
495  void* redData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
496  void* greenData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
497  void* blueData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
498  void* alphaData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
499  QgsRectangle mapRect;
500  int iterLeft = 0, iterTop = 0, iterCols = 0, iterRows = 0;
501  int fileIndex = 0;
502 
503  //create destProvider for whole dataset here
504  QgsRasterDataProvider* destProvider = nullptr;
505  double pixelSize;
506  double geoTransform[6];
507  globalOutputParameters( outputExtent, nCols, nRows, geoTransform, pixelSize );
508 
509  destProvider = initOutput( nCols, nRows, crs, geoTransform, 4, Qgis::Byte );
510 
511  iter->startRasterRead( 1, nCols, nRows, outputExtent );
512 
513  int nParts = 0;
514  if ( progressDialog )
515  {
516  int nPartsX = nCols / iter->maximumTileWidth() + 1;
517  int nPartsY = nRows / iter->maximumTileHeight() + 1;
518  nParts = nPartsX * nPartsY;
519  progressDialog->setMaximum( nParts );
520  progressDialog->show();
521  progressDialog->setLabelText( QObject::tr( "Reading raster part %1 of %2" ).arg( fileIndex + 1 ).arg( nParts ) );
522  }
523 
524  QgsRasterBlock *inputBlock = nullptr;
525  while ( iter->readNextRasterPart( 1, iterCols, iterRows, &inputBlock, iterLeft, iterTop ) )
526  {
527  if ( !inputBlock )
528  {
529  continue;
530  }
531 
532  if ( progressDialog && fileIndex < ( nParts - 1 ) )
533  {
534  progressDialog->setValue( fileIndex + 1 );
535  progressDialog->setLabelText( QObject::tr( "Reading raster part %1 of %2" ).arg( fileIndex + 2 ).arg( nParts ) );
536  QCoreApplication::processEvents( QEventLoop::AllEvents, 1000 );
537  if ( progressDialog->wasCanceled() )
538  {
539  delete inputBlock;
540  break;
541  }
542  }
543 
544  //fill into red/green/blue/alpha channels
545  qgssize nPixels = static_cast< qgssize >( iterCols ) * iterRows;
546  // TODO: should be char not int? we are then copying 1 byte
547  int red = 0;
548  int green = 0;
549  int blue = 0;
550  int alpha = 255;
551  for ( qgssize i = 0; i < nPixels; ++i )
552  {
553  QRgb c = inputBlock->color( i );
554  alpha = qAlpha( c );
555  red = qRed( c );
556  green = qGreen( c );
557  blue = qBlue( c );
558 
559  if ( inputDataType == Qgis::ARGB32_Premultiplied )
560  {
561  double a = alpha / 255.;
562  QgsDebugMsgLevel( QString( "red = %1 green = %2 blue = %3 alpha = %4 p = %5 a = %6" ).arg( red ).arg( green ).arg( blue ).arg( alpha ).arg( static_cast< int >( c ), 0, 16 ).arg( a ), 5 );
563  red /= a;
564  green /= a;
565  blue /= a;
566  }
567  memcpy( reinterpret_cast< char* >( redData ) + i, &red, 1 );
568  memcpy( reinterpret_cast< char* >( greenData ) + i, &green, 1 );
569  memcpy( reinterpret_cast< char* >( blueData ) + i, &blue, 1 );
570  memcpy( reinterpret_cast< char* >( alphaData ) + i, &alpha, 1 );
571  }
572  delete inputBlock;
573 
574  //create output file
575  if ( mTiledMode )
576  {
577  //delete destProvider;
578  QgsRasterDataProvider* partDestProvider = createPartProvider( outputExtent,
579  nCols, iterCols, iterRows,
580  iterLeft, iterTop, mOutputUrl, fileIndex,
581  4, Qgis::Byte, crs );
582 
583  if ( partDestProvider )
584  {
585  //write data to output file
586  partDestProvider->write( redData, 1, iterCols, iterRows, 0, 0 );
587  partDestProvider->write( greenData, 2, iterCols, iterRows, 0, 0 );
588  partDestProvider->write( blueData, 3, iterCols, iterRows, 0, 0 );
589  partDestProvider->write( alphaData, 4, iterCols, iterRows, 0, 0 );
590 
591  addToVRT( partFileName( fileIndex ), 1, iterCols, iterRows, iterLeft, iterTop );
592  addToVRT( partFileName( fileIndex ), 2, iterCols, iterRows, iterLeft, iterTop );
593  addToVRT( partFileName( fileIndex ), 3, iterCols, iterRows, iterLeft, iterTop );
594  addToVRT( partFileName( fileIndex ), 4, iterCols, iterRows, iterLeft, iterTop );
595  delete partDestProvider;
596  }
597  }
598  else if ( destProvider )
599  {
600  destProvider->write( redData, 1, iterCols, iterRows, iterLeft, iterTop );
601  destProvider->write( greenData, 2, iterCols, iterRows, iterLeft, iterTop );
602  destProvider->write( blueData, 3, iterCols, iterRows, iterLeft, iterTop );
603  destProvider->write( alphaData, 4, iterCols, iterRows, iterLeft, iterTop );
604  }
605 
606  ++fileIndex;
607  }
608 
609  if ( destProvider )
610  delete destProvider;
611 
612  qgsFree( redData );
613  qgsFree( greenData );
614  qgsFree( blueData );
615  qgsFree( alphaData );
616 
617  if ( progressDialog )
618  {
619  progressDialog->setValue( progressDialog->maximum() );
620  }
621 
622  if ( mTiledMode )
623  {
624  QString vrtFilePath( mOutputUrl + '/' + vrtFileName() );
625  writeVRT( vrtFilePath );
626  if ( mBuildPyramidsFlag == QgsRaster::PyramidsFlagYes )
627  {
628  buildPyramids( vrtFilePath );
629  }
630  }
631  else
632  {
633  if ( mBuildPyramidsFlag == QgsRaster::PyramidsFlagYes )
634  {
635  buildPyramids( mOutputUrl );
636  }
637  }
638  return NoError;
639 }
640 
641 void QgsRasterFileWriter::addToVRT( const QString& filename, int band, int xSize, int ySize, int xOffset, int yOffset )
642 {
643  QDomElement bandElem = mVRTBands.value( band - 1 );
644 
645  QDomElement simpleSourceElem = mVRTDocument.createElement( QStringLiteral( "SimpleSource" ) );
646 
647  //SourceFilename
648  QDomElement sourceFilenameElem = mVRTDocument.createElement( QStringLiteral( "SourceFilename" ) );
649  sourceFilenameElem.setAttribute( QStringLiteral( "relativeToVRT" ), QStringLiteral( "1" ) );
650  QDomText sourceFilenameText = mVRTDocument.createTextNode( filename );
651  sourceFilenameElem.appendChild( sourceFilenameText );
652  simpleSourceElem.appendChild( sourceFilenameElem );
653 
654  //SourceBand
655  QDomElement sourceBandElem = mVRTDocument.createElement( QStringLiteral( "SourceBand" ) );
656  QDomText sourceBandText = mVRTDocument.createTextNode( QString::number( band ) );
657  sourceBandElem.appendChild( sourceBandText );
658  simpleSourceElem.appendChild( sourceBandElem );
659 
660  //SourceProperties
661  QDomElement sourcePropertiesElem = mVRTDocument.createElement( QStringLiteral( "SourceProperties" ) );
662  sourcePropertiesElem.setAttribute( QStringLiteral( "RasterXSize" ), xSize );
663  sourcePropertiesElem.setAttribute( QStringLiteral( "RasterYSize" ), ySize );
664  sourcePropertiesElem.setAttribute( QStringLiteral( "BlockXSize" ), xSize );
665  sourcePropertiesElem.setAttribute( QStringLiteral( "BlockYSize" ), ySize );
666  sourcePropertiesElem.setAttribute( QStringLiteral( "DataType" ), QStringLiteral( "Byte" ) );
667  simpleSourceElem.appendChild( sourcePropertiesElem );
668 
669  //SrcRect
670  QDomElement srcRectElem = mVRTDocument.createElement( QStringLiteral( "SrcRect" ) );
671  srcRectElem.setAttribute( QStringLiteral( "xOff" ), QStringLiteral( "0" ) );
672  srcRectElem.setAttribute( QStringLiteral( "yOff" ), QStringLiteral( "0" ) );
673  srcRectElem.setAttribute( QStringLiteral( "xSize" ), xSize );
674  srcRectElem.setAttribute( QStringLiteral( "ySize" ), ySize );
675  simpleSourceElem.appendChild( srcRectElem );
676 
677  //DstRect
678  QDomElement dstRectElem = mVRTDocument.createElement( QStringLiteral( "DstRect" ) );
679  dstRectElem.setAttribute( QStringLiteral( "xOff" ), xOffset );
680  dstRectElem.setAttribute( QStringLiteral( "yOff" ), yOffset );
681  dstRectElem.setAttribute( QStringLiteral( "xSize" ), xSize );
682  dstRectElem.setAttribute( QStringLiteral( "ySize" ), ySize );
683  simpleSourceElem.appendChild( dstRectElem );
684 
685  bandElem.appendChild( simpleSourceElem );
686 }
687 
688 #if 0
689 void QgsRasterFileWriter::buildPyramids( const QString& filename )
690 {
691  GDALDatasetH dataSet;
692  GDALAllRegister();
693  dataSet = GDALOpen( filename.toLocal8Bit().data(), GA_Update );
694  if ( !dataSet )
695  {
696  return;
697  }
698 
699  //2,4,8,16,32,64
700  int overviewList[6];
701  overviewList[0] = 2;
702  overviewList[1] = 4;
703  overviewList[2] = 8;
704  overviewList[3] = 16;
705  overviewList[4] = 32;
706  overviewList[5] = 64;
707 
708 #if 0
709  if ( mProgressDialog )
710  {
711  mProgressDialog->setLabelText( QObject::tr( "Building Pyramids..." ) );
712  mProgressDialog->setValue( 0 );
713  mProgressDialog->setWindowModality( Qt::WindowModal );
714  mProgressDialog->show();
715  }
716 #endif
717  GDALBuildOverviews( dataSet, "AVERAGE", 6, overviewList, 0, 0, /*pyramidsProgress*/ 0, /*mProgressDialog*/ 0 );
718 }
719 #endif
720 
721 void QgsRasterFileWriter::buildPyramids( const QString& filename )
722 {
723  QgsDebugMsgLevel( "filename = " + filename, 4 );
724  // open new dataProvider so we can build pyramids with it
725  QgsRasterDataProvider* destProvider = dynamic_cast< QgsRasterDataProvider* >( QgsProviderRegistry::instance()->provider( mOutputProviderKey, filename ) );
726  if ( !destProvider )
727  {
728  return;
729  }
730 
731  // TODO progress report
732  // TODO test mTiledMode - not tested b/c segfault at line # 289
733  // connect( provider, SIGNAL( progressUpdate( int ) ), mPyramidProgress, SLOT( setValue( int ) ) );
734  QList< QgsRasterPyramid> myPyramidList;
735  if ( ! mPyramidsList.isEmpty() )
736  myPyramidList = destProvider->buildPyramidList( mPyramidsList );
737  for ( int myCounterInt = 0; myCounterInt < myPyramidList.count(); myCounterInt++ )
738  {
739  myPyramidList[myCounterInt].build = true;
740  }
741 
742  QgsDebugMsgLevel( QString( "building pyramids : %1 pyramids, %2 resampling, %3 format, %4 options" ).arg( myPyramidList.count() ).arg( mPyramidsResampling ).arg( mPyramidsFormat ).arg( mPyramidsConfigOptions.count() ), 4 );
743  // QApplication::setOverrideCursor( Qt::WaitCursor );
744  QString res = destProvider->buildPyramids( myPyramidList, mPyramidsResampling,
745  mPyramidsFormat, mPyramidsConfigOptions );
746  // QApplication::restoreOverrideCursor();
747 
748  // TODO put this in provider or elsewhere
749  if ( !res.isNull() )
750  {
751  QString title, message;
752  if ( res == QLatin1String( "ERROR_WRITE_ACCESS" ) )
753  {
754  title = QObject::tr( "Building pyramids failed - write access denied" );
755  message = QObject::tr( "Write access denied. Adjust the file permissions and try again." );
756  }
757  else if ( res == QLatin1String( "ERROR_WRITE_FORMAT" ) )
758  {
759  title = QObject::tr( "Building pyramids failed." );
760  message = QObject::tr( "The file was not writable. Some formats do not "
761  "support pyramid overviews. Consult the GDAL documentation if in doubt." );
762  }
763  else if ( res == QLatin1String( "FAILED_NOT_SUPPORTED" ) )
764  {
765  title = QObject::tr( "Building pyramids failed." );
766  message = QObject::tr( "Building pyramid overviews is not supported on this type of raster." );
767  }
768  else if ( res == QLatin1String( "ERROR_JPEG_COMPRESSION" ) )
769  {
770  title = QObject::tr( "Building pyramids failed." );
771  message = QObject::tr( "Building internal pyramid overviews is not supported on raster layers with JPEG compression and your current libtiff library." );
772  }
773  else if ( res == QLatin1String( "ERROR_VIRTUAL" ) )
774  {
775  title = QObject::tr( "Building pyramids failed." );
776  message = QObject::tr( "Building pyramid overviews is not supported on this type of raster." );
777  }
778  QMessageBox::warning( nullptr, title, message );
779  QgsDebugMsgLevel( res + " - " + message, 4 );
780  }
781  delete destProvider;
782 }
783 
784 #if 0
785 int QgsRasterFileWriter::pyramidsProgress( double dfComplete, const char *pszMessage, void* pData )
786 {
787  Q_UNUSED( pszMessage );
788  GDALTermProgress( dfComplete, 0, 0 );
789  QProgressDialog* progressDialog = static_cast<QProgressDialog*>( pData );
790  if ( pData && progressDialog->wasCanceled() )
791  {
792  return 0;
793  }
794 
795  if ( pData )
796  {
797  progressDialog->setRange( 0, 100 );
798  progressDialog->setValue( dfComplete * 100 );
799  }
800  return 1;
801 }
802 #endif
803 
804 void QgsRasterFileWriter::createVRT( int xSize, int ySize, const QgsCoordinateReferenceSystem& crs, double* geoTransform, Qgis::DataType type, const QList<bool>& destHasNoDataValueList, const QList<double>& destNoDataValueList )
805 {
806  mVRTDocument.clear();
807  QDomElement VRTDatasetElem = mVRTDocument.createElement( QStringLiteral( "VRTDataset" ) );
808 
809  //xsize / ysize
810  VRTDatasetElem.setAttribute( QStringLiteral( "rasterXSize" ), xSize );
811  VRTDatasetElem.setAttribute( QStringLiteral( "rasterYSize" ), ySize );
812  mVRTDocument.appendChild( VRTDatasetElem );
813 
814  //CRS
815  QDomElement SRSElem = mVRTDocument.createElement( QStringLiteral( "SRS" ) );
816  QDomText crsText = mVRTDocument.createTextNode( crs.toWkt() );
817  SRSElem.appendChild( crsText );
818  VRTDatasetElem.appendChild( SRSElem );
819 
820  //geotransform
821  if ( geoTransform )
822  {
823  QDomElement geoTransformElem = mVRTDocument.createElement( QStringLiteral( "GeoTransform" ) );
824  QString geoTransformString = QString::number( geoTransform[0] ) + ", " + QString::number( geoTransform[1] ) + ", " + QString::number( geoTransform[2] ) +
825  ", " + QString::number( geoTransform[3] ) + ", " + QString::number( geoTransform[4] ) + ", " + QString::number( geoTransform[5] );
826  QDomText geoTransformText = mVRTDocument.createTextNode( geoTransformString );
827  geoTransformElem.appendChild( geoTransformText );
828  VRTDatasetElem.appendChild( geoTransformElem );
829  }
830 
831  int nBands;
832  if ( mMode == Raw )
833  {
834  nBands = mInput->bandCount();
835  }
836  else
837  {
838  nBands = 4;
839  }
840 
841  QStringList colorInterp;
842  colorInterp << QStringLiteral( "Red" ) << QStringLiteral( "Green" ) << QStringLiteral( "Blue" ) << QStringLiteral( "Alpha" );
843 
844  QMap<Qgis::DataType, QString> dataTypes;
845  dataTypes.insert( Qgis::Byte, QStringLiteral( "Byte" ) );
846  dataTypes.insert( Qgis::UInt16, QStringLiteral( "UInt16" ) );
847  dataTypes.insert( Qgis::Int16, QStringLiteral( "Int16" ) );
848  dataTypes.insert( Qgis::UInt32, QStringLiteral( "Int32" ) );
849  dataTypes.insert( Qgis::Float32, QStringLiteral( "Float32" ) );
850  dataTypes.insert( Qgis::Float64, QStringLiteral( "Float64" ) );
851  dataTypes.insert( Qgis::CInt16, QStringLiteral( "CInt16" ) );
852  dataTypes.insert( Qgis::CInt32, QStringLiteral( "CInt32" ) );
853  dataTypes.insert( Qgis::CFloat32, QStringLiteral( "CFloat32" ) );
854  dataTypes.insert( Qgis::CFloat64, QStringLiteral( "CFloat64" ) );
855 
856  for ( int i = 1; i <= nBands; i++ )
857  {
858  QDomElement VRTBand = mVRTDocument.createElement( QStringLiteral( "VRTRasterBand" ) );
859 
860  VRTBand.setAttribute( QStringLiteral( "band" ), QString::number( i ) );
861  QString dataType = dataTypes.value( type );
862  VRTBand.setAttribute( QStringLiteral( "dataType" ), dataType );
863 
864  if ( mMode == Image )
865  {
866  VRTBand.setAttribute( QStringLiteral( "dataType" ), QStringLiteral( "Byte" ) );
867  QDomElement colorInterpElement = mVRTDocument.createElement( QStringLiteral( "ColorInterp" ) );
868  QDomText interpText = mVRTDocument.createTextNode( colorInterp.value( i - 1 ) );
869  colorInterpElement.appendChild( interpText );
870  VRTBand.appendChild( colorInterpElement );
871  }
872 
873  if ( !destHasNoDataValueList.isEmpty() && destHasNoDataValueList.value( i - 1 ) )
874  {
875  VRTBand.setAttribute( QStringLiteral( "NoDataValue" ), QString::number( destNoDataValueList.value( i - 1 ) ) );
876  }
877 
878  mVRTBands.append( VRTBand );
879  VRTDatasetElem.appendChild( VRTBand );
880  }
881 }
882 
883 bool QgsRasterFileWriter::writeVRT( const QString& file )
884 {
885  QFile outputFile( file );
886  if ( ! outputFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
887  {
888  return false;
889  }
890 
891  QTextStream outStream( &outputFile );
892  mVRTDocument.save( outStream, 2 );
893  return true;
894 }
895 
896 QgsRasterDataProvider* QgsRasterFileWriter::createPartProvider( const QgsRectangle& extent, int nCols, int iterCols,
897  int iterRows, int iterLeft, int iterTop, const QString& outputUrl, int fileIndex, int nBands, Qgis::DataType type,
898  const QgsCoordinateReferenceSystem& crs )
899 {
900  double mup = extent.width() / nCols;
901  double mapLeft = extent.xMinimum() + iterLeft * mup;
902  double mapRight = mapLeft + mup * iterCols;
903  double mapTop = extent.yMaximum() - iterTop * mup;
904  double mapBottom = mapTop - iterRows * mup;
905  QgsRectangle mapRect( mapLeft, mapBottom, mapRight, mapTop );
906 
907  QString outputFile = outputUrl + '/' + partFileName( fileIndex );
908 
909  //geotransform
910  double geoTransform[6];
911  geoTransform[0] = mapRect.xMinimum();
912  geoTransform[1] = mup;
913  geoTransform[2] = 0.0;
914  geoTransform[3] = mapRect.yMaximum();
915  geoTransform[4] = 0.0;
916  geoTransform[5] = -mup;
917 
918  // perhaps we need a separate createOptions for tiles ?
919 
920  QgsRasterDataProvider* destProvider = QgsRasterDataProvider::create( mOutputProviderKey, outputFile, mOutputFormat, nBands, type, iterCols, iterRows, geoTransform, crs, mCreateOptions );
921 
922  // TODO: return provider and report error
923  return destProvider;
924 }
925 
926 QgsRasterDataProvider* QgsRasterFileWriter::initOutput( int nCols, int nRows, const QgsCoordinateReferenceSystem& crs,
927  double* geoTransform, int nBands, Qgis::DataType type,
928  const QList<bool>& destHasNoDataValueList, const QList<double>& destNoDataValueList )
929 {
930  if ( mTiledMode )
931  {
932  createVRT( nCols, nRows, crs, geoTransform, type, destHasNoDataValueList, destNoDataValueList );
933  return nullptr;
934  }
935  else
936  {
937 #if 0
938  // TODO enable "use existing", has no effect for now, because using Create() in gdal provider
939  // should this belong in provider? should also test that source provider is gdal
940  if ( mBuildPyramidsFlag == -4 && mOutputProviderKey == "gdal" && mOutputFormat.toLower() == "gtiff" )
941  mCreateOptions << "COPY_SRC_OVERVIEWS=YES";
942 #endif
943 
944  QgsRasterDataProvider* destProvider = QgsRasterDataProvider::create( mOutputProviderKey, mOutputUrl, mOutputFormat, nBands, type, nCols, nRows, geoTransform, crs, mCreateOptions );
945 
946  if ( !destProvider )
947  {
948  QgsDebugMsg( "No provider created" );
949  }
950 
951  return destProvider;
952  }
953 }
954 
955 void QgsRasterFileWriter::globalOutputParameters( const QgsRectangle& extent, int nCols, int& nRows,
956  double* geoTransform, double& pixelSize )
957 {
958  pixelSize = extent.width() / nCols;
959 
960  //calculate nRows automatically for providers without exact resolution
961  if ( nRows < 0 )
962  {
963  nRows = static_cast< double >( nCols ) / extent.width() * extent.height() + 0.5;
964  }
965  geoTransform[0] = extent.xMinimum();
966  geoTransform[1] = pixelSize;
967  geoTransform[2] = 0.0;
968  geoTransform[3] = extent.yMaximum();
969  geoTransform[4] = 0.0;
970  geoTransform[5] = -( extent.height() / nRows );
971 }
972 
973 QString QgsRasterFileWriter::partFileName( int fileIndex )
974 {
975  // .tif for now
976  QFileInfo outputInfo( mOutputUrl );
977  return QStringLiteral( "%1.%2.tif" ).arg( outputInfo.fileName() ).arg( fileIndex );
978 }
979 
980 QString QgsRasterFileWriter::vrtFileName()
981 {
982  QFileInfo outputInfo( mOutputUrl );
983  return QStringLiteral( "%1.vrt" ).arg( outputInfo.fileName() );
984 }
virtual QgsRectangle extent() const override=0
Returns the extent of the layer.
virtual int bandCount() const =0
Get number of bands.
bool contains(const QgsRectangle &rect) const
return true when rectangle contains other rectangle
static QgsProviderRegistry * instance(const QString &pluginPath=QString::null)
Means of accessing canonical single instance.
A rectangle specified with double values.
Definition: qgsrectangle.h:35
QgsRasterFileWriter(const QString &outputUrl)
Base class for processing modules.
Definition: qgsrasterpipe.h:43
void * qgsMalloc(size_t size)
Allocates size bytes and returns a pointer to the allocated memory.
Definition: qgis.cpp:105
Iterator for sequentially processing raster cells.
static Qgis::DataType typeWithNoDataValue(Qgis::DataType dataType, double *noDataValue)
For given data type returns wider type and sets no data value.
static bool typeIsColor(Qgis::DataType type)
Returns true if data type is color.
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
bool setValue(int row, int column, double value)
Set value on position.
QgsRasterProjector * projector() const
Raster pipe that deals with null values.
Thirty two bit unsigned integer (quint32)
Definition: qgis.h:66
DataType
Raster data types.
Definition: qgis.h:60
QgsRasterInterface * last() const
Definition: qgsrasterpipe.h:88
double maximumValue
The maximum cell value in the raster band.
virtual bool setNoDataValue(int bandNo, double noDataValue)
Set no data value on created dataset.
virtual Qgis::DataType dataType(int bandNo) const =0
Returns data type for the band specified by number.
Thirty two bit floating point (float)
Definition: qgis.h:68
WriterError writeRaster(const QgsRasterPipe *pipe, int nCols, int nRows, const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &crs, QProgressDialog *p=nullptr)
Write raster file.
Sixteen bit signed integer (qint16)
Definition: qgis.h:65
Complex Int16.
Definition: qgis.h:70
virtual double sourceNoDataValue(int bandNo) const
Value representing no data value.
Sixty four bit floating point (double)
Definition: qgis.h:69
virtual Qgis::DataType sourceDataType(int bandNo) const override=0
Returns source data type for the band specified by number, source data type may be shorter than dataT...
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
Definition: qgis.h:75
virtual QgsRasterBandStats bandStatistics(int theBandNo, int theStats=QgsRasterBandStats::All, const QgsRectangle &theExtent=QgsRectangle(), int theSampleSize=0)
Get band statistics.
The RasterBandStats struct is a container for statistics about a single raster band.
virtual Qgis::DataType dataType(int bandNo) const override=0
Returns data type for the band specified by number.
Raster data container.
QgsCoordinateReferenceSystem sourceCrs() const
Get source CRS.
QgsDataProvider * provider(const QString &providerKey, const QString &dataSource)
Create an instance of the provider.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
virtual QList< QgsRasterPyramid > buildPyramidList(QList< int > overviewList=QList< int >())
Accessor for ths raster layers pyramid list.
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:211
Complex Float32.
Definition: qgis.h:72
QRgb color(int row, int column) const
Read a single color.
Complex Int32.
Definition: qgis.h:71
static int typeSize(int dataType)
Sixteen bit unsigned integer (quint16)
Definition: qgis.h:64
QgsRasterRangeList noData(int bandNo) const
bool readNextRasterPart(int bandNumber, int &nCols, int &nRows, QgsRasterBlock **block, int &topLeftCol, int &topLeftRow)
Fetches next part of raster data, caller takes ownership of the block and caller should delete the bl...
void * GDALDatasetH
Base class for processing filters like renderers, reprojector, resampler etc.
virtual bool sourceHasNoDataValue(int bandNo) const
Return true if source band has no data value.
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
Definition: qgis.h:336
void setOutputNoDataValue(int bandNo, double noData)
Set output no data value.
QgsRasterProjector implements approximate projection support for it calculates grid of points in sour...
void setMaximumTileWidth(int w)
Raster namespace.
Definition: qgsraster.h:28
QgsCoordinateReferenceSystem destinationCrs() const
Get destination CRS.
void setMaximumTileHeight(int h)
const QgsRasterInterface * input() const
int maximumTileHeight() const
static double maximumValuePossible(Qgis::DataType)
Helper function that returns the maximum possible value for a GDAL data type.
virtual bool remove()
Remove dataset.
This class represents a coordinate reference system (CRS).
QString toWkt() const
Returns a WKT representation of this CRS.
Class for doing transforms between two map coordinate systems.
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:196
double minimumValue
The minimum cell value in the raster band.
QgsRasterNuller * nuller() const
virtual QString buildPyramids(const QList< QgsRasterPyramid > &thePyramidList, const QString &theResamplingMethod="NEAREST", QgsRaster::RasterPyramidsFormat theFormat=QgsRaster::PyramidsGTiff, const QStringList &theConfigOptions=QStringList())
Create pyramid overviews.
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:201
virtual bool write(void *data, int band, int width, int height, int xOffset, int yOffset)
Writes into the provider datasource.
Complex Float64.
Definition: qgis.h:73
void qgsFree(void *ptr)
Frees the memory space pointed to by ptr.
Definition: qgis.cpp:135
void startRasterRead(int bandNumber, int nCols, int nRows, const QgsRectangle &extent, QgsRasterBlockFeedback *feedback=nullptr)
Start reading of raster band.
Eight bit unsigned integer (quint8)
Definition: qgis.h:63
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32.
Definition: qgis.h:74
virtual const QgsRasterInterface * sourceInput() const
Get source / raw input, the first in pipe, usually provider.
int maximumTileWidth() const
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:216
Base class for raster data providers.
static QgsRasterDataProvider * create(const QString &providerKey, const QString &uri, const QString &format, int nBands, Qgis::DataType type, int width, int height, double *geoTransform, const QgsCoordinateReferenceSystem &crs, const QStringList &createOptions=QStringList())
Creates a new dataset with mDataSourceURI.