QGIS API Documentation 3.41.0-Master (88383c3d16f)
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 auto 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 // round up width and height to the nearest integer as GDAL does (e.g. in gdal_rasterize)
140 // see https://github.com/qgis/QGIS/issues/43547
141 const int rows = static_cast<int>( 0.5 + mExtent.height() / mPixelSize );
142 const int cols = static_cast<int>( 0.5 + mExtent.width() / mPixelSize );
143
144 //build new raster extent based on number of columns and cellsize
145 //this prevents output cellsize being calculated too small
146 const QgsRectangle rasterExtent = QgsRectangle( mExtent.xMinimum(), mExtent.yMaximum() - ( rows * mPixelSize ), mExtent.xMinimum() + ( cols * mPixelSize ), mExtent.yMaximum() );
147
148 QgsRasterFileWriter writer = QgsRasterFileWriter( outputFile );
149 writer.setOutputProviderKey( QStringLiteral( "gdal" ) );
150 writer.setOutputFormat( outputFormat );
151 if ( !createOptions.isEmpty() )
152 {
153 writer.setCreateOptions( createOptions.split( '|' ) );
154 }
155
156 std::unique_ptr<QgsRasterDataProvider> provider( writer.createOneBandRaster( Qgis::DataType::Float32, cols, rows, rasterExtent, mCrs ) );
157 if ( !provider )
158 throw QgsProcessingException( QObject::tr( "Could not create raster output: %1" ).arg( outputFile ) );
159 if ( !provider->isValid() )
160 throw QgsProcessingException( QObject::tr( "Could not create raster output %1: %2" ).arg( outputFile, provider->error().message( QgsErrorMessage::Text ) ) );
161
162 provider->setNoDataValue( 1, -9999 );
163
164 const qgssize totalCellcnt = static_cast<qgssize>( rows ) * cols;
165 int cellcnt = 0;
166
167 auto rasterDataLine = std::make_unique<QgsRasterBlock>( Qgis::DataType::Float32, cols, 1 );
168
169 for ( int row = 0; row < rows; row++ )
170 {
171 for ( int col = 0; col < cols; col++ )
172 {
173 if ( feedback->isCanceled() )
174 {
175 break;
176 }
177
178 if ( col > 0 )
179 mSearchGeometry.translate( mPixelSize, 0 );
180
181 const QList<QgsFeatureId> fids = mIndex.intersects( mSearchGeometry.boundingBox() );
182
183 if ( !fids.isEmpty() )
184 {
185 std::unique_ptr<QgsGeometryEngine> engine( QgsGeometry::createGeometryEngine( mSearchGeometry.constGet() ) );
186 engine->prepareGeometry();
187
188 double absDensity = 0;
189 for ( const QgsFeatureId id : fids )
190 {
191 const QgsGeometry lineGeom = mIndex.geometry( id );
192
193 if ( engine->intersects( lineGeom.constGet() ) )
194 {
195 double analysisLineLength = 0;
196 try
197 {
198 analysisLineLength = mDa.measureLength( QgsGeometry( engine->intersection( mIndex.geometry( id ).constGet() ) ) );
199 }
200 catch ( QgsCsException & )
201 {
202 throw QgsProcessingException( QObject::tr( "An error occurred while calculating feature length" ) );
203 }
204
205 double weight = 1;
206
207 if ( !mWeightField.isEmpty() )
208 {
209 weight = mFeatureWeights.value( id );
210 }
211
212 absDensity += ( analysisLineLength * weight );
213 }
214 }
215
216 double lineDensity = 0;
217 if ( absDensity > 0 )
218 {
219 //only calculate ellipsoidal area if abs density is greater 0
220 double analysisSearchGeometryArea = 0;
221 try
222 {
223 analysisSearchGeometryArea = mDa.measureArea( mSearchGeometry );
224 }
225 catch ( QgsCsException & )
226 {
227 throw QgsProcessingException( QObject::tr( "An error occurred while calculating feature area" ) );
228 }
229
230 lineDensity = absDensity / analysisSearchGeometryArea;
231 }
232 rasterDataLine->setValue( 0, col, lineDensity );
233 }
234 else
235 {
236 //no lines found in search radius
237 rasterDataLine->setValue( 0, col, 0.0 );
238 }
239
240 feedback->setProgress( static_cast<double>( cellcnt ) / static_cast<double>( totalCellcnt ) * 100 );
241 cellcnt++;
242 }
243 provider->writeBlock( rasterDataLine.get(), 1, 0, row );
244
245 //'carriage return and newline' for search geometry
246 mSearchGeometry.translate( ( cols - 1 ) * -mPixelSize, -mPixelSize );
247 }
248
249 QVariantMap outputs;
250 outputs.insert( QStringLiteral( "OUTPUT" ), outputFile );
251 return outputs;
252}
253
254
@ 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:45
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, Qgis::GeosCreationFlags flags=Qgis::GeosCreationFlag::SkipEmptyInteriorRings)
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:6673
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features