port Hub distance algorithm to C++

python/plugins/processing/algs/qgis/HubDistanceLines.py

@@ -39,6 +39,7 @@ from qgis.core import (
     QgsProcessingParameterFeatureSink,
     QgsProcessingException,
     QgsSpatialIndex,
+    QgsProcessingAlgorithm,
 )
 from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
 
@@ -74,6 +75,9 @@ class HubDistanceLines(QgisAlgorithm):
     def __init__(self):
         super().__init__()
 
+    def flags(self):
+        return super().flags() | QgsProcessingAlgorithm.Flag.FlagDeprecated
+
     def initAlgorithm(self, config=None):
         self.units = [
             self.tr("Meters"),

python/plugins/processing/algs/qgis/HubDistancePoints.py

@@ -39,6 +39,7 @@ from qgis.core import (
     QgsProcessingParameterEnum,
     QgsProcessingParameterFeatureSink,
     QgsProcessingException,
+    QgsProcessingAlgorithm,
 )
 from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
 
@@ -71,6 +72,9 @@ class HubDistancePoints(QgisAlgorithm):
     def __init__(self):
         super().__init__()
 
+    def flags(self):
+        return super().flags() | QgsProcessingAlgorithm.Flag.FlagDeprecated
+
     def initAlgorithm(self, config=None):
         self.units = [
             self.tr("Meters"),

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

@@ -1496,8 +1496,8 @@ tests:
         name: expected/gridify_lines.gml
         type: vector
 
-  - algorithm: qgis:distancetonearesthubpoints
-    name: Hub distance points
+  - algorithm: native:distancetonearesthub
+    name: Hub distance
     params:
       INPUT:
         name: points.gml
@@ -1508,25 +1508,12 @@ tests:
       FIELD: name
       UNIT: 0
     results:
-      OUTPUT:
-        name: expected/hub_distance_points.gml
-        type: vector
-
-  - algorithm: qgis:distancetonearesthublinetohub
-    name: Hub distance lines
-    params:
-      INPUT:
-        name: points.gml
-        type: vector
-      HUBS:
-        name: custom/hub_points.gml
-        type: vector
-      FIELD: name
-      UNIT: 0
-    results:
-      OUTPUT:
+      OUTPUT_LINES:
         name: expected/hub_distance_lines.gml
         type: vector
+      OUTPUT_POINTS:
+        name: expected/hub_distance_points.gml
+        type: vector
 
   - algorithm: native:hublines
     name: Hub lines

src/analysis/CMakeLists.txt

@@ -164,6 +164,7 @@ set(QGIS_ANALYSIS_SRCS
   processing/qgsalgorithmgpsbabeltools.cpp
   processing/qgsalgorithmgrid.cpp
   processing/qgsalgorithmhillshade.cpp
+  processing/qgsalgorithmhubdistance.cpp
   processing/qgsalgorithmimportphotos.cpp
   processing/qgsalgorithminterpolatepoint.cpp
   processing/qgsalgorithmintersection.cpp

src/analysis/processing/qgsalgorithmhubdistance.cpp

@@ -0,0 +1,244 @@
+/***************************************************************************
+                         qgsalgorithmhubdistance.cpp
+                         ---------------------
+    begin                : April 2025
+    copyright            : (C) 2025 by Alexander Bruy
+    email                : alexander dot bruy at gmail dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "qgsalgorithmhubdistance.h"
+#include "qgsfeaturerequest.h"
+#include "qgsspatialindex.h"
+
+///@cond PRIVATE
+
+QString QgsHubDistanceAlgorithm::name() const
+{
+  return QStringLiteral( "distancetonearesthub" );
+}
+
+QString QgsHubDistanceAlgorithm::displayName() const
+{
+  return QObject::tr( "Distance to nearest hub" );
+}
+
+QStringList QgsHubDistanceAlgorithm::tags() const
+{
+  return QObject::tr( "lines,points,hub,spoke,geodesic,great,circle,distance" ).split( ',' );
+}
+
+QString QgsHubDistanceAlgorithm::group() const
+{
+  return QObject::tr( "Vector analysis" );
+}
+
+QString QgsHubDistanceAlgorithm::groupId() const
+{
+  return QStringLiteral( "vectoranalysis" );
+}
+
+QString QgsHubDistanceAlgorithm::shortHelpString() const
+{
+  return QObject::tr( "Computes the distance between features from the source layer to the closest feature "
+                      "from the destination layer.\n\n"
+                      "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" );
+}
+
+QgsHubDistanceAlgorithm *QgsHubDistanceAlgorithm::createInstance() const
+{
+  return new QgsHubDistanceAlgorithm();
+}
+
+void QgsHubDistanceAlgorithm::initAlgorithm( const QVariantMap & )
+{
+  addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Source layer (spokes)" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorAnyGeometry ) ) );
+  addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "HUBS" ), QObject::tr( "Destination layer (hubs)" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorAnyGeometry ) ) );
+  addParameter( new QgsProcessingParameterField( QStringLiteral( "FIELD" ), QObject::tr( "Hub layer name attribute" ), QVariant(), QStringLiteral( "HUBS" ) ) );
+
+  const QStringList options = QStringList()
+                              << QObject::tr( "Meters" )
+                              << QObject::tr( "Feets" )
+                              << QObject::tr( "Miles" )
+                              << QObject::tr( "Kilometers" )
+                              << QObject::tr( "Layer units" );
+  addParameter( new QgsProcessingParameterEnum( QStringLiteral( "UNIT" ), QObject::tr( "Measurement unit" ), options, false, 0 ) );
+  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_LINES" ), QObject::tr( "Hub lines" ), Qgis::ProcessingSourceType::VectorLine, QVariant(), true, true ) );
+  addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT_POINTS" ), QObject::tr( "Hub points" ), Qgis::ProcessingSourceType::VectorPoint, QVariant(), true, false ) );
+}
+
+QVariantMap QgsHubDistanceAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
+{
+  if ( parameters.value( QStringLiteral( "INPUT" ) ) == parameters.value( QStringLiteral( "HUBS" ) ) )
+    throw QgsProcessingException( QObject::tr( "Same layer given for both hubs and spokes" ) );
+
+  std::unique_ptr<QgsProcessingFeatureSource> hubSource( parameterAsSource( parameters, QStringLiteral( "HUBS" ), context ) );
+  if ( !hubSource )
+    throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "HUBS" ) ) );
+
+  std::unique_ptr<QgsProcessingFeatureSource> spokeSource( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
+  if ( !spokeSource )
+    throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
+
+  const QString fieldHubName = parameterAsString( parameters, QStringLiteral( "FIELD" ), context );
+  const int hubNameIndex = hubSource->fields().lookupField( fieldHubName );
+
+  const int unitIndex = parameterAsEnum( parameters, QStringLiteral( "UNIT" ), context );
+  Qgis::DistanceUnit unit = Qgis::DistanceUnit::Unknown;
+  switch ( unitIndex )
+  {
+    case 0:
+      unit = Qgis::DistanceUnit::Meters;
+      break;
+    case 1:
+      unit = Qgis::DistanceUnit::Feet;
+      break;
+    case 2:
+      unit = Qgis::DistanceUnit::Miles;
+      break;
+    case 3:
+      unit = Qgis::DistanceUnit::Kilometers;
+      break;
+  }
+
+  QgsFields fields = spokeSource->fields();
+  fields.append( QgsField( QStringLiteral( "hub_name" ), QMetaType::Type::QString ) );
+  fields.append( QgsField( QStringLiteral( "hub_distance" ), QMetaType::Type::Double ) );
+
+  QString linesDest;
+  std::unique_ptr<QgsFeatureSink> linesSink( parameterAsSink( parameters, QStringLiteral( "OUTPUT_LINES" ), context, linesDest, fields, Qgis::WkbType::LineString, hubSource->sourceCrs(), QgsFeatureSink::RegeneratePrimaryKey ) );
+  if ( !linesSink )
+    throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_LINES" ) ) );
+
+  QString pointsDest;
+  std::unique_ptr<QgsFeatureSink> pointsSink( parameterAsSink( parameters, QStringLiteral( "OUTPUT_POINTS" ), context, pointsDest, fields, Qgis::WkbType::Point, hubSource->sourceCrs(), QgsFeatureSink::RegeneratePrimaryKey ) );
+  if ( !pointsSink )
+    throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT_POINTS" ) ) );
+
+  QgsFeatureRequest request;
+  request.setSubsetOfAttributes( QgsAttributeList() << hubNameIndex );
+  request.setDestinationCrs( spokeSource->sourceCrs(), context.transformContext() );
+
+  const QgsSpatialIndex hubsIndex( hubSource->getFeatures( request ), feedback, QgsSpatialIndex::FlagStoreFeatureGeometries );
+  QHash<QgsFeatureId, QVariant> hubsAttributeCache;
+  double step = hubSource->featureCount() > 0 ? 50.0 / hubSource->featureCount() : 1;
+  long long i = 0;
+  const QgsSpatialIndex idx( hubSource->getFeatures( request ), [&]( const QgsFeature &f ) -> bool {
+    if ( feedback-> isCanceled() )
+    {
+      return false;
+    }
+
+    hubsAttributeCache.insert( f.id(), f.attributes().at( hubNameIndex ) );
+
+    i++;
+    feedback->setProgress( i * step );
+    return true; }, QgsSpatialIndex::FlagStoreFeatureGeometries );
+
+  QgsDistanceArea da;
+  da.setSourceCrs( spokeSource->sourceCrs(), context.transformContext() );
+  da.setEllipsoid( context.ellipsoid() );
+
+  // Scan source points, find nearest hub, and write to output file
+  i = 0;
+  QgsFeature spokeFeature, outputFeature;
+  step = spokeSource->featureCount() > 0 ? 50.0 / spokeSource->featureCount() : 1;
+  QgsFeatureIterator features = spokeSource->getFeatures();
+  while ( features.nextFeature( spokeFeature ) )
+  {
+    if ( feedback->isCanceled() )
+    {
+      break;
+    }
+
+    i++;
+    feedback->setProgress( i * step );
+
+    if ( !spokeFeature.hasGeometry() )
+    {
+      if ( linesSink && !linesSink->addFeature( spokeFeature, QgsFeatureSink::Flag::FastInsert ) )
+      {
+        throw QgsProcessingException( writeFeatureError( linesSink.get(), parameters, QStringLiteral( "OUTPUT_LINES" ) ) );
+      }
+      if ( pointsSink && !pointsSink->addFeature( spokeFeature, QgsFeatureSink::Flag::FastInsert ) )
+      {
+        throw QgsProcessingException( writeFeatureError( pointsSink.get(), parameters, QStringLiteral( "OUTPUT_POINTS" ) ) );
+      }
+      continue;
+    }
+
+    QgsPointXY point = spokeFeature.geometry().boundingBox().center();
+
+    const QList<QgsFeatureId> neighbors = hubsIndex.nearestNeighbor( point, 1 );
+    if ( neighbors.isEmpty() )
+    {
+      continue;
+    }
+
+    const QgsPointXY hub = hubsIndex.geometry( neighbors.at( 0 ) ).boundingBox().center();
+    double hubDistance = da.measureLine( point, hub );
+
+    if ( unit != Qgis::DistanceUnit::Unknown )
+    {
+      hubDistance = da.convertLengthMeasurement( hubDistance, unit );
+    }
+
+    QgsAttributes attrs = spokeFeature.attributes();
+    attrs << hubsAttributeCache.value( neighbors.at( 0 ) ) << hubDistance;
+    outputFeature = QgsFeature();
+    outputFeature.setAttributes( attrs );
+
+    if ( linesSink )
+    {
+      outputFeature.setGeometry( QgsGeometry::fromPolylineXY( QgsPolylineXY() << point << hub ) );
+      if ( !linesSink->addFeature( outputFeature, QgsFeatureSink::Flag::FastInsert ) )
+      {
+        throw QgsProcessingException( writeFeatureError( linesSink.get(), parameters, QStringLiteral( "OUTPUT_LINES" ) ) );
+      }
+    }
+
+    if ( pointsSink )
+    {
+      outputFeature.setGeometry( QgsGeometry::fromPointXY( point ) );
+      if ( !pointsSink->addFeature( outputFeature, QgsFeatureSink::Flag::FastInsert ) )
+      {
+        throw QgsProcessingException( writeFeatureError( pointsSink.get(), parameters, QStringLiteral( "OUTPUT_POINTS" ) ) );
+      }
+    }
+  }
+
+  if ( linesSink )
+  {
+    linesSink->finalize();
+  }
+  if ( pointsSink )
+  {
+    pointsSink->finalize();
+  }
+
+  QVariantMap results;
+  if ( linesSink )
+  {
+    results.insert( QStringLiteral( "OUTPUT_LINES" ), linesDest );
+  }
+  if ( pointsSink )
+  {
+    results.insert( QStringLiteral( "OUTPUT_POINTS" ), pointsDest );
+  }
+  return results;
+}
+
+///@endcond

src/analysis/processing/qgsalgorithmhubdistance.h

@@ -0,0 +1,51 @@
+/***************************************************************************
+                         qgsalgorithmhubdistance.h
+                         ---------------------
+    begin                : April 2025
+    copyright            : (C) 2025 by Alexander Bruy
+    email                : alexander dot bruy at gmail dot com
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+
+#ifndef QGSALGORITHMHUBDISTANCE_H
+#define QGSALGORITHMHUBDISTANCE_H
+
+#define SIP_NO_FILE
+
+#include "qgis_sip.h"
+#include "qgsprocessingalgorithm.h"
+
+///@cond PRIVATE
+
+/**
+ * Native hub distance algorithm.
+ */
+class QgsHubDistanceAlgorithm : public QgsProcessingAlgorithm
+{
+  public:
+    QgsHubDistanceAlgorithm() = default;
+    void initAlgorithm( const QVariantMap &configuration = QVariantMap() ) override;
+    QString name() const override;
+    QString displayName() const override;
+    QStringList tags() const override;
+    QString group() const override;
+    QString groupId() const override;
+    QString shortHelpString() const override;
+    QgsHubDistanceAlgorithm *createInstance() const override SIP_FACTORY;
+
+  protected:
+    QVariantMap processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback ) override;
+};
+
+///@endcond PRIVATE
+
+#endif // QGSALGORITHMHUBDISTANCE_H

src/analysis/processing/qgsnativealgorithms.cpp

@@ -148,6 +148,7 @@
 #endif
 #include "qgsalgorithmgrid.h"
 #include "qgsalgorithmhillshade.h"
+#include "qgsalgorithmhubdistance.h"
 #include "qgsalgorithmjoinbyattribute.h"
 #include "qgsalgorithmjoinbylocation.h"
 #include "qgsalgorithmjoinbylocationsummary.h"
@@ -481,6 +482,7 @@ void QgsNativeAlgorithms::loadAlgorithms()
 #endif
   addAlgorithm( new QgsGridAlgorithm() );
   addAlgorithm( new QgsHillshadeAlgorithm() );
+  addAlgorithm( new QgsHubDistanceAlgorithm() );
   addAlgorithm( new QgsImportPhotosAlgorithm() );
   addAlgorithm( new QgsInterpolatePointAlgorithm() );
   addAlgorithm( new QgsIntersectionAlgorithm() );