QGIS API Documentation  2.3.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 }
QList< QDomElement > mVRTBands
virtual int bandCount() const =0
Get number of bands.
A rectangle specified with double values.
Definition: qgsrectangle.h:35
QgsRaster::RasterBuildPyramids mBuildPyramidsFlag
bool mTiledMode
False: Write one file, true: create a directory and add the files numbered.
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:175
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:189
#define QgsDebugMsg(str)
Definition: qgslogger.h:36
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.
const QgsRasterPipe * mPipe
bool writeVRT(const QString &file)
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.
QString partFileName(int fileIndex)
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:37
virtual QList< QgsRasterPyramid > buildPyramidList(QList< int > overviewList=QList< int >())
Accessor for ths raster layers pyramid list.
virtual QGis::DataType srcDataType(int bandNo) const =0
Returns source data type for the band specified by number, source data type may be shorter than dataT...
virtual QgsRectangle extent()=0
Get the extent of the data source.
static bool typeIsColor(QGis::DataType type)
Returns true if data type is color.
void buildPyramids(const QString &filename)
virtual QGis::DataType dataType(int bandNo) const =0
Returns data type for the band specified by number.
QgsRaster::RasterPyramidsFormat mPyramidsFormat
int maximumTileHeight() const
static int typeSize(int dataType)
QgsRasterProjector * projector() const
QgsCoordinateReferenceSystem destCrs() const
Get destination CRS.
QgsRasterDataProvider * initOutput(int nCols, int nRows, const QgsCoordinateReferenceSystem &crs, double *geoTransform, int nBands, QGis::DataType type, QList< bool > destHasNoDataValueList=QList< bool >(), QList< double > destNoDataValueList=QList< double >())
Init VRT (for tiled mode) or create global output provider (single-file mode)
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.
WriterError writeDataRaster(const QgsRasterPipe *pipe, QgsRasterIterator *iter, int nCols, int nRows, const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &crs, QProgressDialog *progressDialog=0)
const QgsRasterInterface * mInput
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:424
QProgressDialog * mProgressDialog
void setOutputNoDataValue(int bandNo, double noData)
Set output no data value.
QString file
Definition: qgssvgcache.cpp:76
void setMaximumTileWidth(int w)
Raster namespace.
Definition: qgsraster.h:28
WriterError writeImageRaster(QgsRasterIterator *iter, int nCols, int nRows, const QgsRectangle &outputExtent, const QgsCoordinateReferenceSystem &crs, QProgressDialog *progressDialog=0)
void setMaximumTileHeight(int h)
QgsCoordinateReferenceSystem srcCrs() const
Get source CRS.
QgsRasterRangeList noData(int bandNo) const
virtual bool remove()
Returns the formats supported by create()
void createVRT(int xSize, int ySize, const QgsCoordinateReferenceSystem &crs, double *geoTransform, QGis::DataType type, QList< bool > destHasNoDataValueList, QList< double > destNoDataValueList)
Initialize vrt member variables.
Class for storing a coordinate reference system (CRS)
DataType
Raster data types.
Definition: qgis.h:204
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.
void globalOutputParameters(const QgsRectangle &extent, int nCols, int &nRows, double *geoTransform, double &pixelSize)
Calculate nRows, geotransform and pixel size for output.
void addToVRT(const QString &filename, int band, int xSize, int ySize, int xOffset, int yOffset)
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:199
QgsRasterDataProvider * createPartProvider(const QgsRectangle &extent, int nCols, int iterCols, int iterRows, int iterLeft, int iterTop, const QString &outputUrl, int fileIndex, int nBands, QGis::DataType type, const QgsCoordinateReferenceSystem &crs)
Create provider and datasource for a part image (vrt mode)
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:184
void qgsFree(void *ptr)
Frees the memory space pointed to by ptr.
Definition: qgis.cpp:205
const QgsRasterInterface * input() const
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:204
QStringList mPyramidsConfigOptions
Base class for raster data providers.
#define tr(sourceText)