QGIS API Documentation 3.38.0-Grenoble (exported)
Loading...
Searching...
No Matches
qgsalgorithmjoinwithlines.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmjoinwithlines.cpp
3 ---------------------
4 begin : April 2017
5 copyright : (C) 2017 by Nyall Dawson
6 email : nyall dot dawson 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 "qgslinestring.h"
20#include "qgsmultilinestring.h"
21#include "qgsdistancearea.h"
22
24
25QString QgsJoinWithLinesAlgorithm::name() const
26{
27 return QStringLiteral( "hublines" );
28}
29
30QString QgsJoinWithLinesAlgorithm::displayName() const
31{
32 return QObject::tr( "Join by lines (hub lines)" );
33}
34
35QStringList QgsJoinWithLinesAlgorithm::tags() const
36{
37 return QObject::tr( "join,connect,lines,points,hub,spoke,geodesic,great,circle" ).split( ',' );
38}
39
40QString QgsJoinWithLinesAlgorithm::group() const
41{
42 return QObject::tr( "Vector analysis" );
43}
44
45QString QgsJoinWithLinesAlgorithm::groupId() const
46{
47 return QStringLiteral( "vectoranalysis" );
48}
49
50void QgsJoinWithLinesAlgorithm::initAlgorithm( const QVariantMap & )
51{
52 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "HUBS" ),
53 QObject::tr( "Hub layer" ) ) );
54 addParameter( new QgsProcessingParameterField( QStringLiteral( "HUB_FIELD" ),
55 QObject::tr( "Hub ID field" ), QVariant(), QStringLiteral( "HUBS" ) ) );
56
57 addParameter( new QgsProcessingParameterField( QStringLiteral( "HUB_FIELDS" ),
58 QObject::tr( "Hub layer fields to copy (leave empty to copy all fields)" ),
59 QVariant(), QStringLiteral( "HUBS" ), Qgis::ProcessingFieldParameterDataType::Any,
60 true, true ) );
61
62 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "SPOKES" ),
63 QObject::tr( "Spoke layer" ) ) );
64 addParameter( new QgsProcessingParameterField( QStringLiteral( "SPOKE_FIELD" ),
65 QObject::tr( "Spoke ID field" ), QVariant(), QStringLiteral( "SPOKES" ) ) );
66
67 addParameter( new QgsProcessingParameterField( QStringLiteral( "SPOKE_FIELDS" ),
68 QObject::tr( "Spoke layer fields to copy (leave empty to copy all fields)" ),
69 QVariant(), QStringLiteral( "SPOKES" ), Qgis::ProcessingFieldParameterDataType::Any,
70 true, true ) );
71
72 addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "GEODESIC" ), QObject::tr( "Create geodesic lines" ), false ) );
73
74 auto distanceParam = std::make_unique< QgsProcessingParameterDistance >( QStringLiteral( "GEODESIC_DISTANCE" ), QObject::tr( "Distance between vertices (geodesic lines only)" ), 1000 );
75 distanceParam->setFlags( distanceParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
76 distanceParam->setDefaultUnit( Qgis::DistanceUnit::Kilometers );
77 distanceParam->setIsDynamic( true );
78 distanceParam->setDynamicPropertyDefinition( QgsPropertyDefinition( QStringLiteral( "Geodesic Distance" ), QObject::tr( "Distance between vertices" ), QgsPropertyDefinition::DoublePositive ) );
79 distanceParam->setDynamicLayerParameterName( QStringLiteral( "HUBS" ) );
80 addParameter( distanceParam.release() );
81
82 auto breakParam = std::make_unique< QgsProcessingParameterBoolean >( QStringLiteral( "ANTIMERIDIAN_SPLIT" ), QObject::tr( "Split lines at antimeridian (±180 degrees longitude)" ), false );
83 breakParam->setFlags( breakParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
84 addParameter( breakParam.release() );
85
86 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Hub lines" ), Qgis::ProcessingSourceType::VectorLine ) );
87}
88
89QString QgsJoinWithLinesAlgorithm::shortHelpString() const
90{
91 return QObject::tr( "This algorithm creates hub and spoke diagrams by connecting lines from points on the Spoke layer to matching points in the Hub layer.\n\n"
92 "Determination of which hub goes with each point is based on a match between the Hub ID field on the hub points and the Spoke ID field on the spoke points.\n\n"
93 "If input layers are not point layers, a point on the surface of the geometries will be taken as the connecting location.\n\n"
94 "Optionally, geodesic lines can be created, which represent the shortest path on the surface of an ellipsoid. When "
95 "geodesic mode is used, it is possible to split the created lines at the antimeridian (±180 degrees longitude), which can improve "
96 "rendering of the lines. Additionally, the distance between vertices can be specified. A smaller distance results in a denser, more "
97 "accurate line." );
98}
99
100QString QgsJoinWithLinesAlgorithm::shortDescription() const
101{
102 return QObject::tr( "Creates lines joining two point layers, based on a common attribute value." );
103}
104
105QgsJoinWithLinesAlgorithm *QgsJoinWithLinesAlgorithm::createInstance() const
106{
107 return new QgsJoinWithLinesAlgorithm();
108}
109
110QVariantMap QgsJoinWithLinesAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
111{
112 if ( parameters.value( QStringLiteral( "SPOKES" ) ) == parameters.value( QStringLiteral( "HUBS" ) ) )
113 throw QgsProcessingException( QObject::tr( "Same layer given for both hubs and spokes" ) );
114
115 std::unique_ptr< QgsProcessingFeatureSource > hubSource( parameterAsSource( parameters, QStringLiteral( "HUBS" ), context ) );
116 if ( !hubSource )
117 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "HUBS" ) ) );
118
119 std::unique_ptr< QgsProcessingFeatureSource > spokeSource( parameterAsSource( parameters, QStringLiteral( "SPOKES" ), context ) );
120 if ( !spokeSource )
121 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "SPOKES" ) ) );
122
123 const QString fieldHubName = parameterAsString( parameters, QStringLiteral( "HUB_FIELD" ), context );
124 const int fieldHubIndex = hubSource->fields().lookupField( fieldHubName );
125 const QStringList hubFieldsToCopy = parameterAsStrings( parameters, QStringLiteral( "HUB_FIELDS" ), context );
126
127 const QString fieldSpokeName = parameterAsString( parameters, QStringLiteral( "SPOKE_FIELD" ), context );
128 const int fieldSpokeIndex = spokeSource->fields().lookupField( fieldSpokeName );
129 const QStringList spokeFieldsToCopy = parameterAsStrings( parameters, QStringLiteral( "SPOKE_FIELDS" ), context );
130
131 if ( fieldHubIndex < 0 || fieldSpokeIndex < 0 )
132 throw QgsProcessingException( QObject::tr( "Invalid ID field" ) );
133
134 const bool geodesic = parameterAsBoolean( parameters, QStringLiteral( "GEODESIC" ), context );
135 const double geodesicDistance = parameterAsDouble( parameters, QStringLiteral( "GEODESIC_DISTANCE" ), context ) * 1000;
136 const bool dynamicGeodesicDistance = QgsProcessingParameters::isDynamic( parameters, QStringLiteral( "GEODESIC_DISTANCE" ) );
137 QgsExpressionContext expressionContext = createExpressionContext( parameters, context, hubSource.get() );
138 QgsProperty geodesicDistanceProperty;
139 if ( dynamicGeodesicDistance )
140 {
141 geodesicDistanceProperty = parameters.value( QStringLiteral( "GEODESIC_DISTANCE" ) ).value< QgsProperty >();
142 }
143
144 const bool splitAntimeridian = parameterAsBoolean( parameters, QStringLiteral( "ANTIMERIDIAN_SPLIT" ), context );
146 da.setSourceCrs( hubSource->sourceCrs(), context.transformContext() );
147 da.setEllipsoid( context.ellipsoid() );
148
149 QgsFields hubOutFields;
150 QgsAttributeList hubFieldIndices;
151 if ( hubFieldsToCopy.empty() )
152 {
153 hubOutFields = hubSource->fields();
154 hubFieldIndices.reserve( hubOutFields.count() );
155 for ( int i = 0; i < hubOutFields.count(); ++i )
156 {
157 hubFieldIndices << i;
158 }
159 }
160 else
161 {
162 hubFieldIndices.reserve( hubOutFields.count() );
163 for ( const QString &field : hubFieldsToCopy )
164 {
165 const int index = hubSource->fields().lookupField( field );
166 if ( index >= 0 )
167 {
168 hubFieldIndices << index;
169 hubOutFields.append( hubSource->fields().at( index ) );
170 }
171 }
172 }
173
174 QgsAttributeList hubFields2Fetch = hubFieldIndices;
175 hubFields2Fetch << fieldHubIndex;
176
177 QgsFields spokeOutFields;
178 QgsAttributeList spokeFieldIndices;
179 if ( spokeFieldsToCopy.empty() )
180 {
181 spokeOutFields = spokeSource->fields();
182 spokeFieldIndices.reserve( spokeOutFields.count() );
183 for ( int i = 0; i < spokeOutFields.count(); ++i )
184 {
185 spokeFieldIndices << i;
186 }
187 }
188 else
189 {
190 for ( const QString &field : spokeFieldsToCopy )
191 {
192 const int index = spokeSource->fields().lookupField( field );
193 if ( index >= 0 )
194 {
195 spokeFieldIndices << index;
196 spokeOutFields.append( spokeSource->fields().at( index ) );
197 }
198 }
199 }
200
201 QgsAttributeList spokeFields2Fetch = spokeFieldIndices;
202 spokeFields2Fetch << fieldSpokeIndex;
203
204
205 const QgsFields fields = QgsProcessingUtils::combineFields( hubOutFields, spokeOutFields );
206
208 bool hasZ = false;
209 if ( !geodesic && ( QgsWkbTypes::hasZ( hubSource->wkbType() ) || QgsWkbTypes::hasZ( spokeSource->wkbType() ) ) )
210 {
211 outType = QgsWkbTypes::addZ( outType );
212 hasZ = true;
213 }
214 bool hasM = false;
215 if ( !geodesic && ( QgsWkbTypes::hasM( hubSource->wkbType() ) || QgsWkbTypes::hasM( spokeSource->wkbType() ) ) )
216 {
217 outType = QgsWkbTypes::addM( outType );
218 hasM = true;
219 }
220
221 QString dest;
222 std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields,
223 outType, hubSource->sourceCrs(), QgsFeatureSink::RegeneratePrimaryKey ) );
224 if ( !sink )
225 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
226
227 auto getPointFromFeature = [hasZ, hasM]( const QgsFeature & feature )->QgsPoint
228 {
229 QgsPoint p;
230 if ( feature.geometry().type() == Qgis::GeometryType::Point && !feature.geometry().isMultipart() )
231 p = *static_cast< const QgsPoint *>( feature.geometry().constGet() );
232 else
233 p = *static_cast< const QgsPoint *>( feature.geometry().pointOnSurface().constGet() );
234 if ( hasZ && !p.is3D() )
235 p.addZValue( 0 );
236 if ( hasM && !p.isMeasure() )
237 p.addMValue( 0 );
238 return p;
239 };
240
242 const double step = hubSource->featureCount() > 0 ? 100.0 / hubSource->featureCount() : 1;
243 int i = 0;
244 QgsFeature hubFeature;
245 while ( hubFeatures.nextFeature( hubFeature ) )
246 {
247 i++;
248 if ( feedback->isCanceled() )
249 {
250 break;
251 }
252
253 feedback->setProgress( i * step );
254
255 if ( !hubFeature.hasGeometry() )
256 continue;
257
258 const QgsPoint hubPoint = getPointFromFeature( hubFeature );
259
260 // only keep selected attributes
261 QgsAttributes hubAttributes;
262 const int attributeCount = hubFeature.attributeCount();
263 for ( int j = 0; j < attributeCount; ++j )
264 {
265 if ( !hubFieldIndices.contains( j ) )
266 continue;
267 hubAttributes << hubFeature.attribute( j );
268 }
269
270 QgsFeatureRequest spokeRequest = QgsFeatureRequest().setDestinationCrs( hubSource->sourceCrs(), context.transformContext() );
271 spokeRequest.setSubsetOfAttributes( spokeFields2Fetch );
272 spokeRequest.setFilterExpression( QgsExpression::createFieldEqualityExpression( fieldSpokeName, hubFeature.attribute( fieldHubIndex ) ) );
273
274 QgsFeatureIterator spokeFeatures = spokeSource->getFeatures( spokeRequest, Qgis::ProcessingFeatureSourceFlag::SkipGeometryValidityChecks );
275 QgsFeature spokeFeature;
276 while ( spokeFeatures.nextFeature( spokeFeature ) )
277 {
278 if ( feedback->isCanceled() )
279 {
280 break;
281 }
282 if ( !spokeFeature.hasGeometry() )
283 continue;
284
285 const QgsPoint spokePoint = getPointFromFeature( spokeFeature );
286 QgsGeometry line;
287 if ( !geodesic )
288 {
289 line = QgsGeometry( new QgsLineString( QVector< QgsPoint >() << hubPoint << spokePoint ) );
290 if ( splitAntimeridian )
291 line = da.splitGeometryAtAntimeridian( line );
292 }
293 else
294 {
295 double distance = geodesicDistance;
296 if ( dynamicGeodesicDistance )
297 {
298 expressionContext.setFeature( hubFeature );
299 distance = geodesicDistanceProperty.valueAsDouble( expressionContext, distance );
300 }
301
302 std::unique_ptr< QgsMultiLineString > ml = std::make_unique< QgsMultiLineString >();
303 std::unique_ptr< QgsLineString > l = std::make_unique< QgsLineString >( QVector< QgsPoint >() << hubPoint );
304 const QVector< QVector< QgsPointXY > > points = da.geodesicLine( QgsPointXY( hubPoint ), QgsPointXY( spokePoint ), distance, splitAntimeridian );
305 QVector< QgsPointXY > points1 = points.at( 0 );
306 points1.pop_front();
307 if ( points.count() == 1 )
308 points1.pop_back();
309
310 const QgsLineString geodesicPoints( points1 );
311 l->append( &geodesicPoints );
312 if ( points.count() == 1 )
313 l->addVertex( spokePoint );
314
315 ml->addGeometry( l.release() );
316 if ( points.count() > 1 )
317 {
318 QVector< QgsPointXY > points2 = points.at( 1 );
319 points2.pop_back();
320 l = std::make_unique< QgsLineString >( points2 );
321 if ( hasZ )
322 l->addZValue( std::numeric_limits<double>::quiet_NaN() );
323 if ( hasM )
324 l->addMValue( std::numeric_limits<double>::quiet_NaN() );
325
326 l->addVertex( spokePoint );
327 ml->addGeometry( l.release() );
328 }
329 line = QgsGeometry( std::move( ml ) );
330 }
331
332 QgsFeature outFeature;
333 QgsAttributes outAttributes = hubAttributes;
334
335 // only keep selected attributes
336 QgsAttributes spokeAttributes;
337 const int attributeCount = spokeFeature.attributeCount();
338 for ( int j = 0; j < attributeCount; ++j )
339 {
340 if ( !spokeFieldIndices.contains( j ) )
341 continue;
342 spokeAttributes << spokeFeature.attribute( j );
343 }
344
345 outAttributes.append( spokeAttributes );
346 outFeature.setAttributes( outAttributes );
347 outFeature.setGeometry( line );
348 if ( !sink->addFeature( outFeature, QgsFeatureSink::FastInsert ) )
349 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
350 }
351 }
352
353 QVariantMap outputs;
354 outputs.insert( QStringLiteral( "OUTPUT" ), dest );
355 return outputs;
356}
357
@ VectorLine
Vector line layers.
@ Kilometers
Kilometers.
@ SkipGeometryValidityChecks
Invalid geometry checks should always be skipped. This flag can be useful for algorithms which always...
WkbType
The WKB type describes the number of dimensions a geometry has.
Definition qgis.h:201
@ LineString
LineString.
@ MultiLineString
MultiLineString.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
bool isMeasure() const
Returns true if the geometry contains m values.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
A vector of attributes.
A general purpose distance and area calculator, capable of performing ellipsoid based calculations.
QVector< QVector< QgsPointXY > > geodesicLine(const QgsPointXY &p1, const QgsPointXY &p2, double interval, bool breakLine=false) const
Calculates the geodesic line between p1 and p2, which represents the shortest path on the ellipsoid b...
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
QgsGeometry splitGeometryAtAntimeridian(const QgsGeometry &geometry) const
Splits a (Multi)LineString geometry at the antimeridian (longitude +/- 180 degrees).
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
Expression contexts are used to encapsulate the parameters around which a QgsExpression should be eva...
void setFeature(const QgsFeature &feature)
Convenience function for setting a feature for the context.
static QString createFieldEqualityExpression(const QString &fieldName, const QVariant &value, QMetaType::Type fieldType=QMetaType::Type::UnknownType)
Create an expression allowing to evaluate if a field is equal to a value.
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.
QgsFeatureRequest & setDestinationCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets the destination crs for feature's geometries.
QgsFeatureRequest & setFilterExpression(const QString &expression)
Set the filter expression.
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
@ RegeneratePrimaryKey
This flag indicates, that a primary key field cannot be guaranteed to be unique and the sink should i...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
int attributeCount() const
Returns the number of attributes attached to the feature.
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Q_INVOKABLE QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
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
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:60
int count
Definition qgsfields.h:50
Q_INVOKABLE int lookupField(const QString &fieldName) const
Looks up field's index from the field name.
A geometry is the spatial representation of a feature.
Line string geometry type, with support for z-dimension and m-values.
A class to represent a 2D point.
Definition qgspointxy.h:60
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
bool addMValue(double mValue=0) override
Adds a measure to the geometry, initialized to a preset value.
Definition qgspoint.cpp:568
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition qgspoint.cpp:557
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 boolean parameter for processing algorithms.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A vector layer or feature source field parameter for processing algorithms.
static bool isDynamic(const QVariantMap &parameters, const QString &name)
Returns true if the parameter with matching name is a dynamic parameter, and must be evaluated once f...
static QgsFields combineFields(const QgsFields &fieldsA, const QgsFields &fieldsB, const QString &fieldsBPrefix=QString())
Combines two field lists, avoiding duplicate field names (in a case-insensitive manner).
Definition for a property.
Definition qgsproperty.h:45
@ DoublePositive
Positive double value (including 0)
Definition qgsproperty.h:56
A store for object properties.
QVariant value(const QgsExpressionContext &context, const QVariant &defaultValue=QVariant(), bool *ok=nullptr) const
Calculates the current value of the property, including any transforms which are set for the property...
double valueAsDouble(const QgsExpressionContext &context, double defaultValue=0.0, bool *ok=nullptr) const
Calculates the current value of the property and interprets it as a double.
static Qgis::WkbType addM(Qgis::WkbType type)
Adds the m dimension to a WKB type and returns the new type.
static Qgis::WkbType addZ(Qgis::WkbType type)
Adds the z dimension to a WKB type and returns the new type.
static bool hasZ(Qgis::WkbType type)
Tests whether a WKB type contains the z-dimension.
static bool hasM(Qgis::WkbType type)
Tests whether a WKB type contains m values.
QList< int > QgsAttributeList
Definition qgsfield.h:27