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