QGIS API Documentation  3.4.15-Madeira (e83d02e274)
qgsrasteranalysisutils.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsrasteranalysisutils.cpp
3  ---------------------
4  Date : June 2018
5  Copyright : (C) 2018 by Nyall Dawson
6  Email : nyall dot dawson at gmail dot com
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include "qgsrasteranalysisutils.h"
17 
18 #include "qgsfeedback.h"
19 #include "qgsrasterblock.h"
20 #include "qgsrasteriterator.h"
21 #include "qgsgeos.h"
23 #include <map>
25 
26 void QgsRasterAnalysisUtils::cellInfoForBBox( const QgsRectangle &rasterBBox, const QgsRectangle &featureBBox, double cellSizeX, double cellSizeY,
27  int &nCellsX, int &nCellsY, int rasterWidth, int rasterHeight, QgsRectangle &rasterBlockExtent )
28 {
29  //get intersecting bbox
30  QgsRectangle intersectBox = rasterBBox.intersect( featureBBox );
31  if ( intersectBox.isEmpty() )
32  {
33  nCellsX = 0;
34  nCellsY = 0;
35  rasterBlockExtent = QgsRectangle();
36  return;
37  }
38 
39  //get offset in pixels in x- and y- direction
40  int offsetX = static_cast< int >( std::floor( ( intersectBox.xMinimum() - rasterBBox.xMinimum() ) / cellSizeX ) );
41  int offsetY = static_cast< int >( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMaximum() ) / cellSizeY ) );
42 
43  int maxColumn = static_cast< int >( std::floor( ( intersectBox.xMaximum() - rasterBBox.xMinimum() ) / cellSizeX ) ) + 1;
44  int maxRow = static_cast< int >( std::floor( ( rasterBBox.yMaximum() - intersectBox.yMinimum() ) / cellSizeY ) ) + 1;
45 
46  nCellsX = maxColumn - offsetX;
47  nCellsY = maxRow - offsetY;
48 
49  //avoid access to cells outside of the raster (may occur because of rounding)
50  nCellsX = std::min( offsetX + nCellsX, rasterWidth ) - offsetX;
51  nCellsY = std::min( offsetY + nCellsY, rasterHeight ) - offsetY;
52 
53  rasterBlockExtent = QgsRectangle( rasterBBox.xMinimum() + offsetX * cellSizeX,
54  rasterBBox.yMaximum() - offsetY * cellSizeY,
55  rasterBBox.xMinimum() + ( nCellsX + offsetX ) * cellSizeX,
56  rasterBBox.yMaximum() - ( nCellsY + offsetY ) * cellSizeY );
57 }
58 
59 void QgsRasterAnalysisUtils::statisticsFromMiddlePointTest( QgsRasterInterface *rasterInterface, int rasterBand, const QgsGeometry &poly, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, const std::function<void( double )> &addValue, bool skipNodata )
60 {
61  std::unique_ptr< QgsGeometryEngine > polyEngine( QgsGeometry::createGeometryEngine( poly.constGet( ) ) );
62  if ( !polyEngine )
63  {
64  return;
65  }
66  polyEngine->prepareGeometry();
67 
68  QgsRasterIterator iter( rasterInterface );
69  iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
70 
71  std::unique_ptr< QgsRasterBlock > block;
72  int iterLeft = 0;
73  int iterTop = 0;
74  int iterCols = 0;
75  int iterRows = 0;
76  QgsRectangle blockExtent;
77  while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
78  {
79  double cellCenterY = blockExtent.yMaximum() - 0.5 * cellSizeY;
80 
81  for ( int row = 0; row < iterRows; ++row )
82  {
83  double cellCenterX = blockExtent.xMinimum() + 0.5 * cellSizeX;
84  for ( int col = 0; col < iterCols; ++col )
85  {
86  double pixelValue = block->value( row, col );
87  if ( validPixel( pixelValue ) && ( !skipNodata || !block->isNoData( row, col ) ) )
88  {
89  QgsPoint cellCenter( cellCenterX, cellCenterY );
90  if ( polyEngine->contains( &cellCenter ) )
91  {
92  addValue( pixelValue );
93  }
94  }
95  cellCenterX += cellSizeX;
96  }
97  cellCenterY -= cellSizeY;
98  }
99  }
100 }
101 
102 void QgsRasterAnalysisUtils::statisticsFromPreciseIntersection( QgsRasterInterface *rasterInterface, int rasterBand, const QgsGeometry &poly, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, const std::function<void( double, double )> &addValue, bool skipNodata )
103 {
104  QgsGeometry pixelRectGeometry;
105 
106  double hCellSizeX = cellSizeX / 2.0;
107  double hCellSizeY = cellSizeY / 2.0;
108  double pixelArea = cellSizeX * cellSizeY;
109  double weight = 0;
110 
111  std::unique_ptr< QgsGeometryEngine > polyEngine( QgsGeometry::createGeometryEngine( poly.constGet( ) ) );
112  if ( !polyEngine )
113  {
114  return;
115  }
116  polyEngine->prepareGeometry();
117 
118  QgsRasterIterator iter( rasterInterface );
119  iter.startRasterRead( rasterBand, nCellsX, nCellsY, rasterBBox );
120 
121  std::unique_ptr< QgsRasterBlock > block;
122  int iterLeft = 0;
123  int iterTop = 0;
124  int iterCols = 0;
125  int iterRows = 0;
126  QgsRectangle blockExtent;
127  while ( iter.readNextRasterPart( rasterBand, iterCols, iterRows, block, iterLeft, iterTop, &blockExtent ) )
128  {
129  double currentY = blockExtent.yMaximum() - 0.5 * cellSizeY;
130  for ( int row = 0; row < iterRows; ++row )
131  {
132  double currentX = blockExtent.xMinimum() + 0.5 * cellSizeX;
133  for ( int col = 0; col < iterCols; ++col )
134  {
135  double pixelValue = block->value( row, col );
136  if ( validPixel( pixelValue ) && ( !skipNodata || !block->isNoData( row, col ) ) )
137  {
138  pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
139  // GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
140  // so we first test to see if there IS an intersection before doing the actual calculation
141  if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
142  {
143  //intersection
144  QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
145  if ( !intersectGeometry.isEmpty() )
146  {
147  double intersectionArea = intersectGeometry.area();
148  if ( intersectionArea > 0.0 )
149  {
150  weight = intersectionArea / pixelArea;
151  addValue( pixelValue, weight );
152  }
153  }
154  }
155  }
156  currentX += cellSizeX;
157  }
158  currentY -= cellSizeY;
159  }
160  }
161 }
162 
163 bool QgsRasterAnalysisUtils::validPixel( double value )
164 {
165  return !std::isnan( value );
166 }
167 
168 static QVector< QPair< QString, Qgis::DataType > > sDataTypes;
169 
170 void populateDataTypes()
171 {
172  if ( sDataTypes.empty() )
173  {
174  sDataTypes.append( qMakePair( QStringLiteral( "Byte" ), Qgis::Byte ) );
175  sDataTypes.append( qMakePair( QStringLiteral( "Int16" ), Qgis::Int16 ) );
176  sDataTypes.append( qMakePair( QStringLiteral( "UInt16" ), Qgis::UInt16 ) );
177  sDataTypes.append( qMakePair( QStringLiteral( "Int32" ), Qgis::Int32 ) );
178  sDataTypes.append( qMakePair( QStringLiteral( "UInt32" ), Qgis::UInt32 ) );
179  sDataTypes.append( qMakePair( QStringLiteral( "Float32" ), Qgis::Float32 ) );
180  sDataTypes.append( qMakePair( QStringLiteral( "Float64" ), Qgis::Float64 ) );
181  sDataTypes.append( qMakePair( QStringLiteral( "CInt16" ), Qgis::CInt16 ) );
182  sDataTypes.append( qMakePair( QStringLiteral( "CInt32" ), Qgis::CInt32 ) );
183  sDataTypes.append( qMakePair( QStringLiteral( "CFloat32" ), Qgis::CFloat32 ) );
184  sDataTypes.append( qMakePair( QStringLiteral( "CFloat64" ), Qgis::CFloat64 ) );
185  }
186 }
187 
188 std::unique_ptr<QgsProcessingParameterDefinition> QgsRasterAnalysisUtils::createRasterTypeParameter( const QString &name, const QString &description, Qgis::DataType defaultType )
189 {
190  populateDataTypes();
191 
192  QStringList names;
193  int defaultChoice = 0;
194  int i = 0;
195  for ( auto it = sDataTypes.constBegin(); it != sDataTypes.constEnd(); ++it )
196  {
197  names.append( it->first );
198  if ( it->second == defaultType )
199  defaultChoice = i;
200  i++;
201  }
202 
203  return qgis::make_unique< QgsProcessingParameterEnum >( name, description, names, false, defaultChoice );
204 }
205 
206 Qgis::DataType QgsRasterAnalysisUtils::rasterTypeChoiceToDataType( int choice )
207 {
208  if ( choice < 0 || choice >= sDataTypes.count() )
209  return Qgis::Float32;
210 
211  return sDataTypes.value( choice ).second;
212 }
213 
215 
A rectangle specified with double values.
Definition: qgsrectangle.h:40
Thirty two bit signed integer (qint32)
Definition: qgis.h:99
bool isEmpty() const
Returns true if the rectangle is empty.
Definition: qgsrectangle.h:425
Iterator for sequentially processing raster cells.
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
double yMaximum() const
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:171
Thirty two bit unsigned integer (quint32)
Definition: qgis.h:98
DataType
Raster data types.
Definition: qgis.h:92
QgsGeometry intersection(const QgsGeometry &geometry) const
Returns a geometry representing the points shared by this geometry and other.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:106
Thirty two bit floating point (float)
Definition: qgis.h:100
Sixteen bit signed integer (qint16)
Definition: qgis.h:97
Complex Int16.
Definition: qgis.h:102
Sixty four bit floating point (double)
Definition: qgis.h:101
static QgsGeometry fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:176
double xMaximum() const
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:161
Complex Float32.
Definition: qgis.h:104
Complex Int32.
Definition: qgis.h:103
Sixteen bit unsigned integer (quint16)
Definition: qgis.h:96
Base class for processing filters like renderers, reprojector, resampler etc.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
QgsRectangle intersect(const QgsRectangle &rect) const
Returns the intersection with the given rectangle.
Definition: qgsrectangle.h:311
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry)
Creates and returns a new geometry engine.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
Complex Float64.
Definition: qgis.h:105
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:166
double area() const
Returns the area of the geometry using GEOS.
Eight bit unsigned integer (quint8)
Definition: qgis.h:95