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