QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
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" ), QgsProcessingParameterField::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" ), QgsProcessingParameterField::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() | QgsProcessingParameterDefinition::FlagAdvanced );
76 distanceParam->setDefaultUnit( QgsUnitTypes::DistanceKilometers );
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() | QgsProcessingParameterDefinition::FlagAdvanced );
84 addParameter( breakParam.release() );
85
86 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Hub lines" ), QgsProcessing::TypeVectorLine ) );
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 = parameterAsFields( 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 = parameterAsFields( 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() == QgsWkbTypes::PointGeometry && !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 for ( int j = 0; j < hubFeature.attributes().count(); ++j )
263 {
264 if ( !hubFieldIndices.contains( j ) )
265 continue;
266 hubAttributes << hubFeature.attribute( j );
267 }
268
269 QgsFeatureRequest spokeRequest = QgsFeatureRequest().setDestinationCrs( hubSource->sourceCrs(), context.transformContext() );
270 spokeRequest.setSubsetOfAttributes( spokeFields2Fetch );
271 spokeRequest.setFilterExpression( QgsExpression::createFieldEqualityExpression( fieldSpokeName, hubFeature.attribute( fieldHubIndex ) ) );
272
273 QgsFeatureIterator spokeFeatures = spokeSource->getFeatures( spokeRequest, QgsProcessingFeatureSource::FlagSkipGeometryValidityChecks );
274 QgsFeature spokeFeature;
275 while ( spokeFeatures.nextFeature( spokeFeature ) )
276 {
277 if ( feedback->isCanceled() )
278 {
279 break;
280 }
281 if ( !spokeFeature.hasGeometry() )
282 continue;
283
284 const QgsPoint spokePoint = getPointFromFeature( spokeFeature );
285 QgsGeometry line;
286 if ( !geodesic )
287 {
288 line = QgsGeometry( new QgsLineString( QVector< QgsPoint >() << hubPoint << spokePoint ) );
289 if ( splitAntimeridian )
290 line = da.splitGeometryAtAntimeridian( line );
291 }
292 else
293 {
294 double distance = geodesicDistance;
295 if ( dynamicGeodesicDistance )
296 {
297 expressionContext.setFeature( hubFeature );
298 distance = geodesicDistanceProperty.valueAsDouble( expressionContext, distance );
299 }
300
301 std::unique_ptr< QgsMultiLineString > ml = std::make_unique< QgsMultiLineString >();
302 std::unique_ptr< QgsLineString > l = std::make_unique< QgsLineString >( QVector< QgsPoint >() << hubPoint );
303 const QVector< QVector< QgsPointXY > > points = da.geodesicLine( QgsPointXY( hubPoint ), QgsPointXY( spokePoint ), distance, splitAntimeridian );
304 QVector< QgsPointXY > points1 = points.at( 0 );
305 points1.pop_front();
306 if ( points.count() == 1 )
307 points1.pop_back();
308
309 const QgsLineString geodesicPoints( points1 );
310 l->append( &geodesicPoints );
311 if ( points.count() == 1 )
312 l->addVertex( spokePoint );
313
314 ml->addGeometry( l.release() );
315 if ( points.count() > 1 )
316 {
317 QVector< QgsPointXY > points2 = points.at( 1 );
318 points2.pop_back();
319 l = std::make_unique< QgsLineString >( points2 );
320 if ( hasZ )
321 l->addZValue( std::numeric_limits<double>::quiet_NaN() );
322 if ( hasM )
323 l->addMValue( std::numeric_limits<double>::quiet_NaN() );
324
325 l->addVertex( spokePoint );
326 ml->addGeometry( l.release() );
327 }
328 line = QgsGeometry( std::move( ml ) );
329 }
330
331 QgsFeature outFeature;
332 QgsAttributes outAttributes = hubAttributes;
333
334 // only keep selected attributes
335 QgsAttributes spokeAttributes;
336 for ( int j = 0; j < spokeFeature.attributes().count(); ++j )
337 {
338 if ( !spokeFieldIndices.contains( j ) )
339 continue;
340 spokeAttributes << spokeFeature.attribute( j );
341 }
342
343 outAttributes.append( spokeAttributes );
344 outFeature.setAttributes( outAttributes );
345 outFeature.setGeometry( line );
346 if ( !sink->addFeature( outFeature, QgsFeatureSink::FastInsert ) )
347 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
348 }
349 }
350
351 QVariantMap outputs;
352 outputs.insert( QStringLiteral( "OUTPUT" ), dest );
353 return outputs;
354}
355
bool is3D() const SIP_HOLDGIL
Returns true if the geometry is 3D and contains a z-value.
bool isMeasure() const SIP_HOLDGIL
Returns true if the geometry contains m values.
A vector of attributes.
Definition: qgsattributes.h:59
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, QVariant::Type fieldType=QVariant::Type::Invalid)
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)
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:56
QgsAttributes attributes
Definition: qgsfeature.h:65
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
Definition: qgsfeature.cpp:160
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:233
QVariant attribute(const QString &name) const
Lookup attribute value by attribute name.
Definition: qgsfeature.cpp:338
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
Definition: qgsfeature.cpp:170
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
bool append(const QgsField &field, FieldOrigin origin=OriginProvider, int originIndex=-1)
Appends a field. The field must have unique name, otherwise it is rejected (returns false)
Definition: qgsfields.cpp:59
int count() const
Returns number of items.
Definition: qgsfields.cpp:133
int lookupField(const QString &fieldName) const
Looks up field's index from the field name.
Definition: qgsfields.cpp:349
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:164
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:45
A class to represent a 2D point.
Definition: qgspointxy.h:59
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:562
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:551
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
@ FlagSkipGeometryValidityChecks
Invalid geometry checks should always be skipped. This flag can be useful for algorithms which always...
Base class for providing feedback from a processing algorithm.
A boolean parameter for processing algorithms.
@ FlagAdvanced
Parameter is an advanced parameter which should be hidden from users by default.
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).
@ TypeVectorLine
Vector line layers.
Definition: qgsprocessing.h:50
Definition for a property.
Definition: qgsproperty.h:46
@ DoublePositive
Positive double value (including 0)
Definition: qgsproperty.h:57
A store for object properties.
Definition: qgsproperty.h:230
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.
@ DistanceKilometers
Kilometers.
Definition: qgsunittypes.h:70
static bool hasM(Type type) SIP_HOLDGIL
Tests whether a WKB type contains m values.
Definition: qgswkbtypes.h:1130
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:70
static Type addZ(Type type) SIP_HOLDGIL
Adds the z dimension to a WKB type and returns the new type.
Definition: qgswkbtypes.h:1176
static bool hasZ(Type type) SIP_HOLDGIL
Tests whether a WKB type contains the z-dimension.
Definition: qgswkbtypes.h:1080
static Type addM(Type type) SIP_HOLDGIL
Adds the m dimension to a WKB type and returns the new type.
Definition: qgswkbtypes.h:1201
QList< int > QgsAttributeList
Definition: qgsfield.h:26
const QgsField & field
Definition: qgsfield.h:463