QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
qgsalgorithmlinedensity.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmlinedensity.cpp
3 ---------------------
4 begin : December 2019
5 copyright : (C) 2019 by Clemens Raffler
6 email : clemens dot raffler at gmail dot com
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
19#include "qgscircle.h"
20#include "qgsgeometryengine.h"
21#include "qgsrasterfilewriter.h"
22
24
25QString QgsLineDensityAlgorithm::name() const
26{
27 return QStringLiteral( "linedensity" );
28}
29
30QString QgsLineDensityAlgorithm::displayName() const
31{
32 return QObject::tr( "Line density" );
33}
34
35QStringList QgsLineDensityAlgorithm::tags() const
36{
37 return QObject::tr( "density,kernel,line,line density,interpolation,weight" ).split( ',' );
38}
39
40QString QgsLineDensityAlgorithm::group() const
41{
42 return QObject::tr( "Interpolation" );
43}
44
45QString QgsLineDensityAlgorithm::groupId() const
46{
47 return QStringLiteral( "interpolation" );
48}
49
50void QgsLineDensityAlgorithm::initAlgorithm( const QVariantMap & )
51{
52 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input line layer" ), QList<int>() << QgsProcessing::TypeVectorLine ) );
53 addParameter( new QgsProcessingParameterField( QStringLiteral( "WEIGHT" ), QObject::tr( "Weight field " ), QVariant(), QStringLiteral( "INPUT" ), QgsProcessingParameterField::Numeric, false, true ) );
54 addParameter( new QgsProcessingParameterDistance( QStringLiteral( "RADIUS" ), QObject::tr( "Search radius" ), 10, QStringLiteral( "INPUT" ), false, 0 ) );
55 addParameter( new QgsProcessingParameterDistance( QStringLiteral( "PIXEL_SIZE" ), QObject::tr( "Pixel size" ), 10, QStringLiteral( "INPUT" ), false ) );
56
57 addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Line density raster" ) ) );
58}
59
60QString QgsLineDensityAlgorithm::shortHelpString() const
61{
62 return QObject::tr( "The line density interpolation algorithm calculates a density measure of linear features "
63 "which is obtained in a circular neighborhood within each raster cell. "
64 "First, the length of the segment of each line that is intersected by the circular neighborhood "
65 "is multiplied with the lines weight factor. In a second step, all length values are summed and "
66 "divided by the area of the circular neighborhood. This process is repeated for all raster cells."
67 );
68}
69
70QgsLineDensityAlgorithm *QgsLineDensityAlgorithm::createInstance() const
71{
72 return new QgsLineDensityAlgorithm();
73}
74
75bool QgsLineDensityAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
76{
77 Q_UNUSED( feedback );
78 mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
79 if ( !mSource )
80 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
81
82 mWeightField = parameterAsString( parameters, QStringLiteral( "WEIGHT" ), context );
83
84 mPixelSize = parameterAsDouble( parameters, QStringLiteral( "PIXEL_SIZE" ), context );
85
86 mSearchRadius = parameterAsDouble( parameters, QStringLiteral( "RADIUS" ), context );
87 if ( mSearchRadius < 0.5 * mPixelSize * std::sqrt( 2 ) )
88 throw QgsProcessingException( QObject::tr( "Raster cells must be fully contained by the search circle. Therefore, "
89 "the search radius must not be smaller than half of the pixel diagonal." ) );
90
91 mExtent = mSource->sourceExtent();
92 mCrs = mSource->sourceCrs();
93 mDa = QgsDistanceArea();
94 mDa.setEllipsoid( context.ellipsoid() );
95 mDa.setSourceCrs( mCrs, context.transformContext() );
96
97 //get cell midpoint from top left cell
98 const QgsPoint firstCellMidpoint = QgsPoint( mExtent.xMinimum() + ( mPixelSize / 2 ), mExtent.yMaximum() - ( mPixelSize / 2 ) );
99 const QgsCircle searchCircle = QgsCircle( firstCellMidpoint, mSearchRadius );
100 mSearchGeometry = QgsGeometry( searchCircle.toPolygon() );
101
102 return true;
103}
104
105QVariantMap QgsLineDensityAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
106{
108
109 const QStringList weightName = QStringList( mWeightField );
110 const QgsFields attrFields = mSource->fields();
111
113 r.setSubsetOfAttributes( weightName, attrFields );
114 QgsFeatureIterator fit = mSource->getFeatures( r );
115 QgsFeature f;
116
117 while ( fit.nextFeature( f ) )
118 {
119 mIndex.addFeature( f, QgsFeatureSink::FastInsert );
120
121 //only populate hash if weight field is given
122 if ( !mWeightField.isEmpty() )
123 {
124 const double analysisWeight = f.attribute( mWeightField ).toDouble();
125 mFeatureWeights.insert( f.id(), analysisWeight );
126 }
127 }
128
129 const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
130 const QFileInfo fi( outputFile );
131 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
132
133 const int rows = std::max( std::ceil( mExtent.height() / mPixelSize ), 1.0 );
134 const int cols = std::max( std::ceil( mExtent.width() / mPixelSize ), 1.0 );
135
136 //build new raster extent based on number of columns and cellsize
137 //this prevents output cellsize being calculated too small
138 const QgsRectangle rasterExtent = QgsRectangle( mExtent.xMinimum(), mExtent.yMaximum() - ( rows * mPixelSize ), mExtent.xMinimum() + ( cols * mPixelSize ), mExtent.yMaximum() );
139
140 QgsRasterFileWriter writer = QgsRasterFileWriter( outputFile );
141 writer.setOutputProviderKey( QStringLiteral( "gdal" ) );
142 writer.setOutputFormat( outputFormat );
143 std::unique_ptr<QgsRasterDataProvider > provider( writer.createOneBandRaster( Qgis::DataType::Float32, cols, rows, rasterExtent, mCrs ) );
144 if ( !provider )
145 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
146 if ( !provider->isValid() )
147 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
148
149 provider->setNoDataValue( 1, -9999 );
150
151 const qgssize totalCellcnt = static_cast<qgssize>( rows ) * cols;
152 int cellcnt = 0;
153
154 std::unique_ptr< QgsRasterBlock > rasterDataLine = std::make_unique< QgsRasterBlock >( Qgis::DataType::Float32, cols, 1 );
155
156 for ( int row = 0; row < rows; row++ )
157 {
158 for ( int col = 0; col < cols; col++ )
159 {
160 if ( feedback->isCanceled() )
161 {
162 break;
163 }
164
165 if ( col > 0 )
166 mSearchGeometry.translate( mPixelSize, 0 );
167
168 const QList<QgsFeatureId> fids = mIndex.intersects( mSearchGeometry.boundingBox() );
169
170 if ( !fids.isEmpty() )
171 {
172 std::unique_ptr< QgsGeometryEngine > engine( QgsGeometry::createGeometryEngine( mSearchGeometry.constGet() ) );
173 engine->prepareGeometry();
174
175 double absDensity = 0;
176 for ( const QgsFeatureId id : fids )
177 {
178 const QgsGeometry lineGeom = mIndex.geometry( id );
179
180 if ( engine->intersects( lineGeom.constGet() ) )
181 {
182 const double analysisLineLength = mDa.measureLength( QgsGeometry( engine->intersection( mIndex.geometry( id ).constGet() ) ) );
183 double weight = 1;
184
185 if ( !mWeightField.isEmpty() )
186 {
187 weight = mFeatureWeights.value( id );
188 }
189
190 absDensity += ( analysisLineLength * weight );
191 }
192 }
193
194 double lineDensity = 0;
195 if ( absDensity > 0 )
196 {
197 //only calculate ellipsoidal area if abs density is greater 0
198 const double analysisSearchGeometryArea = mDa.measureArea( mSearchGeometry );
199 lineDensity = absDensity / analysisSearchGeometryArea;
200 }
201 rasterDataLine->setValue( 0, col, lineDensity );
202 }
203 else
204 {
205 //no lines found in search radius
206 rasterDataLine->setValue( 0, col, 0.0 );
207 }
208
209 feedback->setProgress( static_cast<double>( cellcnt ) / static_cast<double>( totalCellcnt ) * 100 );
210 cellcnt++;
211 }
212 provider->writeBlock( rasterDataLine.get(), 1, 0, row );
213
214 //'carriage return and newline' for search geometry
215 mSearchGeometry.translate( ( cols - 1 ) * -mPixelSize, -mPixelSize );
216 }
217
218 QVariantMap outputs;
219 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
220 return outputs;
221}
222
223
@ Float32
Thirty two bit floating point (float)
Circle geometry type.
Definition: qgscircle.h:44
A general purpose distance and area calculator, capable of performing ellipsoid based calculations.
virtual QgsPolygon * toPolygon(unsigned int segments=36) const
Returns a segmented polygon.
Definition: qgsellipse.cpp:224
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
This class wraps a request for features to a vector layer (or directly its vector data provider).
QgsFeatureRequest & setSubsetOfAttributes(const QgsAttributeList &attrs)
Set a subset of attributes that will be fetched.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition: qgsfeature.h:56
QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Definition: qgsfeature.cpp:338
Q_GADGET QgsFeatureId id
Definition: qgsfeature.h:64
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
Container of fields for a vector layer.
Definition: qgsfields.h:45
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:164
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry)
Creates and returns a new geometry engine representing the specified geometry.
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
Contains information about the context in which a processing algorithm is executed.
QgsCoordinateTransformContext transformContext() const
Returns the coordinate transform context.
QString ellipsoid() const
Returns the ellipsoid to use for distance and area calculations.
Custom exception class for processing related exceptions.
Definition: qgsexception.h:83
Base class for providing feedback from a processing algorithm.
A double numeric parameter for distance values.
An input feature source (such as vector layers) parameter for processing algorithms.
A vector layer or feature source field parameter for processing algorithms.
@ Numeric
Accepts numeric fields.
A raster layer destination parameter, for specifying the destination path for a raster layer created ...
@ TypeVectorLine
Vector line layers.
Definition: qgsprocessing.h:50
The raster file writer which allows you to save a raster to a new file.
static QString driverForExtension(const QString &extension)
Returns the GDAL driver name for a specified file extension.
void setOutputProviderKey(const QString &key)
QgsRasterDataProvider * createOneBandRaster(Qgis::DataType dataType, int width, int height, const QgsRectangle &extent, const QgsCoordinateReferenceSystem &crs)
Create a raster file with one band without initializing the pixel data.
void setOutputFormat(const QString &format)
A rectangle specified with double values.
Definition: qgsrectangle.h:42
A spatial index for QgsFeature objects.
@ FlagStoreFeatureGeometries
Indicates that the spatial index should also store feature geometries. This requires more memory,...
unsigned long long qgssize
Qgssize is used instead of size_t, because size_t is stdlib type, unknown by SIP, and it would be har...
Definition: qgis.h:3032
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features
Definition: qgsfeatureid.h:28