Quantum GIS API Documentation
1.7.4
|
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