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