[FEATURE][processing] Add geodesic mode for "Join by Lines (Hub lines)" algorithm

This allows optional creation of geodesic lines, which represent the
shortest distance between the points based on the ellipsoid.

When geodesic mode is used, it is possible to split the created lines
at the antimeridian (±180 degrees longitude), which can improve
rendering of the lines. Additionally, the distance between vertices
can be specified. A smaller distance results in a denser, more accurate
line.

python/plugins/processing/tests/testdata/custom/geodesic_hub.cpg

@@ -0,0 +1 @@
+UTF-8

python/plugins/processing/tests/testdata/custom/geodesic_hub.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

python/plugins/processing/tests/testdata/custom/geodesic_hub.qpj

@@ -0,0 +1 @@
+PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

python/plugins/processing/tests/testdata/custom/geodesic_spoke.cpg

@@ -0,0 +1 @@
+UTF-8

python/plugins/processing/tests/testdata/custom/geodesic_spoke.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

python/plugins/processing/tests/testdata/custom/geodesic_spoke.qpj

@@ -0,0 +1 @@
+PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

python/plugins/processing/tests/testdata/expected/geodesic_hub_no_split.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

python/plugins/processing/tests/testdata/expected/geodesic_hub_no_split.qpj

@@ -0,0 +1 @@
+PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

python/plugins/processing/tests/testdata/expected/geodesic_hub_split.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

python/plugins/processing/tests/testdata/expected/geodesic_hub_split.qpj

@@ -0,0 +1 @@
+PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]

python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml

@@ -2988,6 +2988,48 @@ tests:
         name: expected/hub_lines_field_subset.gml
         type: vector
 
+  - algorithm: native:hublines
+    name: Hub lines geodesic
+    ellipsoid: GRS80
+    project_crs: EPSG:4326
+    params:
+      ANTIMERIDIAN_SPLIT: false
+      GEODESIC: true
+      GEODESIC_DISTANCE: 1000.0
+      HUBS:
+        name: custom/geodesic_hub.shp
+        type: vector
+      HUB_FIELD: id
+      SPOKES:
+        name: custom/geodesic_spoke.shp
+        type: vector
+      SPOKE_FIELD: id
+    results:
+      OUTPUT:
+        name: expected/geodesic_hub_no_split.shp
+        type: vector
+
+  - algorithm: native:hublines
+    name: Hub lines, geodesic, split
+    ellipsoid: GRS80
+    project_crs: EPSG:4326
+    params:
+      ANTIMERIDIAN_SPLIT: true
+      GEODESIC: true
+      GEODESIC_DISTANCE: 1000.0
+      HUBS:
+        name: custom/geodesic_hub.shp
+        type: vector
+      HUB_FIELD: id
+      SPOKES:
+        name: custom/geodesic_spoke.shp
+        type: vector
+      SPOKE_FIELD: id
+    results:
+      OUTPUT:
+        name: expected/geodesic_hub_split.shp
+        type: vector
+
   - algorithm: qgis:pointstopath
     name: Points to path (non grouped)
     params:

src/analysis/processing/qgsalgorithmjoinwithlines.cpp

@@ -17,6 +17,7 @@
 
 #include "qgsalgorithmjoinwithlines.h"
 #include "qgslinestring.h"
+#include "qgsmultilinestring.h"
 
 ///@cond PRIVATE
 
@@ -32,7 +33,7 @@ QString QgsJoinWithLinesAlgorithm::displayName() const
 
 QStringList QgsJoinWithLinesAlgorithm::tags() const
 {
-  return QObject::tr( "join,connect,lines,points,hub,spoke" ).split( ',' );
+  return QObject::tr( "join,connect,lines,points,hub,spoke,geodesic,great,circle" ).split( ',' );
 }
 
 QString QgsJoinWithLinesAlgorithm::group() const
@@ -67,6 +68,20 @@ void QgsJoinWithLinesAlgorithm::initAlgorithm( const QVariantMap & )
                 QVariant(), QStringLiteral( "SPOKES" ), QgsProcessingParameterField::Any,
                 true, true ) );
 
