QGIS API Documentation  3.27.0-Master (bef583a8ef)
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 "qgsgdalutils.h"
19 #include "qgsninecellfilter.h"
20 #include "qgslogger.h"
21 #include "cpl_string.h"
22 #include "qgsfeedback.h"
23 #include "qgsogrutils.h"
24 #include "qgsmessagelog.h"
25 #include "qgsconfig.h"
26 
27 #ifdef HAVE_OPENCL
28 #include "qgsopenclutils.h"
29 #endif
30 
31 #include <QFile>
32 #include <QDebug>
33 #include <QFileInfo>
34 #include <iterator>
35 
36 
37 
38 QgsNineCellFilter::QgsNineCellFilter( const QString &inputFile, const QString &outputFile, const QString &outputFormat )
39  : mInputFile( inputFile )
40  , mOutputFile( outputFile )
41  , mOutputFormat( outputFormat )
42 {
43 
44 }
45 
46 // TODO: return an anum instead of an int
48 {
49 #ifdef HAVE_OPENCL
50  if ( QgsOpenClUtils::enabled() && QgsOpenClUtils::available() && ! openClProgramBaseName( ).isEmpty() )
51  {
52  // Load the program sources
53  const QString source( QgsOpenClUtils::sourceFromBaseName( openClProgramBaseName( ) ) );
54  if ( ! source.isEmpty() )
55  {
56  try
57  {
58  QgsDebugMsg( QObject::tr( "Running OpenCL program: %1" ).arg( openClProgramBaseName( ) ) );
59  return processRasterGPU( source, feedback );
60  }
61  catch ( cl::Error &e )
62  {
63  const QString err = QObject::tr( "Error running OpenCL program: %1 - %2" ).arg( e.what( ), QgsOpenClUtils::errorText( e.err( ) ) );
64  QgsMessageLog::logMessage( err, QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::MessageLevel::Critical );
65  throw QgsProcessingException( err );
66  }
67  }
68  else
69  {
70  const QString err = QObject::tr( "Error loading OpenCL program sources" );
72  QgsOpenClUtils::LOGMESSAGE_TAG, Qgis::MessageLevel::Critical );
73  throw QgsProcessingException( err );
74  }
75  }
76  else
77  {
78  return processRasterCPU( feedback );
79  }
80 #ifndef _MSC_VER
81  return 1;
82 #endif
83 #else
84  return processRasterCPU( feedback );
85 #endif
86 }
87 
88 gdal::dataset_unique_ptr QgsNineCellFilter::openInputFile( int &nCellsX, int &nCellsY )
89 {
90  gdal::dataset_unique_ptr inputDataset( GDALOpen( mInputFile.toUtf8().constData(), GA_ReadOnly ) );
91  if ( inputDataset )
92  {
93  nCellsX = GDALGetRasterXSize( inputDataset.get() );
94  nCellsY = GDALGetRasterYSize( inputDataset.get() );
95 
96  //we need at least one band
97  if ( GDALGetRasterCount( inputDataset.get() ) < 1 )
98  {
99  return nullptr;
100  }
101  }
102  return inputDataset;
103 }
104 
105 GDALDriverH QgsNineCellFilter::openOutputDriver()
106 {
107  //open driver
108  GDALDriverH outputDriver = GDALGetDriverByName( mOutputFormat.toLocal8Bit().data() );
109 
110  if ( !outputDriver )
111  {
112  return outputDriver; //return nullptr, driver does not exist
113  }
114 
115  if ( !QgsGdalUtils::supportsRasterCreate( outputDriver ) )
116  {
117  return nullptr; //driver exist, but it does not support the create operation
118  }
119 
120  return outputDriver;
121 }
122 
123 
124 gdal::dataset_unique_ptr QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver )
125 {
126  if ( !inputDataset )
127  {
128  return nullptr;
129  }
130 
131  const int xSize = GDALGetRasterXSize( inputDataset );
132  const int ySize = GDALGetRasterYSize( inputDataset );
133 
134  //open output file
135  char **papszOptions = nullptr;
136  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 1, GDT_Float32, papszOptions ) );
137  if ( !outputDataset )
138  {
139  return outputDataset;
140  }
141 
142  //get geotransform from inputDataset
143  double geotransform[6];
144  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )
145  {
146  return nullptr;
147  }
148  GDALSetGeoTransform( outputDataset.get(), geotransform );
149 
150  //make sure mCellSizeX and mCellSizeY are always > 0
151  mCellSizeX = geotransform[1];
152  if ( mCellSizeX < 0 )
153  {
155  }
156  mCellSizeY = geotransform[5];
157  if ( mCellSizeY < 0 )
158  {
160  }
161 
162  const char *projection = GDALGetProjectionRef( inputDataset );
163  GDALSetProjection( outputDataset.get(), projection );
164 
165  return outputDataset;
166 }
167 
168 #ifdef HAVE_OPENCL
169 
170 // TODO: return an anum instead of an int
171 int QgsNineCellFilter::processRasterGPU( const QString &source, QgsFeedback *feedback )
172 {
173 
174  GDALAllRegister();
175 
176  //open input file
177  int xSize, ySize;
178  const gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) );
179  if ( !inputDataset )
180  {
181  return 1; //opening of input file failed
182  }
183 
184  //output driver
185  GDALDriverH outputDriver = openOutputDriver();
186  if ( !outputDriver )
187  {
188  return 2;
189  }
190 
191  gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) );
192  if ( !outputDataset )
193  {
194  return 3; //create operation on output file failed
195  }
196 
197  //open first raster band for reading (operation is only for single band raster)
198  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
199  if ( !rasterBand )
200  {
201  return 4;
202  }
203  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
204 
205  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
206  if ( !outputRasterBand )
207  {
208  return 5;
209  }
210  //try to set -9999 as nodata value
211  GDALSetRasterNoDataValue( outputRasterBand, -9999 );
212  mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
213 
214  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
215  {
216  return 6;
217  }
218 
219  // Prepare context and queue
220  const cl::Context ctx = QgsOpenClUtils::context();
221  cl::CommandQueue queue = QgsOpenClUtils::commandQueue();
222 
223  //keep only three scanlines in memory at a time, make room for initial and final nodata
224  QgsOpenClUtils::CPLAllocator<float> scanLine( xSize + 2 );
225  QgsOpenClUtils::CPLAllocator<float> resultLine( xSize );
226 
227  // Cast to float (because double just crashes on some GPUs)
228  std::vector<float> rasterParams;
229 
230  rasterParams.push_back( mInputNodataValue ); // 0
231  rasterParams.push_back( mOutputNodataValue ); // 1
232  rasterParams.push_back( mZFactor ); // 2
233  rasterParams.push_back( mCellSizeX ); // 3
234  rasterParams.push_back( mCellSizeY ); // 4
235 
236  // Allow subclasses to add extra params needed for computation:
237  // used to pass additional args to opencl program
238  addExtraRasterParams( rasterParams );
239 
240  const std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) );
241  const std::size_t inputSize( sizeof( float ) * ( xSize ) );
242 
243  cl::Buffer rasterParamsBuffer( queue, rasterParams.begin(), rasterParams.end(), true, false, nullptr );
244  cl::Buffer scanLine1Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
245  cl::Buffer scanLine2Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
246  cl::Buffer scanLine3Buffer( ctx, CL_MEM_READ_ONLY, bufferSize, nullptr, nullptr );
247  cl::Buffer *scanLineBuffer[3] = {&scanLine1Buffer, &scanLine2Buffer, &scanLine3Buffer};
248  cl::Buffer resultLineBuffer( ctx, CL_MEM_WRITE_ONLY, inputSize, nullptr, nullptr );
249 
250  // Create a program from the kernel source
251  const cl::Program program( QgsOpenClUtils::buildProgram( source, QgsOpenClUtils::ExceptionBehavior::Throw ) );
252 
253  // Create the OpenCL kernel
254  auto kernel = cl::KernelFunctor <
255  cl::Buffer &,
256  cl::Buffer &,
257  cl::Buffer &,
258  cl::Buffer &,
259  cl::Buffer &
260  > ( program, "processNineCellWindow" );
261 
262  // Rotate buffer index
263  std::vector<int> rowIndex = {0, 1, 2};
264 
265  // values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
266  for ( int i = 0; i < ySize; ++i )
267  {
268  if ( feedback && feedback->isCanceled() )
269  {
270  break;
271  }
272 
273  if ( feedback )
274  {
275  feedback->setProgress( 100.0 * static_cast< double >( i ) / ySize );
276  }
277 
278  if ( i == 0 )
279  {
280  // Fill scanline 1 with (input) nodata for the values above the first row and
281  // feed scanline2 with the first actual data row
282  for ( int a = 0; a < xSize + 2 ; ++a )
283  {
284  scanLine[a] = mInputNodataValue;
285  }
286  queue.enqueueWriteBuffer( scanLine1Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
287 
288  // Read scanline2: first real raster row
289  if ( GDALRasterIO( rasterBand, GF_Read, 0, i, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
290  {
291  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
292  }
293  queue.enqueueWriteBuffer( scanLine2Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
294 
295  // Read scanline3: second real raster row
296  if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
297  {
298  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
299  }
300  queue.enqueueWriteBuffer( scanLine3Buffer, CL_TRUE, 0, bufferSize, scanLine.get() );
301  }
302  else
303  {
304  // Normally fetch only scanLine3 and move forward one row
305  // Read scanline 3, fill the last row with nodata values if it's the last iteration
306  if ( i == ySize - 1 ) //fill the row below the bottom with nodata values
307  {
308  for ( int a = 0; a < xSize + 2; ++a )
309  {
310  scanLine[a] = mInputNodataValue;
311  }
312  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
313  }
314  else // Read line i + 1 and put it into scanline 3
315  // Overwrite from input, skip first and last
316  {
317  if ( GDALRasterIO( rasterBand, GF_Read, 0, i + 1, xSize, 1, &scanLine[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
318  {
319  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
320  }
321  queue.enqueueWriteBuffer( *scanLineBuffer[rowIndex[2]], CL_TRUE, 0, bufferSize, scanLine.get() ); // row 0
322  }
323  }
324 
325  kernel( cl::EnqueueArgs(
326  queue,
327  cl::NDRange( xSize )
328  ),
329  *scanLineBuffer[rowIndex[0]],
330  *scanLineBuffer[rowIndex[1]],
331  *scanLineBuffer[rowIndex[2]],
332  resultLineBuffer,
333  rasterParamsBuffer
334  );
335 
336  queue.enqueueReadBuffer( resultLineBuffer, CL_TRUE, 0, inputSize, resultLine.get() );
337 
338  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, i, xSize, 1, resultLine.get(), xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
339  {
340  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
341  }
342  std::rotate( rowIndex.begin(), rowIndex.begin() + 1, rowIndex.end() );
343  }
344 
345  if ( feedback && feedback->isCanceled() )
346  {
347  //delete the dataset without closing (because it is faster)
348  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
349  return 7;
350  }
351  return 0;
352 }
353 #endif
354 
355 
356 // TODO: return an anum instead of an int
357 int QgsNineCellFilter::processRasterCPU( QgsFeedback *feedback )
358 {
359 
360  GDALAllRegister();
361 
362  //open input file
363  int xSize, ySize;
364  const gdal::dataset_unique_ptr inputDataset( openInputFile( xSize, ySize ) );
365  if ( !inputDataset )
366  {
367  return 1; //opening of input file failed
368  }
369 
370  //output driver
371  GDALDriverH outputDriver = openOutputDriver();
372  if ( !outputDriver )
373  {
374  return 2;
375  }
376 
377  gdal::dataset_unique_ptr outputDataset( openOutputFile( inputDataset.get(), outputDriver ) );
378  if ( !outputDataset )
379  {
380  return 3; //create operation on output file failed
381  }
382 
383  //open first raster band for reading (operation is only for single band raster)
384  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset.get(), 1 );
385  if ( !rasterBand )
386  {
387  return 4;
388  }
389  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, nullptr );
390 
391  GDALRasterBandH outputRasterBand = GDALGetRasterBand( outputDataset.get(), 1 );
392  if ( !outputRasterBand )
393  {
394  return 5;
395  }
396  //try to set -9999 as nodata value
397  GDALSetRasterNoDataValue( outputRasterBand, -9999 );
398  mOutputNodataValue = GDALGetRasterNoDataValue( outputRasterBand, nullptr );
399 
400  if ( ySize < 3 ) //we require at least three rows (should be true for most datasets)
401  {
402  return 6;
403  }
404 
405  //keep only three scanlines in memory at a time, make room for initial and final nodata
406  const std::size_t bufferSize( sizeof( float ) * ( xSize + 2 ) );
407  float *scanLine1 = ( float * ) CPLMalloc( bufferSize );
408  float *scanLine2 = ( float * ) CPLMalloc( bufferSize );
409  float *scanLine3 = ( float * ) CPLMalloc( bufferSize );
410 
411  float *resultLine = ( float * ) CPLMalloc( sizeof( float ) * xSize );
412 
413  //values outside the layer extent (if the 3x3 window is on the border) are sent to the processing method as (input) nodata values
414  for ( int yIndex = 0; yIndex < ySize; ++yIndex )
415  {
416  if ( feedback && feedback->isCanceled() )
417  {
418  break;
419  }
420 
421  if ( feedback )
422  {
423  feedback->setProgress( 100.0 * static_cast< double >( yIndex ) / ySize );
424  }
425 
426  if ( yIndex == 0 )
427  {
428  //fill scanline 1 with (input) nodata for the values above the first row and feed scanline2 with the first row
429  for ( int a = 0; a < xSize + 2 ; ++a )
430  {
431  scanLine1[a] = mInputNodataValue;
432  }
433  // Read scanline2
434  if ( GDALRasterIO( rasterBand, GF_Read, 0, 0, xSize, 1, &scanLine2[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
435  {
436  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
437  }
438  }
439  else
440  {
441  //normally fetch only scanLine3 and release scanline 1 if we move forward one row
442  CPLFree( scanLine1 );
443  scanLine1 = scanLine2;
444  scanLine2 = scanLine3;
445  scanLine3 = ( float * ) CPLMalloc( bufferSize );
446  }
447 
448  // Read scanline 3
449  if ( yIndex == ySize - 1 ) //fill the row below the bottom with nodata values
450  {
451  for ( int a = 0; a < xSize + 2; ++a )
452  {
453  scanLine3[a] = mInputNodataValue;
454  }
455  }
456  else
457  {
458  if ( GDALRasterIO( rasterBand, GF_Read, 0, yIndex + 1, xSize, 1, &scanLine3[1], xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
459  {
460  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
461  }
462  }
463  // Set first and last extra columns to nodata
464  scanLine1[0] = scanLine1[xSize + 1] = mInputNodataValue;
465  scanLine2[0] = scanLine2[xSize + 1] = mInputNodataValue;
466  scanLine3[0] = scanLine3[xSize + 1] = mInputNodataValue;
467 
468 
469 
470  // j is the x axis index, skip 0 and last cell that have been filled with nodata
471  for ( int xIndex = 0; xIndex < xSize ; ++xIndex )
472  {
473  // cells(x, y) x11, x21, x31, x12, x22, x32, x13, x23, x33
474  resultLine[ xIndex ] = processNineCellWindow( &scanLine1[ xIndex ], &scanLine1[ xIndex + 1 ], &scanLine1[ xIndex + 2 ],
475  &scanLine2[ xIndex ], &scanLine2[ xIndex + 1 ], &scanLine2[ xIndex + 2 ],
476  &scanLine3[ xIndex ], &scanLine3[ xIndex + 1 ], &scanLine3[ xIndex + 2 ] );
477 
478  }
479 
480  if ( GDALRasterIO( outputRasterBand, GF_Write, 0, yIndex, xSize, 1, resultLine, xSize, 1, GDT_Float32, 0, 0 ) != CE_None )
481  {
482  QgsDebugMsg( QStringLiteral( "Raster IO Error" ) );
483  }
484  }
485 
486  CPLFree( resultLine );
487  CPLFree( scanLine1 );
488  CPLFree( scanLine2 );
489  CPLFree( scanLine3 );
490 
491  if ( feedback && feedback->isCanceled() )
492  {
493  //delete the dataset without closing (because it is faster)
494  gdal::fast_delete_and_close( outputDataset, outputDriver, mOutputFile );
495  return 7;
496  }
497  return 0;
498 }
Base class for feedback objects to be used for cancellation of something running in a worker thread.
Definition: qgsfeedback.h:45
bool isCanceled() const SIP_HOLDGIL
Tells whether the operation has been canceled already.
Definition: qgsfeedback.h:54
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition: qgsfeedback.h:63
static bool supportsRasterCreate(GDALDriverH driver)
Reads whether a driver supports GDALCreate() for raster purposes.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
int processRaster(QgsFeedback *feedback=nullptr)
Starts the calculation, reads from mInputFile and stores the result in mOutputFile.
float mInputNodataValue
The nodata value of the input layer.
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.
QgsNineCellFilter(const QString &inputFile, const QString &outputFile, const QString &outputFormat)
Constructor that takes input file, output file and output format (GDAL string)
float mOutputNodataValue
The nodata value of the output layer.
double mZFactor
Scale factor for z-value if x-/y- units are different to z-units (111120 for degree->meters and 37040...
static Q_DECL_DEPRECATED cl::Program buildProgram(const cl::Context &context, const QString &source, ExceptionBehavior exceptionBehavior=Catch)
Build the program from source in the given context and depending on exceptionBehavior can throw or ca...
static cl::Context context()
Context factory.
static bool enabled()
Returns true if OpenCL is enabled in the user settings.
static bool available()
Checks whether a suitable OpenCL platform and device is available on this system and initialize the Q...
static QString errorText(const int errorCode)
Returns a string representation from an OpenCL errorCode.
static cl::CommandQueue commandQueue()
Create an OpenCL command queue from the default context.
static QString sourceFromBaseName(const QString &baseName)
Returns the full path to a an OpenCL source file from the baseName ('.cl' extension is automatically ...
static QLatin1String LOGMESSAGE_TAG
OpenCL string for message logs.
Custom exception class for processing related exceptions.
Definition: qgsexception.h:83
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:80
std::unique_ptr< std::remove_pointer< GDALDatasetH >::type, GDALDatasetCloser > dataset_unique_ptr
Scoped GDAL dataset.
Definition: qgsogrutils.h:139
void * GDALDatasetH
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
Tiny smart-pointer-like wrapper around CPLMalloc and CPLFree: this is needed because OpenCL C++ API m...