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