+  addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "GEODESIC" ), QObject::tr( "Create geodesic lines" ), false ) );
+
+  auto distanceParam = qgis::make_unique< QgsProcessingParameterDistance >( QStringLiteral( "GEODESIC_DISTANCE" ), QObject::tr( "Distance between vertices (geodesic lines only)" ), 1000 );
+  distanceParam->setFlags( distanceParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
+  distanceParam->setDefaultUnit( QgsUnitTypes::DistanceKilometers );
+  distanceParam->setIsDynamic( true );
+  distanceParam->setDynamicPropertyDefinition( QgsPropertyDefinition( QStringLiteral( "Geodesic Distance" ), QObject::tr( "Distance between vertices" ), QgsPropertyDefinition::DoublePositive ) );
+  distanceParam->setDynamicLayerParameterName( QStringLiteral( "HUBS" ) );
+  addParameter( distanceParam.release() );
+
+  auto breakParam = qgis::make_unique< QgsProcessingParameterBoolean >( QStringLiteral( "ANTIMERIDIAN_SPLIT" ), QObject::tr( "Split lines at antimeridian (±180 degrees longitude)" ), false );
+  breakParam->setFlags( breakParam->flags() | QgsProcessingParameterDefinition::FlagAdvanced );
+  addParameter( breakParam.release() );
+
   addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Hub lines" ), QgsProcessing::TypeVectorLine ) );
 }
 
@@ -74,7 +89,16 @@ QString QgsJoinWithLinesAlgorithm::shortHelpString() const
 {
   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"
                       "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"
-                      "If input layers are not point layers, a point on the surface of the geometries will be taken as the connecting location." );
+                      "If input layers are not point layers, a point on the surface of the geometries will be taken as the connecting location.\n\n"
+                      "Optionally, geodesic lines can be created, which represent the shortest path on the surface of an ellipsoid. When "
+                      "geodesic mode is used, it is possible to split the created lines at the antimeridian (±180 degrees longitude), which can improve "
+                      "rendering of the lines. Additionally, the distance between vertices can be specified. A smaller distance results in a denser, more "
+                      "accurate line." );
+}
+
+QString QgsJoinWithLinesAlgorithm::shortDescription() const
+{
+  return QObject::tr( "Creates lines joining two point layers, based on a common attribute value." );
 }
 
 QgsJoinWithLinesAlgorithm *QgsJoinWithLinesAlgorithm::createInstance() const
