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