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