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