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