QGIS API Documentation  2.99.0-Master (53aba61)
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 <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  GDALAllRegister();
49 
50  //open input file
51  int xSize, ySize;
52  GDALDatasetH inputDataset = openInputFile( xSize, ySize );
53  if ( !inputDataset )
54  {
55  return 1; //opening of input file failed
56  }
57 
58  //output driver
59  GDALDriverH outputDriver = openOutputDriver();
60  if ( !outputDriver )
61  {
62  return 2;
63  }
64 
65  GDALDatasetH outputDataset = openOutputFile( inputDataset, outputDriver );
66  if ( !outputDataset )
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 )
74  {
75  GDALClose( inputDataset );
76  GDALClose( outputDataset );
77  return 4;
78  }
79  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
80 
81  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset, 1 );
82  if ( !outputRasterBand )
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, nullptr );
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  //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
107  for ( int i = 0; i < ySize; ++i )
108  {
109  if ( feedback && feedback->isCanceled() )
110  {
111  break;
112  }
113 
114  if ( feedback )
115  {
116  feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize );
117  }
118 
119  if ( i == 0 )
120  {
121  //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
122  for ( int a = 0; a < xSize; ++a )
123  {
124  scanLine1[a] = mInputNodataValue;
125  }
126  if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, scanLine2, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
127  {
128  QgsDebugMsg( "Raster IO Error" );
129  }
130  }
131  else
132  {
133  //normally fetch only scanLine3 and release scanline 1 if we move forward one row
134  CPLFree( scanLine1 );
135  scanLine1 = scanLine2;
136  scanLine2 = scanLine3;
137  scanLine3 = ( float * ) CPLMalloc( sizeof( float ) * xSize );
138  }
139 
140  if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
141  {
142  for ( int a = 0; a < xSize; ++a )
143  {
144  scanLine3[a] = mInputNodataValue;
145  }
146  }
147  else
148  {
149  if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, scanLine3, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
150  {
151  QgsDebugMsg( "Raster IO Error" );
152  }
153  }
154 
155  for ( int j = 0; j < xSize; ++j )
156  {
157  if ( j == 0 )
158  {
159  resultLine[j] = processNineCellWindow( &mInputNodataValue, &scanLine1[j], &scanLine1[j + 1], &mInputNodataValue, &scanLine2[j],
160  &scanLine2[j + 1], &mInputNodataValue, &scanLine3[j], &scanLine3[j + 1] );
161  }
162  else if ( j == xSize - 1 )
163  {
164  resultLine[j] = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &mInputNodataValue, &scanLine2[j - 1], &scanLine2[j],
165  &mInputNodataValue, &scanLine3[j - 1], &scanLine3[j], &mInputNodataValue );
166  }
167  else
168  {
169  resultLine[j] = processNineCellWindow( &scanLine1[j - 1], &scanLine1[j], &scanLine1[j + 1], &scanLine2[j - 1], &scanLine2[j],
170  &scanLine2[j + 1], &scanLine3[j - 1], &scanLine3[j], &scanLine3[j + 1] );
171  }
172  }
173 
174  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
175  {
176  QgsDebugMsg( "Raster IO Error" );
177  }
178  }
179 
180  CPLFree( resultLine );
181  CPLFree( scanLine1 );
182  CPLFree( scanLine2 );
183  CPLFree( scanLine3 );
184 
185  GDALClose( inputDataset );
186 
187  if ( feedback && feedback->isCanceled() )
188  {
189  //delete the dataset without closing (because it is faster)
190  GDALDeleteDataset( outputDriver, mOutputFile.toUtf8().constData() );
191  return 7;
192  }
193  GDALClose( outputDataset );
194 
195  return 0;
196 }
197 
198 GDALDatasetH QgsNineCellFilter::openInputFile( int &nCellsX, int &nCellsY )
199 {
200  GDALDatasetH inputDataset = GDALOpen( mInputFile.toUtf8().constData(), GA_ReadOnly );
201  if ( inputDataset )
202  {
203  nCellsX = GDALGetRasterXSize( inputDataset );
204  nCellsY = GDALGetRasterYSize( inputDataset );
205 
206  //we need at least one band
207  if ( GDALGetRasterCount( inputDataset ) < 1 )
208  {
209  GDALClose( inputDataset );
210  return nullptr;
211  }
212  }
213  return inputDataset;
214 }
215 
216 GDALDriverH QgsNineCellFilter::openOutputDriver()
217 {
218  char **driverMetadata = nullptr;
219 
220  //open driver
221  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
222 
223  if ( !outputDriver )
224  {
225  return outputDriver; //return nullptr, driver does not exist
226  }
227 
228  driverMetadata = GDALGetMetadata( outputDriver, nullptr );
229  if ( !CSLFetchBoolean( driverMetadata, GDAL_DCAP_CREATE, false ) )
230  {
231  return nullptr; //driver exist, but it does not support the create operation
232  }
233 
234  return outputDriver;
235 }
236 
237 GDALDatasetH QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
238 {
239  if ( !inputDataset )
240  {
241  return nullptr;
242  }
243 
244  int xSize = GDALGetRasterXSize( inputDataset );
245  int ySize = GDALGetRasterYSize( inputDataset );
246 
247  //open output file
248  char **papszOptions = nullptr;
249  GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 1, GDT_Float32, papszOptions );
250  if ( !outputDataset )
251  {
252  return outputDataset;
253  }
254 
255  //get geotransform from inputDataset
256  double geotransform[6];
257  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
258  {
259  GDALClose( outputDataset );
260  return nullptr;
261  }
262  GDALSetGeoTransform( outputDataset, geotransform );
263 
264  //make sure mCellSizeX and mCellSizeY are always > 0
265  mCellSizeX = geotransform[1];
266  if ( mCellSizeX < 0 )
267  {
269  }
270  mCellSizeY = geotransform[5];
271  if ( mCellSizeY < 0 )
272  {
274  }
275 
276  const char *projection = GDALGetProjectionRef( inputDataset );
277  GDALSetProjection( outputDataset, projection );
278 
279  return outputDataset;
280 }
281 
#define QgsDebugMsg(str)
Definition: qgslogger.h:37
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
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:54
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.