Use GDALAutoCreateWarpedVRT to handle south-up (or rotated) rasters a…

…lso for calculator. Fixes bug #3281

git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@14842 c8812cc2-4d05-0410-92ff-de0c093fc19c

src/analysis/raster/qgsrastercalculator.cpp

@@ -22,9 +22,11 @@
 #include "cpl_string.h"
 #include <QProgressDialog>
 
+#include "gdalwarper.h"
+
 QgsRasterCalculator::QgsRasterCalculator( const QString& formulaString, const QString& outputFile, const QString& outputFormat,
     const QgsRectangle& outputExtent, int nOutputColumns, int nOutputRows, const QVector<QgsRasterCalculatorEntry>& rasterEntries ): mFormulaString( formulaString ), mOutputFile( outputFile ), mOutputFormat( outputFormat ),
-  mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries )
+    mOutputRectangle( outputExtent ), mNumOutputColumns( nOutputColumns ), mNumOutputRows( nOutputRows ), mRasterEntries( rasterEntries )
 {
 }
 
@@ -37,7 +39,7 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
   //prepare search string / tree
   QString errorString;
   QgsRasterCalcNode* calcNode = QgsRasterCalcNode::parseRasterCalcString( mFormulaString, errorString );
-  if( !calcNode )
+  if ( !calcNode )
   {
     //error
   }
@@ -51,34 +53,53 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
   QVector< GDALDatasetH > mInputDatasets; //raster references and corresponding dataset
 
   QVector<QgsRasterCalculatorEntry>::const_iterator it = mRasterEntries.constBegin();
-  for( ; it != mRasterEntries.constEnd(); ++it )
+  for ( ; it != mRasterEntries.constEnd(); ++it )
   {
-    if( !it->raster )  // no raster layer in entry
+    if ( !it->raster ) // no raster layer in entry
     {
       return 2;
     }
     GDALDatasetH inputDataset = GDALOpen( it->raster->source().toLocal8Bit().data(), GA_ReadOnly );
-    if( inputDataset == NULL )
+    if ( inputDataset == NULL )
     {
       return 2;
     }
+
+    //check if the input dataset is south up or rotated. If yes, use GDALAutoCreateWarpedVRT to create a north up raster
+    double inputGeoTransform[6];
+    if ( GDALGetGeoTransform( inputDataset, inputGeoTransform ) == CE_None
+         && ( inputGeoTransform[1] < 0.0
+              || inputGeoTransform[2] != 0.0
+              || inputGeoTransform[4] != 0.0
+              || inputGeoTransform[5] > 0.0 ) )
+    {
+      GDALDatasetH vDataset = GDALAutoCreateWarpedVRT( inputDataset, NULL, NULL, GRA_NearestNeighbour, 0.2, NULL );
+      mInputDatasets.push_back( vDataset );
+      mInputDatasets.push_back( inputDataset );
+      inputDataset = vDataset;
+    }
+    else
+    {
+      mInputDatasets.push_back( inputDataset );
+    }
+
+
     GDALRasterBandH inputRasterBand = GDALGetRasterBand( inputDataset, it->bandNumber );
-    if( inputRasterBand == NULL )
+    if ( inputRasterBand == NULL )
     {
       return 2;
     }
 
     int nodataSuccess;
     double nodataValue = GDALGetRasterNoDataValue( inputRasterBand, &nodataSuccess );
 
-    mInputDatasets.push_back( inputDataset );
     mInputRasterBands.insert( it->ref, inputRasterBand );
     inputScanLineData.insert( it->ref, new QgsRasterMatrix( mNumOutputColumns, 1, new float[mNumOutputColumns], nodataValue ) );
   }
 
   //open output dataset for writing
   GDALDriverH outputDriver = openOutputDriver();
-  if( outputDriver == NULL )
+  if ( outputDriver == NULL )
   {
     return 1;
   }
@@ -90,29 +111,29 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
 
   float* resultScanLine = ( float * ) CPLMalloc( sizeof( float ) * mNumOutputColumns );
 
-  if( p )
+  if ( p )
   {
     p->setMaximum( mNumOutputRows );
   }
 
   QgsRasterMatrix resultMatrix;
 
   //read / write line by line
-  for( int i = 0; i < mNumOutputRows; ++i )
+  for ( int i = 0; i < mNumOutputRows; ++i )
   {
-    if( p )
+    if ( p )
     {
       p->setValue( i );
     }
 
-    if( p && p->wasCanceled() )
+    if ( p && p->wasCanceled() )
     {
       break;
     }
 
     //fill buffers
     QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
-    for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
+    for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
     {
       double sourceTransformation[6];
       GDALRasterBandH sourceRasterBand = mInputRasterBands[bufferIt.key()];
@@ -121,15 +142,15 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
       readRasterPart( targetGeoTransform, 0, i, mNumOutputColumns, 1, sourceTransformation, sourceRasterBand, bufferIt.value()->data() );
     }
 
-    if( calcNode->calculate( inputScanLineData, resultMatrix ) )
+    if ( calcNode->calculate( inputScanLineData, resultMatrix ) )
     {
       bool resultIsNumber = resultMatrix.isNumber();
       float* calcData;
 
-      if( resultIsNumber )  //scalar result. Insert number for every pixel
+      if ( resultIsNumber ) //scalar result. Insert number for every pixel
       {
         calcData = new float[mNumOutputColumns];
-        for( int j = 0; j < mNumOutputColumns; ++j )
+        for ( int j = 0; j < mNumOutputColumns; ++j )
         {
           calcData[j] = resultMatrix.number();
         }
@@ -140,49 +161,49 @@ int QgsRasterCalculator::processCalculation( QProgressDialog* p )
       }
 
       //replace all matrix nodata values with output nodatas
-      for( int j = 0; j < mNumOutputColumns; ++j )
+      for ( int j = 0; j < mNumOutputColumns; ++j )
       {
-        if( calcData[j] == resultMatrix.nodataValue() )
+        if ( calcData[j] == resultMatrix.nodataValue() )
         {
           calcData[j] = outputNodataValue;
         }
       }
 
       //write scanline to the dataset
-      if( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
+      if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, mNumOutputColumns, 1, calcData, mNumOutputColumns, 1, GDT_Float32, 0, 0 ) != CE_None )
       {
         qWarning( "RasterIO error!" );
       }
 
-      if( resultIsNumber )
+      if ( resultIsNumber )
       {
         delete[] calcData;
       }
     }
 
   }
 
