Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
qgsninecellfilter.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsninecellfilter.h - description
3  -------------------
4  begin : August 6th, 2009
5  copyright : (C) 2009 by Marco Hugentobler
6  email : marco dot hugentobler at karto dot baug dot ethz dot ch
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #include "qgsninecellfilter.h"
19 #include "cpl_string.h"
20 #include <QProgressDialog>
21 
22 #if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
23 #define TO8(x) (x).toUtf8().constData()
24 #else
25 #define TO8(x) (x).toLocal8Bit().constData()
26 #endif
27 
28 QgsNineCellFilter::QgsNineCellFilter( const QString& inputFile, const QString& outputFile, const QString& outputFormat ): \
29  mInputFile( inputFile ), mOutputFile( outputFile ), mOutputFormat( outputFormat ), mCellSizeX( -1 ), mCellSizeY( -1 ), mInputNodataValue( -1 ), mOutputNodataValue( -1 )
30 {
31 
32 }
33 
35 {
36 
37 }
38 
40 {
41 
42 }
43 
44 int QgsNineCellFilter::processRaster( QProgressDialog* p )
45 {
46  GDALAllRegister();
47 
48  //open input file
49  int xSize, ySize;
50  GDALDatasetH inputDataset = openInputFile( xSize, ySize );
51  if ( inputDataset == NULL )
52  {
53  return 1; //opening of input file failed
54  }
55 
56  //output driver
57  GDALDriverH outputDriver = openOutputDriver();
58  if ( outputDriver == 0 )
59  {
60  return 2;
61  }
62 
63  GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
64  if ( outputDataset == NULL )
65  {
66  return 3; //create operation on output file failed
67  }
68 
69  //open first raster band for reading (operation is only for single band raster)
70  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
71  if ( rasterBand == NULL )
72  {
73  GDALClose( inputDataset );
74  GDALClose( outputDataset );
75  return 4;
76  }
77  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );
78 
79  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
80  if ( outputRasterBand == NULL )
81  {
82  GDALClose( inputDataset );
83  GDALClose( outputDataset );
84  return 5;
85  }
86  //try to set -9999 as nodata value
87  GDALSetRasterNoDataValue( outputRasterBand, -9999 );
88  mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, NULL );
89 
90  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
91  {
92  GDALClose( inputDataset );
93  GDALClose( outputDataset );
94  return 6;
95  }
96 
97  //keep only three scanlines in memory at a time
98  float* scanLine1 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
99  float* scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
100  float* scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
101 
102  float* resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
103 
104  if ( p )
105  {
106  p->setMaximum( ySize );
107  }
108 
109  //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
110  for ( int i = 0; i < ySize; ++i )
111  {
112  if ( p )
113  {
114  p->setValue( i );
115  }
116 
117  if ( p && p->wasCanceled() )
118  {
119  break;
120  }
121 
122  if ( i == 0 )
123  {
124  //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
125  for ( int a = 0; a < xSize; ++a )
126  {
127  scanLine1[a] = mInputNodataValue;
128  }
129  GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 );
130  }
131  else
132  {
133  //normally fetch only scanLine3 and release scanline 1 if we move forward one row
134  CPLFree( scanLine1 );
135  scanLine1 = scanLine2;
136  scanLine2 = scanLine3;
137  scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
138  }
139 
140  if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
141  {
142  for ( int a = 0; a < xSize; ++a )
143  {
144  scanLine3[a] = mInputNodataValue;
145  }
146  }
147  else
148  {
149  GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 );
150  }
151 
152  for ( int j = 0; j < xSize; ++j )
153  {
154  if ( j == 0 )
155  {
156  resultLine[j] = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j+1], &mInputNodataValue, &scanLine2[j], \
157  &scanLine2[j+1], &mInputNodataValue, &scanLine3[j], &scanLine3[j+1] );
158  }
159  else if ( j == xSize - 1 )
160  {
161  resultLine[j] = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &mInputNodataValue, &scanLine2[j-1], &scanLine2[j], \
162  &mInputNodataValue, &scanLine3[j-1], &scanLine3[j], &mInputNodataValue );
163  }
164  else
165  {
166  resultLine[j] = processNineCellWindow( &scanLine1[j-1], &scanLine1[j], &scanLine1[j+1], &scanLine2[j-1], &scanLine2[j], \
167  &scanLine2[j+1], &scanLine3[j-1], &scanLine3[j], &scanLine3[j+1] );
168  }
169  }
170 
171  GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 );
172  }
173 
174  if ( p )
175  {
176  p->setValue( ySize );
177  }
178 
179  CPLFree( resultLine );
180  CPLFree( scanLine1 );
181  CPLFree( scanLine2 );
182  CPLFree( scanLine3 );
183 
184  GDALClose( inputDataset );
185 
186  if ( p && p->wasCanceled() )
187  {
188  //delete the dataset without closing (because it is faster)
189  GDALDeleteDataset( outputDriver, mOutputFile.toLocal8Bit().data() );
190  return 7;
191  }
192  GDALClose( outputDataset );
193 
194  return 0;
195 }
196 
197 GDALDatasetH QgsNineCellFilter::openInputFile( int& nCellsX, int& nCellsY )
198 {
199  GDALDatasetH inputDataset = GDALOpen( TO8( mInputFile ), GA_ReadOnly );
200  if ( inputDataset != NULL )
201  {
202  nCellsX = GDALGetRasterXSize( inputDataset );
203  nCellsY = GDALGetRasterYSize( inputDataset );
204 
205  //we need at least one band
206  if ( GDALGetRasterCount( inputDataset ) < 1 )
207  {
208  GDALClose( inputDataset );
209  return NULL;
210  }
211  }
212  return inputDataset;
213 }
214 
216 {
217  char **driverMetadata;
218 
219  //open driver
220  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
221 
222  if ( outputDriver == NULL )
223  {
224  return outputDriver; //return NULL, driver does not exist
225  }
226 
227  driverMetadata = GDALGetMetadata( outputDriver, NULL );
228  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
229  {
230  return NULL; //driver exist, but it does not support the create operation
231  }
232 
233  return outputDriver;
234 }
235 
236 GDALDatasetH QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
237 {
238  if ( inputDataset == NULL )
239  {
240  return NULL;
241  }
242 
243  int xSize = GDALGetRasterXSize( inputDataset );
244  int ySize = GDALGetRasterYSize( inputDataset );;
245 
246  //open output file
247  char **papszOptions = NULL;
248  GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), xSize, ySize, 1, GDT_Float32, papszOptions );
249  if ( outputDataset == NULL )
250  {
251  return outputDataset;
252  }
253 
254  //get geotransform from inputDataset
255  double geotransform[6];
256  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
257  {
258  GDALClose( outputDataset );
259  return NULL;
260  }
261  GDALSetGeoTransform( outputDataset, geotransform );
262 
263  //make sure mCellSizeX and mCellSizeY are always > 0
264  mCellSizeX = geotransform[1];
265  if ( mCellSizeX < 0 )
266  {
268  }
269  mCellSizeY = geotransform[5];
270  if ( mCellSizeY < 0 )
271  {
273  }
274 
275  const char* projection = GDALGetProjectionRef( inputDataset );
276  GDALSetProjection( outputDataset, projection );
277 
278  return outputDataset;
279 }
280