QGIS API Documentation  2.99.0-Master (7d4f81d)
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 "qgsfeedback.h"
22 #include <QProgressDialog>
23 #include <QFile>
24 
25 QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
26  : mInputFile( inputFile )
27  , mOutputFile( outputFile )
28  , mOutputFormat( outputFormat )
29  , mCellSizeX( -1.0 )
30  , mCellSizeY( -1.0 )
31  , mInputNodataValue( -1.0 )
32  , mOutputNodataValue( -1.0 )
33  , mZFactor( 1.0 )
34 {
35 
36 }
37 
39  : mCellSizeX( -1.0 )
40  , mCellSizeY( -1.0 )
41  , mInputNodataValue( -1.0 )
42  , mOutputNodataValue( -1.0 )
43  , mZFactor( 1.0 )
44 {
45 }
46 
48 {
49  GDALAllRegister();
50 
51  //open input file
52  int xSize, ySize;
53  GDALDatasetH inputDataset = openInputFile( xSize, ySize );
54  if ( !inputDataset )
55  {
56  return 1; //opening of input file failed
57  }
58 
59  //output driver
60  GDALDriverH outputDriver = openOutputDriver();
61  if ( !outputDriver )
62  {
63  return 2;
64  }
65 
66  GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
67  if ( !outputDataset )
68  {
69  return 3; //create operation on output file failed
70  }
71 
72  //open first raster band for reading (operation is only for single band raster)
73  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, 1 );
74  if ( !rasterBand )
75  {
76  GDALClose( inputDataset );
77  GDALClose( outputDataset );
78  return 4;
79  }
80  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
81 
82  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
83  if ( !outputRasterBand )
84  {
85  GDALClose( inputDataset );
86  GDALClose( outputDataset );
87  return 5;
88  }
89  //try to set -9999 as nodata value
90  GDALSetRasterNoDataValue( outputRasterBand, -9999 );
91  mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
92 
93  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
94  {
95  GDALClose( inputDataset );
96  GDALClose( outputDataset );
97  return 6;
98  }
99 
100  //keep only three scanlines in memory at a time
101  float *scanLine1 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
102  float *scanLine2 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
103  float *scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
104 
105  float *resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
106 
107  //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
108  for ( int i = 0; i < ySize; ++i )
109  {
110  if ( feedback && feedback->isCanceled() )
111  {
112  break;
113  }
114 
115  if ( feedback )
116  {
117  feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize );
118  }
119 
120  if ( i == 0 )
121  {
122  //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
123  for ( int a = 0; a < xSize; ++a )
124  {
125  scanLine1[a] = mInputNodataValue;
126  }
127  if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
128  {
129  QgsDebugMsg( "Raster IO Error" );
130  }
131  }
132  else
133  {
134  //normally fetch only scanLine3 and release scanline 1 if we move forward one row
135  CPLFree( scanLine1 );
136  scanLine1 = scanLine2;
137  scanLine2 = scanLine3;
138  scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
139  }
140 
141  if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
142  {
143  for ( int a = 0; a < xSize; ++a )
144  {
145  scanLine3[a] = mInputNodataValue;
146  }
147  }
148  else
149  {
150  if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
151  {
152  QgsDebugMsg( "Raster IO Error" );
153  }
154  }
155 
156  for ( int j = 0; j < xSize; ++j )
157  {
158  if ( j == 0 )
159  {
160  resultLine[j] = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j + 1], &mInputNodataValue, &scanLine2[j],
161  &scanLine2[j + 1], &mInputNodataValue, &scanLine3[j], &scanLine3[j + 1] );
162  }
163  else if ( j == xSize - 1 )
164  {
165  resultLine[j] = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &mInputNodataValue, &scanLine2[j - 1], &scanLine2[j],
166  &mInputNodataValue, &scanLine3[j - 1], &scanLine3[j], &mInputNodataValue );
167  }
168  else
169  {
170  resultLine[j] = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &scanLine1[j + 1], &scanLine2[j - 1], &scanLine2[j],
171  &scanLine2[j + 1], &scanLine3[j - 1], &scanLine3[j], &scanLine3[j + 1] );
172  }
173  }
174 
175  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
176  {
177  QgsDebugMsg( "Raster IO Error" );
178  }
179  }
180 
181  CPLFree( resultLine );
182  CPLFree( scanLine1 );
183  CPLFree( scanLine2 );
184  CPLFree( scanLine3 );
185 
186  GDALClose( inputDataset );
187 
188  if ( feedback && feedback->isCanceled() )
189  {
190  //delete the dataset without closing (because it is faster)
191  GDALDeleteDataset( outputDriver, mOutputFile.toUtf8().constData() );
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( mInputFile.toUtf8().constData(), GA_ReadOnly );
202  if ( inputDataset )
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 nullptr;
212  }
213  }
214  return inputDataset;
215 }
216 
217 GDALDriverH QgsNineCellFilter::openOutputDriver()
218 {
219  char **driverMetadata = nullptr;
220 
221  //open driver
222  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
223 
224  if ( !outputDriver )
225  {
226  return outputDriver; //return nullptr, driver does not exist
227  }
228 
229  driverMetadata = GDALGetMetadata( outputDriver, nullptr );
230  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
231  {
232  return nullptr; //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 )
241  {
242  return nullptr;
243  }
244 
245  int xSize = GDALGetRasterXSize( inputDataset );
246  int ySize = GDALGetRasterYSize( inputDataset );
247 
248  //open output file
249  char **papszOptions = nullptr;
250  GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 1, GDT_Float32, papszOptions );
251  if ( !outputDataset )
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 nullptr;
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 
#define QgsDebugMsg(str)
Definition: qgslogger.h:37
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:72
int processRaster(QgsFeedback *feedback=nullptr)
Starts the calculation, reads from mInputFile and stores the result in mOutputFile.
Base class for feedback objects to be used for cancelation of something running in a worker thread...
Definition: qgsfeedback.h:43
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 * GDALDatasetH
float mOutputNodataValue
The nodata value of the output layer.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:63
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.