-  if( p )
+  if ( p )
   {
     p->setValue( mNumOutputRows );
   }
 
   //close datasets and release memory
   delete calcNode;
   QMap< QString, QgsRasterMatrix* >::iterator bufferIt = inputScanLineData.begin();
-  for( ; bufferIt != inputScanLineData.end(); ++bufferIt )
+  for ( ; bufferIt != inputScanLineData.end(); ++bufferIt )
   {
     delete bufferIt.value();
   }
   inputScanLineData.clear();
 
   QVector< GDALDatasetH >::iterator datasetIt = mInputDatasets.begin();
-  for( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
+  for ( ; datasetIt != mInputDatasets.end(); ++ datasetIt )
   {
     GDALClose( *datasetIt );
   }
 
-  if( p && p->wasCanceled() )
+  if ( p && p->wasCanceled() )
   {
     //delete the dataset without closing (because it is faster)
     GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() );
@@ -204,13 +225,13 @@ GDALDriverH QgsRasterCalculator::openOutputDriver()
   //open driver
   GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
 
-  if( outputDriver == NULL )
+  if ( outputDriver == NULL )
   {
     return outputDriver; //return NULL, driver does not exist
   }
 
   driverMetadata = GDALGetMetadata( outputDriver, NULL );
-  if( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
+  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
   {
     return NULL; //driver exist, but it does not support the create operation
   }
@@ -223,7 +244,7 @@ GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
   //open output file
   char **papszOptions = NULL;
   GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );
-  if( outputDataset == NULL )
+  if ( outputDataset == NULL )
   {
     return outputDataset;
   }
@@ -239,7 +260,7 @@ GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver )
 void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffset, int yOffset, int nCols, int nRows, double* sourceTransform, GDALRasterBandH sourceBand, float* rasterBuffer )
 {
   //If dataset transform is the same as the requested transform, do a normal GDAL raster io
-  if( transformationsEqual( targetGeotransform, sourceTransform ) )
+  if ( transformationsEqual( targetGeotransform, sourceTransform ) )
   {
     GDALRasterIO( sourceBand, GF_Read, xOffset, yOffset, nCols, nRows, rasterBuffer, nCols, nRows, GDT_Float32, 0, 0 );
     return;
@@ -255,10 +276,10 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse
   QgsRectangle intersection = targetRect.intersect( &sourceRect );
 
   //no intersection, fill all the pixels with nodata values
-  if( intersection.isEmpty() )
+  if ( intersection.isEmpty() )
   {
     int nPixels = nCols * nRows;
-    for( int i = 0; i < nPixels; ++i )
+    for ( int i = 0; i < nPixels; ++i )
     {
       rasterBuffer[i] = nodataValue;
     }
@@ -284,17 +305,17 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse
   double targetPixelY = targetGeotransform[3] + targetGeotransform[5] * yOffset + targetGeotransform[5] / 2.0; //coordinates of current target pixel
   int sourceIndexX, sourceIndexY; //current raster index in  source pixels
   double sx, sy;
-  for( int i = 0; i < nRows; ++i )
+  for ( int i = 0; i < nRows; ++i )
   {
     targetPixelX = targetPixelXMin;
-    for( int j = 0; j < nCols; ++j )
+    for ( int j = 0; j < nCols; ++j )
     {
       sx = ( targetPixelX - sourceRasterXMin ) / sourceTransform[1];
       sourceIndexX = sx > 0 ? sx : floor( sx );
       sy = ( targetPixelY - sourceRasterYMax ) / sourceTransform[5];
       sourceIndexY = sy > 0 ? sy : floor( sy );
-      if( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
-          && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
+      if ( sourceIndexX >= 0 && sourceIndexX < nSourcePixelsX
+           && sourceIndexY >= 0 && sourceIndexY < nSourcePixelsY )
       {
         rasterBuffer[j + i*nRows] = sourceRaster[ sourceIndexX  + nSourcePixelsX * sourceIndexY ];
       }
@@ -313,9 +334,9 @@ void QgsRasterCalculator::readRasterPart( double* targetGeotransform, int xOffse
 
 bool QgsRasterCalculator::transformationsEqual( double* t1, double* t2 ) const
 {
-  for( int i = 0; i < 6; ++i )
+  for ( int i = 0; i < 6; ++i )
   {
-    if( !doubleNear( t1[i], t2[i], 0.00001 ) )
+    if ( !doubleNear( t1[i], t2[i], 0.00001 ) )
     {
       return false;
     }