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