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