Branch data Line data Source code
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 : :
18 : : #include "qgsalgorithmjoinwithlines.h"
19 : : #include "qgslinestring.h"
20 : : #include "qgsmultilinestring.h"
21 : :
22 : : ///@cond PRIVATE
23 : :
24 : 0 : QString QgsJoinWithLinesAlgorithm::name() const
25 : : {
26 : 0 : return QStringLiteral( "hublines" );
27 : : }
28 : :
29 : 0 : QString QgsJoinWithLinesAlgorithm::displayName() const
30 : : {
31 : 0 : return QObject::tr( "Join by lines (hub lines)" );
32 : : }
33 : :
34 : 0 : QStringList QgsJoinWithLinesAlgorithm::tags() const
35 : : {
36 : 0 : return QObject::tr( "join,connect,lines,points,hub,spoke,geodesic,great,circle" ).split( ',' );
37 : 0 : }
38 : :
39 : 0 : QString QgsJoinWithLinesAlgorithm::group() const
40 : : {
41 : 0 : return QObject::tr( "Vector analysis" );
42 : : }
43 : :
44 : 0 : QString QgsJoinWithLinesAlgorithm::groupId() const
45 : : {
46 : 0 : return QStringLiteral( "vectoranalysis" );
47 : : }
48 : :
49 : 0 : void QgsJoinWithLinesAlgorithm::initAlgorithm( const QVariantMap & )
50 : : {
51 : 0 : addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "HUBS" ),
52 : 0 : QObject::tr( "Hub layer" ) ) );
53 : 0 : addParameter( new QgsProcessingParameterField( QStringLiteral( "HUB_FIELD" ),
54 : 0 : QObject::tr( "Hub ID field" ), QVariant(), QStringLiteral( "HUBS" ) ) );
55 : :
56 : 0 : addParameter( new QgsProcessingParameterField( QStringLiteral( "HUB_FIELDS" ),
57 : 0 : QObject::tr( "Hub layer fields to copy (leave empty to copy all fields)" ),
58 : 0 : QVariant(), QStringLiteral( "HUBS" ), QgsProcessingParameterField::Any,
59 : : true, true ) );
60 : :
61 : 0 : addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "SPOKES" ),
62 : 0 : QObject::tr( "Spoke layer" ) ) );
63 : 0 : addParameter( new QgsProcessingParameterField( QStringLiteral( "SPOKE_FIELD" ),
64 : 0 : QObject::tr( "Spoke ID field" ), QVariant(), QStringLiteral( "SPOKES" ) ) );
65 : :
66 : 0 : addParameter( new QgsProcessingParameterField( QStringLiteral( "SPOKE_FIELDS" ),
67 : 0 : QObject::tr( "Spoke layer fields to copy (leave empty to copy all fields)" ),
68 : 0 : QVariant(), QStringLiteral( "SPOKES" ), QgsProcessingParameterField::Any,
69 : : true, true ) );
70 : :
71 : 0 : addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "GEODESIC" ), QObject::tr( "Create geodesic lines" ), false ) );
72 : :
73 : 0 : auto distanceParam = std::make_unique< QgsProcessingParameterDistance >( QStringLiteral( "GEODESIC_DISTANCE" ), QObject::tr( "Distance between vertices (geodesic lines only)" ), 1000 );
74 : 0 : distanceParam->setFlags( distanceParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
75 : 0 : distanceParam->setDefaultUnit( QgsUnitTypes::DistanceKilometers );
76 : 0 : distanceParam->setIsDynamic( true );
77 : 0 : distanceParam->setDynamicPropertyDefinition( QgsPropertyDefinition( QStringLiteral( "Geodesic Distance" ), QObject::tr( "Distance between vertices" ), QgsPropertyDefinition::DoublePositive ) );
78 : 0 : distanceParam->setDynamicLayerParameterName( QStringLiteral( "HUBS" ) );
79 : 0 : addParameter( distanceParam.release() );
80 : :
81 : 0 : auto breakParam = std::make_unique< QgsProcessingParameterBoolean >( QStringLiteral( "ANTIMERIDIAN_SPLIT" ), QObject::tr( "Split lines at antimeridian (±180 degrees longitude)" ), false );
82 : 0 : breakParam->setFlags( breakParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
83 : 0 : addParameter( breakParam.release() );
84 : :
85 : 0 : addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Hub lines" ), QgsProcessing::TypeVectorLine ) );
86 : 0 : }
87 : :
88 : 0 : QString QgsJoinWithLinesAlgorithm::shortHelpString() const
89 : : {
90 : 0 : 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"
91 : : "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"
92 : : "If input layers are not point layers, a point on the surface of the geometries will be taken as the connecting location.\n\n"
93 : : "Optionally, geodesic lines can be created, which represent the shortest path on the surface of an ellipsoid. When "
94 : : "geodesic mode is used, it is possible to split the created lines at the antimeridian (±180 degrees longitude), which can improve "
95 : : "rendering of the lines. Additionally, the distance between vertices can be specified. A smaller distance results in a denser, more "
96 : : "accurate line." );
97 : : }
98 : :
99 : 0 : QString QgsJoinWithLinesAlgorithm::shortDescription() const
100 : : {
101 : 0 : return QObject::tr( "Creates lines joining two point layers, based on a common attribute value." );
102 : : }
103 : :
104 : 0 : QgsJoinWithLinesAlgorithm *QgsJoinWithLinesAlgorithm::createInstance() const
105 : : {
106 : 0 : return new QgsJoinWithLinesAlgorithm();
107 : : }
108 : :
109 : 0 : QVariantMap QgsJoinWithLinesAlgorithm::processAlgorithm( const QVariantMap ¶meters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
110 : : {
111 : 0 : if ( parameters.value( QStringLiteral( "SPOKES" ) ) == parameters.value( QStringLiteral( "HUBS" ) ) )
112 : 0 : throw QgsProcessingException( QObject::tr( "Same layer given for both hubs and spokes" ) );
113 : :
114 : 0 : std::unique_ptr< QgsProcessingFeatureSource > hubSource( parameterAsSource( parameters, QStringLiteral( "HUBS" ), context ) );
115 : 0 : if ( !hubSource )
116 : 0 : throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "HUBS" ) ) );
117 : :
118 : 0 : std::unique_ptr< QgsProcessingFeatureSource > spokeSource( parameterAsSource( parameters, QStringLiteral( "SPOKES" ), context ) );
119 : 0 : if ( !spokeSource )
120 : 0 : throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "SPOKES" ) ) );
121 : :
122 : 0 : QString fieldHubName = parameterAsString( parameters, QStringLiteral( "HUB_FIELD" ), context );
123 : 0 : int fieldHubIndex = hubSource->fields().lookupField( fieldHubName );
124 : 0 : const QStringList hubFieldsToCopy = parameterAsFields( parameters, QStringLiteral( "HUB_FIELDS" ), context );
125 : :
126 : 0 : QString fieldSpokeName = parameterAsString( parameters, QStringLiteral( "SPOKE_FIELD" ), context );
127 : 0 : int fieldSpokeIndex = spokeSource->fields().lookupField( fieldSpokeName );
128 : 0 : const QStringList spokeFieldsToCopy = parameterAsFields( parameters, QStringLiteral( "SPOKE_FIELDS" ), context );
129 : :
130 : 0 : if ( fieldHubIndex < 0 || fieldSpokeIndex < 0 )
131 : 0 : throw QgsProcessingException( QObject::tr( "Invalid ID field" ) );
132 : :
133 : 0 : const bool geodesic = parameterAsBoolean( parameters, QStringLiteral( "GEODESIC" ), context );
134 : 0 : const double geodesicDistance = parameterAsDouble( parameters, QStringLiteral( "GEODESIC_DISTANCE" ), context ) * 1000;
135 : 0 : bool dynamicGeodesicDistance = QgsProcessingParameters::isDynamic( parameters, QStringLiteral( "GEODESIC_DISTANCE" ) );
136 : 0 : QgsExpressionContext expressionContext = createExpressionContext( parameters, context, hubSource.get() );
137 : 0 : QgsProperty geodesicDistanceProperty;
138 : 0 : if ( dynamicGeodesicDistance )
139 : : {
140 : 0 : geodesicDistanceProperty = parameters.value( QStringLiteral( "GEODESIC_DISTANCE" ) ).value< QgsProperty >();
141 : 0 : }
142 : :
143 : 0 : const bool splitAntimeridian = parameterAsBoolean( parameters, QStringLiteral( "ANTIMERIDIAN_SPLIT" ), context );
144 : 0 : QgsDistanceArea da;
145 : 0 : da.setSourceCrs( hubSource->sourceCrs(), context.transformContext() );
146 : 0 : da.setEllipsoid( context.ellipsoid() );
147 : :
148 : 0 : QgsFields hubOutFields;
149 : 0 : QgsAttributeList hubFieldIndices;
150 : 0 : if ( hubFieldsToCopy.empty() )
151 : : {
152 : 0 : hubOutFields = hubSource->fields();
153 : 0 : hubFieldIndices.reserve( hubOutFields.count() );
154 : 0 : for ( int i = 0; i < hubOutFields.count(); ++i )
155 : : {
156 : 0 : hubFieldIndices << i;
157 : 0 : }
158 : 0 : }
159 : : else
160 : : {
161 : 0 : hubFieldIndices.reserve( hubOutFields.count() );
162 : 0 : for ( const QString &field : hubFieldsToCopy )
163 : : {
164 : 0 : int index = hubSource->fields().lookupField( field );
165 : 0 : if ( index >= 0 )
166 : : {
167 : 0 : hubFieldIndices << index;
168 : 0 : hubOutFields.append( hubSource->fields().at( index ) );
169 : 0 : }
170 : : }
171 : : }
172 : :
173 : 0 : QgsAttributeList hubFields2Fetch = hubFieldIndices;
174 : 0 : hubFields2Fetch << fieldHubIndex;
175 : :
176 : 0 : QgsFields spokeOutFields;
177 : 0 : QgsAttributeList spokeFieldIndices;
178 : 0 : if ( spokeFieldsToCopy.empty() )
179 : : {
180 : 0 : spokeOutFields = spokeSource->fields();
181 : 0 : spokeFieldIndices.reserve( spokeOutFields.count() );
182 : 0 : for ( int i = 0; i < spokeOutFields.count(); ++i )
183 : : {
184 : 0 : spokeFieldIndices << i;
185 : 0 : }
186 : 0 : }
187 : : else
188 : : {
189 : 0 : for ( const QString &field : spokeFieldsToCopy )
190 : : {
191 : 0 : int index = spokeSource->fields().lookupField( field );
192 : 0 : if ( index >= 0 )
193 : : {
194 : 0 : spokeFieldIndices << index;
195 : 0 : spokeOutFields.append( spokeSource->fields().at( index ) );
196 : 0 : }
197 : : }
198 : : }
199 : :
200 : 0 : QgsAttributeList spokeFields2Fetch = spokeFieldIndices;
201 : 0 : spokeFields2Fetch << fieldSpokeIndex;
202 : :
203 : :
204 : 0 : QgsFields fields = QgsProcessingUtils::combineFields( hubOutFields, spokeOutFields );
205 : :
206 : 0 : QgsWkbTypes::Type outType = geodesic ? QgsWkbTypes::MultiLineString : QgsWkbTypes::LineString;
207 : 0 : bool hasZ = false;
208 : 0 : if ( QgsWkbTypes::hasZ( hubSource->wkbType() ) || QgsWkbTypes::hasZ( spokeSource->wkbType() ) )
209 : : {
210 : 0 : outType = QgsWkbTypes::addZ( outType );
211 : 0 : hasZ = true;
212 : 0 : }
213 : 0 : bool hasM = false;
214 : 0 : if ( QgsWkbTypes::hasM( hubSource->wkbType() ) || QgsWkbTypes::hasM( spokeSource->wkbType() ) )
215 : : {
216 : 0 : outType = QgsWkbTypes::addM( outType );
217 : 0 : hasM = true;
218 : 0 : }
219 : :
220 : 0 : QString dest;
221 : 0 : std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields,
222 : 0 : outType, hubSource->sourceCrs(), QgsFeatureSink::RegeneratePrimaryKey ) );
223 : 0 : if ( !sink )
224 : 0 : throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
225 : :
226 : 0 : auto getPointFromFeature = [hasZ, hasM]( const QgsFeature & feature )->QgsPoint
227 : : {
228 : 0 : QgsPoint p;
229 : 0 : if ( feature.geometry().type() == QgsWkbTypes::PointGeometry && !feature.geometry().isMultipart() )
230 : 0 : p = *static_cast< const QgsPoint *>( feature.geometry().constGet() );
231 : : else
232 : 0 : p = *static_cast< const QgsPoint *>( feature.geometry().pointOnSurface().constGet() );
233 : 0 : if ( hasZ && !p.is3D() )
234 : 0 : p.addZValue( 0 );
235 : 0 : if ( hasM && !p.isMeasure() )
236 : 0 : p.addMValue( 0 );
237 : 0 : return p;
238 : 0 : };
239 : :
240 : 0 : QgsFeatureIterator hubFeatures = hubSource->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( hubFields2Fetch ), QgsProcessingFeatureSource::FlagSkipGeometryValidityChecks );
241 : 0 : double step = hubSource->featureCount() > 0 ? 100.0 / hubSource->featureCount() : 1;
242 : 0 : int i = 0;
243 : 0 : QgsFeature hubFeature;
244 : 0 : while ( hubFeatures.nextFeature( hubFeature ) )
245 : : {
246 : 0 : i++;
247 : 0 : if ( feedback->isCanceled() )
248 : : {
249 : 0 : break;
250 : : }
251 : :
252 : 0 : feedback->setProgress( i * step );
253 : :
254 : 0 : if ( !hubFeature.hasGeometry() )
255 : 0 : continue;
256 : :
257 : 0 : QgsPoint hubPoint = getPointFromFeature( hubFeature );
258 : :
259 : : // only keep selected attributes
260 : 0 : QgsAttributes hubAttributes;
261 : 0 : for ( int j = 0; j < hubFeature.attributes().count(); ++j )
262 : : {
263 : 0 : if ( !hubFieldIndices.contains( j ) )
264 : 0 : continue;
265 : 0 : hubAttributes << hubFeature.attribute( j );
266 : 0 : }
267 : :
268 : 0 : QgsFeatureRequest spokeRequest = QgsFeatureRequest().setDestinationCrs( hubSource->sourceCrs(), context.transformContext() );
269 : 0 : spokeRequest.setSubsetOfAttributes( spokeFields2Fetch );
270 : 0 : spokeRequest.setFilterExpression( QgsExpression::createFieldEqualityExpression( fieldSpokeName, hubFeature.attribute( fieldHubIndex ) ) );
271 : :
272 : 0 : QgsFeatureIterator spokeFeatures = spokeSource->getFeatures( spokeRequest, QgsProcessingFeatureSource::FlagSkipGeometryValidityChecks );
273 : 0 : QgsFeature spokeFeature;
274 : 0 : while ( spokeFeatures.nextFeature( spokeFeature ) )
275 : : {
276 : 0 : if ( feedback->isCanceled() )
277 : : {
278 : 0 : break;
279 : : }
280 : 0 : if ( !spokeFeature.hasGeometry() )
281 : 0 : continue;
282 : :
283 : 0 : QgsPoint spokePoint = getPointFromFeature( spokeFeature );
284 : 0 : QgsGeometry line;
285 : 0 : if ( !geodesic )
286 : : {
287 : 0 : line = QgsGeometry( new QgsLineString( QVector< QgsPoint >() << hubPoint << spokePoint ) );
288 : 0 : if ( splitAntimeridian )
289 : 0 : line = da.splitGeometryAtAntimeridian( line );
290 : 0 : }
291 : : else
292 : : {
293 : 0 : double distance = geodesicDistance;
294 : 0 : if ( dynamicGeodesicDistance )
295 : : {
296 : 0 : expressionContext.setFeature( hubFeature );
297 : 0 : distance = geodesicDistanceProperty.valueAsDouble( expressionContext, distance );
298 : 0 : }
299 : :
300 : 0 : std::unique_ptr< QgsMultiLineString > ml = std::make_unique< QgsMultiLineString >();
301 : 0 : std::unique_ptr< QgsLineString > l = std::make_unique< QgsLineString >( QVector< QgsPoint >() << hubPoint );
302 : 0 : QVector< QVector< QgsPointXY > > points = da.geodesicLine( QgsPointXY( hubPoint ), QgsPointXY( spokePoint ), distance, splitAntimeridian );
303 : 0 : QVector< QgsPointXY > points1 = points.at( 0 );
304 : 0 : points1.pop_front();
305 : 0 : if ( points.count() == 1 )
306 : 0 : points1.pop_back();
307 : :
308 : 0 : QgsLineString geodesicPoints( points1 );
309 : 0 : l->append( &geodesicPoints );
310 : 0 : if ( points.count() == 1 )
311 : 0 : l->addVertex( spokePoint );
312 : :
313 : 0 : ml->addGeometry( l.release() );
314 : 0 : if ( points.count() > 1 )
315 : : {
316 : 0 : QVector< QgsPointXY > points2 = points.at( 1 );
317 : 0 : points2.pop_back();
318 : 0 : l = std::make_unique< QgsLineString >( points2 );
319 : 0 : if ( hasZ )
320 : 0 : l->addZValue( std::numeric_limits<double>::quiet_NaN() );
321 : 0 : if ( hasM )
322 : 0 : l->addMValue( std::numeric_limits<double>::quiet_NaN() );
323 : :
324 : 0 : l->addVertex( spokePoint );
325 : 0 : ml->addGeometry( l.release() );
326 : 0 : }
327 : 0 : line = QgsGeometry( std::move( ml ) );
328 : 0 : }
329 : :
330 : 0 : QgsFeature outFeature;
331 : 0 : QgsAttributes outAttributes = hubAttributes;
332 : :
333 : : // only keep selected attributes
334 : 0 : QgsAttributes spokeAttributes;
335 : 0 : for ( int j = 0; j < spokeFeature.attributes().count(); ++j )
336 : : {
337 : 0 : if ( !spokeFieldIndices.contains( j ) )
338 : 0 : continue;
339 : 0 : spokeAttributes << spokeFeature.attribute( j );
340 : 0 : }
341 : :
342 : 0 : outAttributes.append( spokeAttributes );
343 : 0 : outFeature.setAttributes( outAttributes );
344 : 0 : outFeature.setGeometry( line );
345 : 0 : sink->addFeature( outFeature, QgsFeatureSink::FastInsert );
346 : 0 : }
347 : 0 : }
348 : :
349 : 0 : QVariantMap outputs;
350 : 0 : outputs.insert( QStringLiteral( "OUTPUT" ), dest );
351 : 0 : return outputs;
352 : 0 : }
353 : :
354 : : ///@endcond
|