@@ -106,6 +130,21 @@ QVariantMap QgsJoinWithLinesAlgorithm::processAlgorithm( const QVariantMap &para
   if ( fieldHubIndex < 0 || fieldSpokeIndex < 0 )
     throw QgsProcessingException( QObject::tr( "Invalid ID field" ) );
 
+  const bool geodesic = parameterAsBool( parameters, QStringLiteral( "GEODESIC" ), context );
+  const double geodesicDistance = parameterAsDouble( parameters, QStringLiteral( "GEODESIC_DISTANCE" ), context ) * 1000;
+  bool dynamicGeodesicDistance = QgsProcessingParameters::isDynamic( parameters, QStringLiteral( "GEODESIC_DISTANCE" ) );
+  QgsExpressionContext expressionContext = createExpressionContext( parameters, context, hubSource.get() );
+  QgsProperty geodesicDistanceProperty;
+  if ( dynamicGeodesicDistance )
+  {
+    geodesicDistanceProperty = parameters.value( QStringLiteral( "GEODESIC_DISTANCE" ) ).value< QgsProperty >();
+  }
+
+  const bool splitAntimeridian = parameterAsBool( parameters, QStringLiteral( "ANTIMERIDIAN_SPLIT" ), context );
+  QgsDistanceArea da;
+  da.setSourceCrs( hubSource->sourceCrs(), context.transformContext() );
+  da.setEllipsoid( context.project()->ellipsoid() );
+
   QgsFields hubOutFields;
   QgsAttributeList hubFieldIndices;
   if ( hubFieldsToCopy.empty() )
@@ -139,6 +178,7 @@ QVariantMap QgsJoinWithLinesAlgorithm::processAlgorithm( const QVariantMap &para
   if ( spokeFieldsToCopy.empty() )
   {
     spokeOutFields = spokeSource->fields();
+    spokeFieldIndices.reserve( spokeOutFields.count() );
     for ( int i = 0; i < spokeOutFields.count(); ++i )
     {
       spokeFieldIndices << i;
@@ -163,7 +203,7 @@ QVariantMap QgsJoinWithLinesAlgorithm::processAlgorithm( const QVariantMap &para
 
   QgsFields fields = QgsProcessingUtils::combineFields( hubOutFields, spokeOutFields );
 
-  QgsWkbTypes::Type outType = QgsWkbTypes::LineString;
+  QgsWkbTypes::Type outType = geodesic ? QgsWkbTypes::MultiLineString : QgsWkbTypes::LineString;
   bool hasZ = false;
   if ( QgsWkbTypes::hasZ( hubSource->wkbType() ) || QgsWkbTypes::hasZ( spokeSource->wkbType() ) )
   {
@@ -241,7 +281,51 @@ QVariantMap QgsJoinWithLinesAlgorithm::processAlgorithm( const QVariantMap &para
         continue;
 
       QgsPoint spokePoint = getPointFromFeature( spokeFeature );
-      QgsGeometry line( new QgsLineString( QVector< QgsPoint >() << hubPoint << spokePoint ) );
+      QgsGeometry line;
+      if ( !geodesic )
+      {
+        line = QgsGeometry( new QgsLineString( QVector< QgsPoint >() << hubPoint << spokePoint ) );
+        if ( splitAntimeridian )
+          line = da.splitGeometryAtAntimeridian( line );
+      }
+      else
+      {
+        double distance = geodesicDistance;
+        if ( dynamicGeodesicDistance )
+        {
+          expressionContext.setFeature( hubFeature );
+          distance = geodesicDistanceProperty.valueAsDouble( expressionContext, distance );
+        }
+
+        std::unique_ptr< QgsMultiLineString > ml = qgis::make_unique< QgsMultiLineString >();
+        std::unique_ptr< QgsLineString > l = qgis::make_unique< QgsLineString >( QVector< QgsPoint >() << hubPoint );
+        QVector< QVector< QgsPointXY > > points = da.geodesicLine( QgsPointXY( hubPoint ), QgsPointXY( spokePoint ), distance, splitAntimeridian );
+        QVector< QgsPointXY > points1 = points.at( 0 );
+        points1.pop_front();
+        if ( points.count() == 1 )
+          points1.pop_back();
+
+        QgsLineString geodesicPoints( points1 );
+        l->append( &geodesicPoints );
+        if ( points.count() == 1 )
+          l->addVertex( spokePoint );
+
+        ml->addGeometry( l.release() );
+        if ( points.count() > 1 )
+        {
+          QVector< QgsPointXY > points2 = points.at( 1 );
+          points2.pop_back();
+          l = qgis::make_unique< QgsLineString >( points2 );
+          if ( hasZ )
+            l->addZValue( std::numeric_limits<double>::quiet_NaN() );
+          if ( hasM )
+            l->addMValue( std::numeric_limits<double>::quiet_NaN() );
+
+          l->addVertex( spokePoint );
+          ml->addGeometry( l.release() );
+        }
+        line = QgsGeometry( std::move( ml ) );
+      }
 
       QgsFeature outFeature;
       QgsAttributes outAttributes = hubAttributes;

src/analysis/processing/qgsalgorithmjoinwithlines.h

@@ -41,6 +41,7 @@ class QgsJoinWithLinesAlgorithm : public QgsProcessingAlgorithm
     QString group() const override;
     QString groupId() const override;
     QString shortHelpString() const override;
+    QString shortDescription() const override;
     QgsJoinWithLinesAlgorithm *createInstance() const override SIP_FACTORY;
 
   protected: