QGIS API Documentation  2.11.0-Master
qgsrasterfilewriter.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrasterfilewriter.cpp
3  ---------------------
4  begin : July 2012
5  copyright : (C) 2012 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 #include <typeinfo>
16 
17 #include "qgsrasterfilewriter.h"
18 #include "qgsproviderregistry.h"
19 #include "qgsrasterinterface.h"
20 #include "qgsrasteriterator.h"
21 #include "qgsrasterlayer.h"
22 #include "qgsrasterprojector.h"
23 
24 #include <QCoreApplication>
25 #include <QProgressDialog>
26 #include <QTextStream>
27 #include <QMessageBox>
28 
30  mMode( Raw ), mOutputUrl( outputUrl ), mOutputProviderKey( "gdal" ), mOutputFormat( "GTiff" ),
31  mTiledMode( false ), mMaxTileWidth( 500 ), mMaxTileHeight( 500 ),
32  mBuildPyramidsFlag( QgsRaster::PyramidsFlagNo ),
33  mPyramidsFormat( QgsRaster::PyramidsGTiff ),
34  mProgressDialog( 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 
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.
Sixty four bit floating point (double)
Definition: qgis.h:213
Sixteen bit signed integer (qint16)
Definition: qgis.h:209
A rectangle specified with double values.
Definition: qgsrectangle.h:35
QgsRasterFileWriter(const QString &outputUrl)
Complex Float64.
Definition: qgis.h:217
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:202
Iterator for sequentially processing raster cells.
QDomNode appendChild(const QDomNode &newChild)
void startRasterRead(int bandNumber, int nCols, int nRows, const QgsRectangle &extent)
Start reading of raster band.
void setMaximum(int maximum)
void push_back(const T &value)
Eight bit unsigned integer (quint8)
Definition: qgis.h:207
void setWindowModality(Qt::WindowModality windowModality)
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:192
Sixteen bit unsigned integer (quint16)
Definition: qgis.h:208
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
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
A helper to get an wkt representation of this srs.
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.
QString tr(const char *sourceText, const char *disambiguation, int n)
int size() const
bool isNull() const
T value(int i) const
int maximumTileWidth() const
virtual QgsRasterBandStats bandStatistics(int theBandNo, int theStats=QgsRasterBandStats::All, const QgsRectangle &theExtent=QgsRectangle(), int theSampleSize=0)
Get band statistics.
void setValue(int progress)
static QGis::DataType typeWithNoDataValue(QGis::DataType dataType, double *noDataValue)
For given data type returns wider type and sets no data value.
QString number(int n, int base)
int count(const T &value) const
The RasterBandStats struct is a container for statistics about a single raster band.
void processEvents(QFlags< QEventLoop::ProcessEventsFlag > flags)
void append(const T &value)
Thirty two bit floating point (float)
Definition: qgis.h:212
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.
QString fileName() const
void setAttribute(const QString &name, const QString &value)
bool isEmpty() const
static bool typeIsColor(QGis::DataType type)
Returns true if data type is color.
int maximumTileHeight() const
static int typeSize(int dataType)
QDir dir() const
DataType
Raster data types.
Definition: qgis.h: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...
void * GDALDatasetH
void setRange(int minimum, int maximum)
Base class for processing filters like renderers, reprojector, resampler etc.
QDomText createTextNode(const QString &value)
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
Definition: qgis.h:449
QString toLower() const
QByteArray toLocal8Bit() const
bool exists() const
void setOutputNoDataValue(int bandNo, double noData)
Set output no data value.
Complex Int16.
Definition: qgis.h:214
virtual QgsRectangle extent() override=0
Get the extent of the data source.
void setMaximumTileWidth(int w)
virtual QGis::DataType dataType(int bandNo) const override=0
Returns data type for the band specified by number.
Raster namespace.
Definition: qgsraster.h:28
void clear()
void setMaximumTileHeight(int h)
void save(QTextStream &str, int indent) const
QgsCoordinateReferenceSystem srcCrs() const
Get source CRS.
QgsRasterRangeList noData(int bandNo) const
QString absolutePath() const
virtual bool remove()
Remove dataset.
bool mkdir(const QString &dirName) const
Class for storing a coordinate reference system (CRS)
Class for doing transforms between two map coordinate systems.
char * data()
double minimumValue
The minimum cell value in the raster band.
virtual QString buildPyramids(const QList< QgsRasterPyramid > &thePyramidList, const QString &theResamplingMethod="NEAREST", QgsRaster::RasterPyramidsFormat theFormat=QgsRaster::PyramidsGTiff, const QStringList &theConfigOptions=QStringList())
Create pyramid overviews.
StandardButton warning(QWidget *parent, const QString &title, const QString &text, QFlags< QMessageBox::StandardButton > buttons, StandardButton defaultButton)
iterator insert(const Key &key, const T &value)
void show()
virtual bool write(void *data, int band, int width, int height, int xOffset, int yOffset)
Writes into the provider datasource.
Complex Int32.
Definition: qgis.h:215
virtual bool srcHasNoDataValue(int bandNo) const
Return true if source band has no data value.
QRgb color(int row, int column) const
Read a single color.
QDomElement createElement(const QString &tagName)
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32.
Definition: qgis.h:219
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:202
Color, alpha, red, green, blue, 4 bytes the same as QImage::Format_ARGB32_Premultiplied.
Definition: qgis.h:221
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:187
void qgsFree(void *ptr)
Frees the memory space pointed to by ptr.
Definition: qgis.cpp:232
const QgsRasterInterface * input() const
Thirty two bit unsigned integer (quint32)
Definition: qgis.h:210
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:207
Complex Float32.
Definition: qgis.h:216
const T value(const Key &key) const
Base class for raster data providers.
void replace(int i, const T &value)