QGIS API Documentation
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 
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.
Eight bit unsigned integer (quint8)
Definition: qgis.h:132
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:141
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:140
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:137
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:133
Sixty four bit floating point (double)
Definition: qgis.h:138
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
Definition: qgis.h:146
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:135
int maximumTileHeight() const
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32.
Definition: qgis.h:144
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:134
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:459
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:129
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:139
const QgsRasterInterface * input() const
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212
Complex Float64.
Definition: qgis.h:142
const T value(const Key &key) const
Base class for raster data providers.
void replace(int i, const T &value)