Quantum GIS API Documentation  1.7.4
src/analysis/raster/qgsninecellfilter.cpp
Go to the documentation of this file.
00001 /***************************************************************************
00002                           qgsninecellfilter.h  -  description
00003                              -------------------
00004     begin                : August 6th, 2009
00005     copyright            : (C) 2009 by Marco Hugentobler
00006     email                : marco dot hugentobler at karto dot baug dot ethz dot ch
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #include "qgsninecellfilter.h"
00019 #include "cpl_string.h"
00020 #include <QProgressDialog>
00021 
00022 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
00023 #define TO8(x) (x).toUtf8().constData()
00024 #else
00025 #define TO8(x) (x).toLocal8Bit().constData()
00026 #endif
00027 
00028 QgsNineCellFilter::QgsNineCellFilter( const QString& inputFile, const QString& outputFile, const QString& outputFormat ): \
00029     mInputFile( inputFile ), mOutputFile( outputFile ), mOutputFormat( outputFormat ), mCellSizeX( -1 ), mCellSizeY( -1 ), mInputNodataValue( -1 ), mOutputNodataValue( -1 )
00030 {
00031 
00032 }
00033 
00034 QgsNineCellFilter::QgsNineCellFilter()
00035 {
00036 
00037 }
00038 
00039 QgsNineCellFilter::~QgsNineCellFilter()
00040 {
00041 
00042 }
00043 
00044 int QgsNineCellFilter::processRaster( QProgressDialog* p )
00045 {
00046   GDALAllRegister();
00047 
00048   //open input file
00049   int xSize, ySize;
00050   GDALDatasetH  inputDataset = openInputFile( xSize, ySize );
00051   if ( inputDataset == NULL )
00052   {
00053     return 1; //opening of input file failed
00054   }
00055 
00056   //output driver
00057   GDALDriverH outputDriver = openOutputDriver();
00058   if ( outputDriver == 0 )
00059   {
00060     return 2;
00061   }
00062 
00063   GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
00064   if ( outputDataset == NULL )
00065   {
00066     return 3; //create operation on output file failed
00067   }
00068 
00069   //open first raster band for reading (operation is only for single band raster)
00070   GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
00071   if ( rasterBand == NULL )
00072   {
00073     GDALClose( inputDataset );
00074     GDALClose( outputDataset );
00075     return 4;
00076   }
00077   mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
00078 
00079   GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
00080   if ( outputRasterBand == NULL )
00081   {
00082     GDALClose( inputDataset );
00083     GDALClose( outputDataset );
00084     return 5;
00085   }
00086   //try to set -9999 as nodata value
00087   GDALSetRasterNoDataValue( outputRasterBand, -9999 );
00088   mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, NULL );
00089 
00090   if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
00091   {
00092     GDALClose( inputDataset );
00093     GDALClose( outputDataset );
00094     return 6;
00095   }
00096 
00097   //keep only three scanlines in memory at a time
00098   float* scanLine1 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
00099   float* scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
00100   float* scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
00101 
00102   float* resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
00103 
00104   if ( p )
00105   {
00106     p->setMaximum( ySize );
00107   }
00108 
00109   //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
00110   for ( int i = 0; i < ySize; ++i )
00111   {
00112     if ( p )
00113     {
00114       p->setValue( i );
00115     }
00116 
00117     if ( p && p->wasCanceled() )
00118     {
00119       break;
00120     }
00121 
00122     if ( i == 0 )
00123     {
00124       //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
00125       for ( int a = 0; a < xSize; ++a )
00126       {
00127         scanLine1[a] = mInputNodataValue;
00128       }
00129       GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 );
00130     }
00131     else
00132     {
00133       //normally fetch only scanLine3 and release scanline 1 if we move forward one row
00134       CPLFree( scanLine1 );
00135       scanLine1 = scanLine2;
00136       scanLine2 = scanLine3;
00137       scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
00138     }
00139 
00140     if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
00141     {
00142       for ( int a = 0; a < xSize; ++a )
00143       {
00144         scanLine3[a] = mInputNodataValue;
00145       }
00146     }
00147     else
00148     {
00149       GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 );
00150     }
00151 
00152     for ( int j = 0; j < xSize; ++j )
00153     {
00154       if ( j == 0 )
00155       {
00156         resultLine[j] = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j+1], &mInputNodataValue, &scanLine2[j], \
00157                                                &scanLine2[j+1], &mInputNodataValue, &scanLine3[j], &scanLine3[j+1] );
00158       }
00159       else if ( j == xSize - 1 )
00160       {
00161         resultLine[j] = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &mInputNodataValue, &scanLine2[j-1], &scanLine2[j], \
00162                                                &mInputNodataValue, &scanLine3[j-1], &scanLine3[j], &mInputNodataValue );
00163       }
00164       else
00165       {
00166         resultLine[j] = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &scanLine1[j+1], &scanLine2[j-1], &scanLine2[j], \
00167                                                &scanLine2[j+1], &scanLine3[j-1], &scanLine3[j], &scanLine3[j+1] );
00168       }
00169     }
00170 
00171     GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 );
00172   }
00173 
00174   if ( p )
00175   {
00176     p->setValue( ySize );
00177   }
00178 
00179   CPLFree( resultLine );
00180   CPLFree( scanLine1 );
00181   CPLFree( scanLine2 );
00182   CPLFree( scanLine3 );
00183 
00184   GDALClose( inputDataset );
00185 
00186   if ( p && p->wasCanceled() )
00187   {
00188     //delete the dataset without closing (because it is faster)
00189     GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() );
00190     return 7;
00191   }
00192   GDALClose( outputDataset );
00193 
00194   return 0;
00195 }
00196 
00197 GDALDatasetH QgsNineCellFilter::openInputFile( int& nCellsX, int& nCellsY )
00198 {
00199   GDALDatasetH inputDataset = GDALOpen( TO8( mInputFile ), GA_ReadOnly );
00200   if ( inputDataset != NULL )
00201   {
00202     nCellsX = GDALGetRasterXSize( inputDataset );
00203     nCellsY = GDALGetRasterYSize( inputDataset );
00204 
00205     //we need at least one band
00206     if ( GDALGetRasterCount( inputDataset ) < 1 )
00207     {
00208       GDALClose( inputDataset );
00209       return NULL;
00210     }
00211   }
00212   return inputDataset;
00213 }
00214 
00215 GDALDriverH QgsNineCellFilter::openOutputDriver()
00216 {
00217   char **driverMetadata;
00218 
00219   //open driver
00220   GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
00221 
00222   if ( outputDriver == NULL )
00223   {
00224     return outputDriver; //return NULL, driver does not exist
00225   }
00226 
00227   driverMetadata = GDALGetMetadata( outputDriver, NULL );
00228   if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
00229   {
00230     return NULL; //driver exist, but it does not support the create operation
00231   }
00232 
00233   return outputDriver;
00234 }
00235 
00236 GDALDatasetH QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
00237 {
00238   if ( inputDataset == NULL )
00239   {
00240     return NULL;
00241   }
00242 
00243   int xSize = GDALGetRasterXSize( inputDataset );
00244   int ySize = GDALGetRasterYSize( inputDataset );;
00245 
00246   //open output file
00247   char **papszOptions = NULL;
00248   GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), xSize, ySize, 1, GDT_Float32, papszOptions );
00249   if ( outputDataset == NULL )
00250   {
00251     return outputDataset;
00252   }
00253 
00254   //get geotransform from inputDataset
00255   double geotransform[6];
00256   if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
00257   {
00258     GDALClose( outputDataset );
00259     return NULL;
00260   }
00261   GDALSetGeoTransform( outputDataset, geotransform );
00262 
00263   //make sure mCellSizeX and mCellSizeY are always > 0
00264   mCellSizeX = geotransform[1];
00265   if ( mCellSizeX < 0 )
00266   {
00267     mCellSizeX = -mCellSizeX;
00268   }
00269   mCellSizeY = geotransform[5];
00270   if ( mCellSizeY < 0 )
00271   {
00272     mCellSizeY = -mCellSizeY;
00273   }
00274 
00275   const char* projection = GDALGetProjectionRef( inputDataset );
00276   GDALSetProjection( outputDataset, projection );
00277 
00278   return outputDataset;
00279 }
00280 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines