QGIS API Documentation  2.15.0-Master (af20121)
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 ) // 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  for ( ;; )
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  partDestProvider->setNoDataValue( i, destNoDataValueList.value( i - 1 ) );
440  partDestProvider->write( destBlockList[i - 1]->bits( 0 ), i, iterCols, iterRows, 0, 0 );
441  delete destBlockList[i - 1];
442  addToVRT( partFileName( fileIndex ), i, iterCols, iterRows, iterLeft, iterTop );
443  }
444  delete partDestProvider;
445  }
446  }
447  else if ( destProvider )
448  {
449  //loop over data
450  for ( int i = 1; i <= nBands; ++i )
451  {
452  destProvider->write( destBlockList[i - 1]->bits( 0 ), i, iterCols, iterRows, iterLeft, iterTop );
453  delete destBlockList[i - 1];
454  }
455  }
456  ++fileIndex;
457  }
458 
459  QgsDebugMsgLevel( "Done", 4 );
460  return NoError;
461 }
462 
463 QgsRasterFileWriter::WriterError QgsRasterFileWriter::writeImageRaster( QgsRasterIterator* iter, int nCols, int nRows, const QgsRectangle& outputExtent,
464  const QgsCoordinateReferenceSystem& crs, QProgressDialog* progressDialog )
465 {
466  QgsDebugMsgLevel( "Entered", 4 );
467  if ( !iter )
468  {
469  return SourceProviderError;
470  }
471 
472  const QgsRasterInterface* iface = iter->input();
473  if ( !iface )
474  return SourceProviderError;
475 
476  QGis::DataType inputDataType = iface->dataType( 1 );
477  if ( inputDataType != QGis::ARGB32 && inputDataType != QGis::ARGB32_Premultiplied )
478  {
479  return SourceProviderError;
480  }
481 
482  iter->setMaximumTileWidth( mMaxTileWidth );
483  iter->setMaximumTileHeight( mMaxTileHeight );
484 
485  void* redData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
486  void* greenData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
487  void* blueData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
488  void* alphaData = qgsMalloc( mMaxTileWidth * mMaxTileHeight );
489  QgsRectangle mapRect;
490  int iterLeft = 0, iterTop = 0, iterCols = 0, iterRows = 0;
491  int fileIndex = 0;
492 
493  //create destProvider for whole dataset here
494  QgsRasterDataProvider* destProvider = nullptr;
495  double pixelSize;
496  double geoTransform[6];
497  globalOutputParameters( outputExtent, nCols, nRows, geoTransform, pixelSize );
498 
499  destProvider = initOutput( nCols, nRows, crs, geoTransform, 4, QGis::Byte );
500 
501  iter->startRasterRead( 1, nCols, nRows, outputExtent );
502 
503  int nParts = 0;
504  if ( progressDialog )
505  {
506  int nPartsX = nCols / iter->maximumTileWidth() + 1;
507  int nPartsY = nRows / iter->maximumTileHeight() + 1;
508  nParts = nPartsX * nPartsY;
509  progressDialog->setMaximum( nParts );
510  progressDialog->show();
511  progressDialog->setLabelText( QObject::tr( "Reading raster part %1 of %2" ).arg( fileIndex + 1 ).arg( nParts ) );
512  }
513 
514  QgsRasterBlock *inputBlock = nullptr;
515  while ( iter->readNextRasterPart( 1, iterCols, iterRows, &inputBlock, iterLeft, iterTop ) )
516  {
517  if ( !inputBlock )
518  {
519  continue;
520  }
521 
522  if ( progressDialog && fileIndex < ( nParts - 1 ) )
523  {
524  progressDialog->setValue( fileIndex + 1 );
525  progressDialog->setLabelText( QObject::tr( "Reading raster part %1 of %2" ).arg( fileIndex + 2 ).arg( nParts ) );
526  QCoreApplication::processEvents( QEventLoop::AllEvents, 1000 );
527  if ( progressDialog->wasCanceled() )
528  {
529  delete inputBlock;
530  break;
531  }
532  }
533 
534  //fill into red/green/blue/alpha channels
535  qgssize nPixels = static_cast< qgssize >( iterCols ) * iterRows;
536  // TODO: should be char not int? we are then copying 1 byte
537  int red = 0;
538  int green = 0;
539  int blue = 0;
540  int alpha = 255;
541  for ( qgssize i = 0; i < nPixels; ++i )
542  {
543  QRgb c = inputBlock->color( i );
544  alpha = qAlpha( c );
545  red = qRed( c );
546  green = qGreen( c );
547  blue = qBlue( c );
548 
549  if ( inputDataType == QGis::ARGB32_Premultiplied )
550  {
551  double a = alpha / 255.;
552  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 );
553  red /= a;
554  green /= a;
555  blue /= a;
556  }
557  memcpy( reinterpret_cast< char* >( redData ) + i, &red, 1 );
558  memcpy( reinterpret_cast< char* >( greenData ) + i, &green, 1 );
559  memcpy( reinterpret_cast< char* >( blueData ) + i, &blue, 1 );
560  memcpy( reinterpret_cast< char* >( alphaData ) + i, &alpha, 1 );
561  }
562  delete inputBlock;
563 
564  //create output file
565  if ( mTiledMode )
566  {
567  //delete destProvider;
568  QgsRasterDataProvider* partDestProvider = createPartProvider( outputExtent,
569  nCols, iterCols, iterRows,
570  iterLeft, iterTop, mOutputUrl, fileIndex,
571  4, QGis::Byte, crs );
572 
573  if ( partDestProvider )
574  {
575  //write data to output file
576  partDestProvider->write( redData, 1, iterCols, iterRows, 0, 0 );
577  partDestProvider->write( greenData, 2, iterCols, iterRows, 0, 0 );
578  partDestProvider->write( blueData, 3, iterCols, iterRows, 0, 0 );
579  partDestProvider->write( alphaData, 4, iterCols, iterRows, 0, 0 );
580 
581  addToVRT( partFileName( fileIndex ), 1, iterCols, iterRows, iterLeft, iterTop );
582  addToVRT( partFileName( fileIndex ), 2, iterCols, iterRows, iterLeft, iterTop );
583  addToVRT( partFileName( fileIndex ), 3, iterCols, iterRows, iterLeft, iterTop );
584  addToVRT( partFileName( fileIndex ), 4, iterCols, iterRows, iterLeft, iterTop );
585  delete partDestProvider;
586  }
587  }
588  else if ( destProvider )
589  {
590  destProvider->write( redData, 1, iterCols, iterRows, iterLeft, iterTop );
591  destProvider->write( greenData, 2, iterCols, iterRows, iterLeft, iterTop );
592  destProvider->write( blueData, 3, iterCols, iterRows, iterLeft, iterTop );
593  destProvider->write( alphaData, 4, iterCols, iterRows, iterLeft, iterTop );
594  }
595 
596  ++fileIndex;
597  }
598 
599  if ( destProvider )
600  delete destProvider;
601 
602  qgsFree( redData );
603  qgsFree( greenData );
604  qgsFree( blueData );
605  qgsFree( alphaData );
606 
607  if ( progressDialog )
608  {
609  progressDialog->setValue( progressDialog->maximum() );
610  }
611 
612  if ( mTiledMode )
613  {
614  QString vrtFilePath( mOutputUrl + '/' + vrtFileName() );
615  writeVRT( vrtFilePath );
616  if ( mBuildPyramidsFlag == QgsRaster::PyramidsFlagYes )
617  {
618  buildPyramids( vrtFilePath );
619  }
620  }
621  else
622  {
623  if ( mBuildPyramidsFlag == QgsRaster::PyramidsFlagYes )
624  {
625  buildPyramids( mOutputUrl );
626  }
627  }
628  return NoError;
629 }
630 
631 void QgsRasterFileWriter::addToVRT( const QString& filename, int band, int xSize, int ySize, int xOffset, int yOffset )
632 {
633  QDomElement bandElem = mVRTBands.value( band - 1 );
634 
635  QDomElement simpleSourceElem = mVRTDocument.createElement( "SimpleSource" );
636 
637  //SourceFilename
638  QDomElement sourceFilenameElem = mVRTDocument.createElement( "SourceFilename" );
639  sourceFilenameElem.setAttribute( "relativeToVRT", "1" );
640  QDomText sourceFilenameText = mVRTDocument.createTextNode( filename );
641  sourceFilenameElem.appendChild( sourceFilenameText );
642  simpleSourceElem.appendChild( sourceFilenameElem );
643 
644  //SourceBand
645  QDomElement sourceBandElem = mVRTDocument.createElement( "SourceBand" );
646  QDomText sourceBandText = mVRTDocument.createTextNode( QString::number( band ) );
647  sourceBandElem.appendChild( sourceBandText );
648  simpleSourceElem.appendChild( sourceBandElem );
649 
650  //SourceProperties
651  QDomElement sourcePropertiesElem = mVRTDocument.createElement( "SourceProperties" );
652  sourcePropertiesElem.setAttribute( "RasterXSize", xSize );
653  sourcePropertiesElem.setAttribute( "RasterYSize", ySize );
654  sourcePropertiesElem.setAttribute( "BlockXSize", xSize );
655  sourcePropertiesElem.setAttribute( "BlockYSize", ySize );
656  sourcePropertiesElem.setAttribute( "DataType", "Byte" );
657  simpleSourceElem.appendChild( sourcePropertiesElem );
658 
659  //SrcRect
660  QDomElement srcRectElem = mVRTDocument.createElement( "SrcRect" );
661  srcRectElem.setAttribute( "xOff", "0" );
662  srcRectElem.setAttribute( "yOff", "0" );
663  srcRectElem.setAttribute( "xSize", xSize );
664  srcRectElem.setAttribute( "ySize", ySize );
665  simpleSourceElem.appendChild( srcRectElem );
666 
667  //DstRect
668  QDomElement dstRectElem = mVRTDocument.createElement( "DstRect" );
669  dstRectElem.setAttribute( "xOff", xOffset );
670  dstRectElem.setAttribute( "yOff", yOffset );
671  dstRectElem.setAttribute( "xSize", xSize );
672  dstRectElem.setAttribute( "ySize", ySize );
673  simpleSourceElem.appendChild( dstRectElem );
674 
675  bandElem.appendChild( simpleSourceElem );
676 }
677 
678 #if 0
679 void QgsRasterFileWriter::buildPyramids( const QString& filename )
680 {
681  GDALDatasetH dataSet;
682  GDALAllRegister();
683  dataSet = GDALOpen( filename.toLocal8Bit().data(), GA_Update );
684  if ( !dataSet )
685  {
686  return;
687  }
688 
689  //2,4,8,16,32,64
690  int overviewList[6];
691  overviewList[0] = 2;
692  overviewList[1] = 4;
693  overviewList[2] = 8;
694  overviewList[3] = 16;
695  overviewList[4] = 32;
696  overviewList[5] = 64;
697 
698 #if 0
699  if ( mProgressDialog )
700  {
701  mProgressDialog->setLabelText( QObject::tr( "Building Pyramids..." ) );
702  mProgressDialog->setValue( 0 );
703  mProgressDialog->setWindowModality( Qt::WindowModal );
704  mProgressDialog->show();
705  }
706 #endif
707  GDALBuildOverviews( dataSet, "AVERAGE", 6, overviewList, 0, 0, /*pyramidsProgress*/ 0, /*mProgressDialog*/ 0 );
708 }
709 #endif
710 
711 void QgsRasterFileWriter::buildPyramids( const QString& filename )
712 {
713  QgsDebugMsgLevel( "filename = " + filename, 4 );
714  // open new dataProvider so we can build pyramids with it
715  QgsRasterDataProvider* destProvider = dynamic_cast< QgsRasterDataProvider* >( QgsProviderRegistry::instance()->provider( mOutputProviderKey, filename ) );
716  if ( !destProvider )
717  {
718  return;
719  }
720 
721  // TODO progress report
722  // TODO test mTiledMode - not tested b/c segfault at line # 289
723  // connect( provider, SIGNAL( progressUpdate( int ) ), mPyramidProgress, SLOT( setValue( int ) ) );
724  QList< QgsRasterPyramid> myPyramidList;
725  if ( ! mPyramidsList.isEmpty() )
726  myPyramidList = destProvider->buildPyramidList( mPyramidsList );
727  for ( int myCounterInt = 0; myCounterInt < myPyramidList.count(); myCounterInt++ )
728  {
729  myPyramidList[myCounterInt].build = true;
730  }
731 
732  QgsDebugMsgLevel( QString( "building pyramids : %1 pyramids, %2 resampling, %3 format, %4 options" ).arg( myPyramidList.count() ).arg( mPyramidsResampling ).arg( mPyramidsFormat ).arg( mPyramidsConfigOptions.count() ), 4 );
733  // QApplication::setOverrideCursor( Qt::WaitCursor );
734  QString res = destProvider->buildPyramids( myPyramidList, mPyramidsResampling,
735  mPyramidsFormat, mPyramidsConfigOptions );
736  // QApplication::restoreOverrideCursor();
737 
738  // TODO put this in provider or elsewhere
739  if ( !res.isNull() )
740  {
741  QString title, message;
742  if ( res == "ERROR_WRITE_ACCESS" )
743  {
744  title = QObject::tr( "Building pyramids failed - write access denied" );
745  message = QObject::tr( "Write access denied. Adjust the file permissions and try again." );
746  }
747  else if ( res == "ERROR_WRITE_FORMAT" )
748  {
749  title = QObject::tr( "Building pyramids failed." );
750  message = QObject::tr( "The file was not writable. Some formats do not "
751  "support pyramid overviews. Consult the GDAL documentation if in doubt." );
752  }
753  else if ( res == "FAILED_NOT_SUPPORTED" )
754  {
755  title = QObject::tr( "Building pyramids failed." );
756  message = QObject::tr( "Building pyramid overviews is not supported on this type of raster." );
757  }
758  else if ( res == "ERROR_JPEG_COMPRESSION" )
759  {
760  title = QObject::tr( "Building pyramids failed." );
761  message = QObject::tr( "Building internal pyramid overviews is not supported on raster layers with JPEG compression and your current libtiff library." );
762  }
763  else if ( res == "ERROR_VIRTUAL" )
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  QMessageBox::warning( nullptr, title, message );
769  QgsDebugMsgLevel( res + " - " + message, 4 );
770  }
771  delete destProvider;
772 }
773 
774 #if 0
775 int QgsRasterFileWriter::pyramidsProgress( double dfComplete, const char *pszMessage, void* pData )
776 {
777  Q_UNUSED( pszMessage );
778  GDALTermProgress( dfComplete, 0, 0 );
779  QProgressDialog* progressDialog = static_cast<QProgressDialog*>( pData );
780  if ( pData && progressDialog->wasCanceled() )
781  {
782  return 0;
783  }
784 
785  if ( pData )
786  {
787  progressDialog->setRange( 0, 100 );
788  progressDialog->setValue( dfComplete * 100 );
789  }
790  return 1;
791 }
792 #endif
793 
794 void QgsRasterFileWriter::createVRT( int xSize, int ySize, const QgsCoordinateReferenceSystem& crs, double* geoTransform, QGis::DataType type, const QList<bool>& destHasNoDataValueList, const QList<double>& destNoDataValueList )
795 {
796  mVRTDocument.clear();
797  QDomElement VRTDatasetElem = mVRTDocument.createElement( "VRTDataset" );
798 
799  //xsize / ysize
800  VRTDatasetElem.setAttribute( "rasterXSize", xSize );
801  VRTDatasetElem.setAttribute( "rasterYSize", ySize );
802  mVRTDocument.appendChild( VRTDatasetElem );
803 
804  //CRS
805  QDomElement SRSElem = mVRTDocument.createElement( "SRS" );
806  QDomText crsText = mVRTDocument.createTextNode( crs.toWkt() );
807  SRSElem.appendChild( crsText );
808  VRTDatasetElem.appendChild( SRSElem );
809 
810  //geotransform
811  if ( geoTransform )
812  {
813  QDomElement geoTransformElem = mVRTDocument.createElement( "GeoTransform" );
814  QString geoTransformString = QString::number( geoTransform[0] ) + ", " + QString::number( geoTransform[1] ) + ", " + QString::number( geoTransform[2] ) +
815  ", " + QString::number( geoTransform[3] ) + ", " + QString::number( geoTransform[4] ) + ", " + QString::number( geoTransform[5] );
816  QDomText geoTransformText = mVRTDocument.createTextNode( geoTransformString );
817  geoTransformElem.appendChild( geoTransformText );
818  VRTDatasetElem.appendChild( geoTransformElem );
819  }
820 
821  int nBands;
822  if ( mMode == Raw )
823  {
824  nBands = mInput->bandCount();
825  }
826  else
827  {
828  nBands = 4;
829  }
830 
831  QStringList colorInterp;
832  colorInterp << "Red" << "Green" << "Blue" << "Alpha";
833 
835  dataTypes.insert( QGis::Byte, "Byte" );
836  dataTypes.insert( QGis::UInt16, "UInt16" );
837  dataTypes.insert( QGis::Int16, "Int16" );
838  dataTypes.insert( QGis::UInt32, "Int32" );
839  dataTypes.insert( QGis::Float32, "Float32" );
840  dataTypes.insert( QGis::Float64, "Float64" );
841  dataTypes.insert( QGis::CInt16, "CInt16" );
842  dataTypes.insert( QGis::CInt32, "CInt32" );
843  dataTypes.insert( QGis::CFloat32, "CFloat32" );
844  dataTypes.insert( QGis::CFloat64, "CFloat64" );
845 
846  for ( int i = 1; i <= nBands; i++ )
847  {
848  QDomElement VRTBand = mVRTDocument.createElement( "VRTRasterBand" );
849 
850  VRTBand.setAttribute( "band", QString::number( i ) );
851  QString dataType = dataTypes.value( type );
852  VRTBand.setAttribute( "dataType", dataType );
853 
854  if ( mMode == Image )
855  {
856  VRTBand.setAttribute( "dataType", "Byte" );
857  QDomElement colorInterpElement = mVRTDocument.createElement( "ColorInterp" );
858  QDomText interpText = mVRTDocument.createTextNode( colorInterp.value( i - 1 ) );
859  colorInterpElement.appendChild( interpText );
860  VRTBand.appendChild( colorInterpElement );
861  }
862 
863  if ( !destHasNoDataValueList.isEmpty() && destHasNoDataValueList.value( i - 1 ) )
864  {
865  VRTBand.setAttribute( "NoDataValue", QString::number( destNoDataValueList.value( i - 1 ) ) );
866  }
867 
868  mVRTBands.append( VRTBand );
869  VRTDatasetElem.appendChild( VRTBand );
870  }
871 }
872 
873 bool QgsRasterFileWriter::writeVRT( const QString& file )
874 {
875  QFile outputFile( file );
876  if ( ! outputFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
877  {
878  return false;
879  }
880 
881  QTextStream outStream( &outputFile );
882  mVRTDocument.save( outStream, 2 );
883  return true;
884 }
885 
886 QgsRasterDataProvider* QgsRasterFileWriter::createPartProvider( const QgsRectangle& extent, int nCols, int iterCols,
887  int iterRows, int iterLeft, int iterTop, const QString& outputUrl, int fileIndex, int nBands, QGis::DataType type,
888  const QgsCoordinateReferenceSystem& crs )
889 {
890  double mup = extent.width() / nCols;
891  double mapLeft = extent.xMinimum() + iterLeft * mup;
892  double mapRight = mapLeft + mup * iterCols;
893  double mapTop = extent.yMaximum() - iterTop * mup;
894  double mapBottom = mapTop - iterRows * mup;
895  QgsRectangle mapRect( mapLeft, mapBottom, mapRight, mapTop );
896 
897  QString outputFile = outputUrl + '/' + partFileName( fileIndex );
898 
899  //geotransform
900  double geoTransform[6];
901  geoTransform[0] = mapRect.xMinimum();
902  geoTransform[1] = mup;
903  geoTransform[2] = 0.0;
904  geoTransform[3] = mapRect.yMaximum();
905  geoTransform[4] = 0.0;
906  geoTransform[5] = -mup;
907 
908  // perhaps we need a separate createOptions for tiles ?
909 
910  QgsRasterDataProvider* destProvider = QgsRasterDataProvider::create( mOutputProviderKey, outputFile, mOutputFormat, nBands, type, iterCols, iterRows, geoTransform, crs, mCreateOptions );
911 
912  // TODO: return provider and report error
913  return destProvider;
914 }
915 
916 QgsRasterDataProvider* QgsRasterFileWriter::initOutput( int nCols, int nRows, const QgsCoordinateReferenceSystem& crs,
917  double* geoTransform, int nBands, QGis::DataType type,
918  const QList<bool>& destHasNoDataValueList, const QList<double>& destNoDataValueList )
919 {
920  if ( mTiledMode )
921  {
922  createVRT( nCols, nRows, crs, geoTransform, type, destHasNoDataValueList, destNoDataValueList );
923  return nullptr;
924  }
925  else
926  {
927 #if 0
928  // TODO enable "use existing", has no effect for now, because using Create() in gdal provider
929  // should this belong in provider? should also test that source provider is gdal
930  if ( mBuildPyramidsFlag == -4 && mOutputProviderKey == "gdal" && mOutputFormat.toLower() == "gtiff" )
931  mCreateOptions << "COPY_SRC_OVERVIEWS=YES";
932 #endif
933 
934  QgsRasterDataProvider* destProvider = QgsRasterDataProvider::create( mOutputProviderKey, mOutputUrl, mOutputFormat, nBands, type, nCols, nRows, geoTransform, crs, mCreateOptions );
935 
936  if ( !destProvider )
937  {
938  QgsDebugMsg( "No provider created" );
939  }
940 
941  return destProvider;
942  }
943 }
944 
945 void QgsRasterFileWriter::globalOutputParameters( const QgsRectangle& extent, int nCols, int& nRows,
946  double* geoTransform, double& pixelSize )
947 {
948  pixelSize = extent.width() / nCols;
949 
950  //calculate nRows automatically for providers without exact resolution
951  if ( nRows < 0 )
952  {
953  nRows = static_cast< double >( nCols ) / extent.width() * extent.height() + 0.5;
954  }
955  geoTransform[0] = extent.xMinimum();
956  geoTransform[1] = pixelSize;
957  geoTransform[2] = 0.0;
958  geoTransform[3] = extent.yMaximum();
959  geoTransform[4] = 0.0;
960  geoTransform[5] = -( extent.height() / nRows );
961 }
962 
963 QString QgsRasterFileWriter::partFileName( int fileIndex )
964 {
965  // .tif for now
966  QFileInfo outputInfo( mOutputUrl );
967  return QString( "%1.%2.tif" ).arg( outputInfo.fileName() ).arg( fileIndex );
968 }
969 
970 QString QgsRasterFileWriter::vrtFileName()
971 {
972  QFileInfo outputInfo( mOutputUrl );
973  return QString( "%1.vrt" ).arg( outputInfo.fileName() );
974 }
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:228
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:489
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:258
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)