QGIS API Documentation 3.40.0-Bratislava (b56115d8743)
Loading...
Searching...
No Matches
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>() << static_cast< int >( Qgis::ProcessingSourceType::VectorLine ) ) );
53 addParameter( new QgsProcessingParameterField( QStringLiteral( "WEIGHT" ), QObject::tr( "Weight field " ), QVariant(), QStringLiteral( "INPUT" ), Qgis::ProcessingFieldParameterDataType::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 std::unique_ptr< QgsProcessingParameterString > createOptsParam = std::make_unique< QgsProcessingParameterString >( QStringLiteral( "CREATE_OPTIONS" ), QObject::tr( "Creation options" ), QVariant(), false, true );
58 createOptsParam->setMetadata( QVariantMap( {{QStringLiteral( "widget_wrapper" ), QVariantMap( {{QStringLiteral( "widget_type" ), QStringLiteral( "rasteroptions" ) }} ) }} ) );
59 createOptsParam->setFlags( createOptsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
60 addParameter( createOptsParam.release() );
61
62 addParameter( new QgsProcessingParameterRasterDestination( QStringLiteral( "OUTPUT" ), QObject::tr( "Line density raster" ) ) );
63}
64
65QString QgsLineDensityAlgorithm::shortHelpString() const
66{
67 return QObject::tr( "The line density interpolation algorithm calculates a density measure of linear features "
68 "which is obtained in a circular neighborhood within each raster cell. "
69 "First, the length of the segment of each line that is intersected by the circular neighborhood "
70 "is multiplied with the lines weight factor. In a second step, all length values are summed and "
71 "divided by the area of the circular neighborhood. This process is repeated for all raster cells."
72 );
73}
74
75QgsLineDensityAlgorithm *QgsLineDensityAlgorithm::createInstance() const
76{
77 return new QgsLineDensityAlgorithm();
78}
79
80bool QgsLineDensityAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
81{
82 Q_UNUSED( feedback );
83 mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
84 if ( !mSource )
85 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
86
87 mWeightField = parameterAsString( parameters, QStringLiteral( "WEIGHT" ), context );
88
89 mPixelSize = parameterAsDouble( parameters, QStringLiteral( "PIXEL_SIZE" ), context );
90
91 mSearchRadius = parameterAsDouble( parameters, QStringLiteral( "RADIUS" ), context );
92 if ( mSearchRadius < 0.5 * mPixelSize * std::sqrt( 2 ) )
93 throw QgsProcessingException( QObject::tr( "Raster cells must be fully contained by the search circle. Therefore, "
94 "the search radius must not be smaller than half of the pixel diagonal." ) );
95
96 mExtent = mSource->sourceExtent();
97 mCrs = mSource->sourceCrs();
98 mDa = QgsDistanceArea();
99 mDa.setEllipsoid( context.ellipsoid() );
100 mDa.setSourceCrs( mCrs, context.transformContext() );
101
102 //get cell midpoint from top left cell
103 const QgsPoint firstCellMidpoint = QgsPoint( mExtent.xMinimum() + ( mPixelSize / 2 ), mExtent.yMaximum() - ( mPixelSize / 2 ) );
104 const QgsCircle searchCircle = QgsCircle( firstCellMidpoint, mSearchRadius );
105 mSearchGeometry = QgsGeometry( searchCircle.toPolygon() );
106
107 return true;
108}
109
110QVariantMap QgsLineDensityAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
111{
113
114 const QStringList weightName = QStringList( mWeightField );
115 const QgsFields attrFields = mSource->fields();
116
118 r.setSubsetOfAttributes( weightName, attrFields );
119 QgsFeatureIterator fit = mSource->getFeatures( r );
120 QgsFeature f;
121
122 while ( fit.nextFeature( f ) )
123 {
124 mIndex.addFeature( f, QgsFeatureSink::FastInsert );
125
126 //only populate hash if weight field is given
127 if ( !mWeightField.isEmpty() )
128 {
129 const double analysisWeight = f.attribute( mWeightField ).toDouble();
130 mFeatureWeights.insert( f.id(), analysisWeight );
131 }
132 }
133
134 const QString createOptions = parameterAsString( parameters, QStringLiteral( "CREATE_OPTIONS" ), context ).trimmed();
135 const QString outputFile = parameterAsOutputLayer( parameters, QStringLiteral( "OUTPUT" ), context );
136 const QFileInfo fi( outputFile );
137 const QString outputFormat = QgsRasterFileWriter::driverForExtension( fi.suffix() );
138
139 const int rows = std::max( std::ceil( mExtent.height() / mPixelSize ), 1.0 );
140 const int cols = std::max( std::ceil( mExtent.width() / mPixelSize ), 1.0 );
141
142 //build new raster extent based on number of columns and cellsize
143 //this prevents output cellsize being calculated too small
144 const QgsRectangle rasterExtent = QgsRectangle( mExtent.xMinimum(), mExtent.yMaximum() - ( rows * mPixelSize ), mExtent.xMinimum() + ( cols * mPixelSize ), mExtent.yMaximum() );
145
146 QgsRasterFileWriter writer = QgsRasterFileWriter( outputFile );
147 writer.setOutputProviderKey( QStringLiteral( "gdal" ) );
148 writer.setOutputFormat( outputFormat );
149 if ( !createOptions.isEmpty() )
150 {
151 writer.setCreateOptions( createOptions.split( '|' ) );
152 }
153
154 std::unique_ptr<QgsRasterDataProvider > provider( writer.createOneBandRaster( Qgis::DataType::Float32, cols, rows, rasterExtent, mCrs ) );
155 if ( !provider )
156 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
157 if ( !provider->isValid() )
158 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
159
160 provider->setNoDataValue( 1, -9999 );
161
162 const qgssize totalCellcnt = static_cast<qgssize>( rows ) * cols;
163 int cellcnt = 0;
164
165 std::unique_ptr< QgsRasterBlock > rasterDataLine = std::make_unique< QgsRasterBlock >( Qgis::DataType::Float32, cols, 1 );
166
167 for ( int row = 0; row < rows; row++ )
168 {
169 for ( int col = 0; col < cols; col++ )
170 {
171 if ( feedback->isCanceled() )
172 {
173 break;
174 }
175
176 if ( col > 0 )
177 mSearchGeometry.translate( mPixelSize, 0 );
178
179 const QList<QgsFeatureId> fids = mIndex.intersects( mSearchGeometry.boundingBox() );
180
181 if ( !fids.isEmpty() )
182 {
183 std::unique_ptr< QgsGeometryEngine > engine( QgsGeometry::createGeometryEngine( mSearchGeometry.constGet() ) );
184 engine->prepareGeometry();
185
186 double absDensity = 0;
187 for ( const QgsFeatureId id : fids )
188 {
189 const QgsGeometry lineGeom = mIndex.geometry( id );
190
191 if ( engine->intersects( lineGeom.constGet() ) )
192 {
193 double analysisLineLength = 0;
194 try
195 {
196 analysisLineLength = mDa.measureLength( QgsGeometry( engine->intersection( mIndex.geometry( id ).constGet() ) ) );
197 }
198 catch ( QgsCsException & )
199 {
200 throw QgsProcessingException( QObject::tr( "An error occurred while calculating feature length" ) );
201 }
202
203 double weight = 1;
204
205 if ( !mWeightField.isEmpty() )
206 {
207 weight = mFeatureWeights.value( id );
208 }
209
210 absDensity += ( analysisLineLength * weight );
211 }
212 }
213
214 double lineDensity = 0;
215 if ( absDensity > 0 )
216 {
217 //only calculate ellipsoidal area if abs density is greater 0
218 double analysisSearchGeometryArea = 0;
219 try
220 {
221 analysisSearchGeometryArea = mDa.measureArea( mSearchGeometry );
222 }
223 catch ( QgsCsException & )
224 {
225 throw QgsProcessingException( QObject::tr( "An error occurred while calculating feature area" ) );
226 }
227
228 lineDensity = absDensity / analysisSearchGeometryArea;
229 }
230 rasterDataLine->setValue( 0, col, lineDensity );
231 }
232 else
233 {
234 //no lines found in search radius
235 rasterDataLine->setValue( 0, col, 0.0 );
236 }
237
238 feedback->setProgress( static_cast<double>( cellcnt ) / static_cast<double>( totalCellcnt ) * 100 );
239 cellcnt++;
240 }
241 provider->writeBlock( rasterDataLine.get(), 1, 0, row );
242
243 //'carriage return and newline' for search geometry
244 mSearchGeometry.translate( ( cols - 1 ) * -mPixelSize, -mPixelSize );
245 }
246
247 QVariantMap outputs;
248 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
249 return outputs;
250}
251
252
@ VectorLine
Vector line layers.
@ Numeric
Accepts numeric fields.
@ Float32
Thirty two bit floating point (float)
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
Circle geometry type.
Definition qgscircle.h:43
Custom exception class for Coordinate Reference System related exceptions.
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.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
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:58
QgsFeatureId id
Definition qgsfeature.h:66
Q_INVOKABLE QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Container of fields for a vector layer.
Definition qgsfields.h:46
A geometry is the spatial representation of a feature.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry, double precision=0.0)
Creates and returns a new geometry engine representing the specified geometry using precision on a gr...
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.
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.
A raster layer destination parameter, for specifying the destination path for a raster layer created ...
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)
Sets the name of the data provider for the raster output.
void setOutputFormat(const QString &format)
Sets the output format.
void setCreateOptions(const QStringList &list)
Sets a list of data source creation options to use when creating the output raster file.
QgsRasterDataProvider * createOneBandRaster(Qgis::DataType dataType, int width, int height, const QgsRectangle &extent, const QgsCoordinateReferenceSystem &crs) SIP_FACTORY
Create a raster file with one band without initializing the pixel data.
A rectangle specified with double values.
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:6465
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features