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