Skip to content
Permalink
Browse files

[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.
  • Loading branch information
nyalldawson committed Jan 16, 2019
1 parent 1248bea commit f89d061ba0c36474c6a89e246cf1eb4bb80b8133
Showing with 141 additions and 4 deletions.
  1. +1 −0 python/plugins/processing/tests/testdata/custom/geodesic_hub.cpg
  2. BIN python/plugins/processing/tests/testdata/custom/geodesic_hub.dbf
  3. +1 −0 python/plugins/processing/tests/testdata/custom/geodesic_hub.prj
  4. +1 −0 python/plugins/processing/tests/testdata/custom/geodesic_hub.qpj
  5. BIN python/plugins/processing/tests/testdata/custom/geodesic_hub.shp
  6. BIN python/plugins/processing/tests/testdata/custom/geodesic_hub.shx
  7. +1 −0 python/plugins/processing/tests/testdata/custom/geodesic_spoke.cpg
  8. BIN python/plugins/processing/tests/testdata/custom/geodesic_spoke.dbf
  9. +1 −0 python/plugins/processing/tests/testdata/custom/geodesic_spoke.prj
  10. +1 −0 python/plugins/processing/tests/testdata/custom/geodesic_spoke.qpj
  11. BIN python/plugins/processing/tests/testdata/custom/geodesic_spoke.shp
  12. BIN python/plugins/processing/tests/testdata/custom/geodesic_spoke.shx
  13. BIN python/plugins/processing/tests/testdata/expected/geodesic_hub_no_split.dbf
  14. +1 −0 python/plugins/processing/tests/testdata/expected/geodesic_hub_no_split.prj
  15. +1 −0 python/plugins/processing/tests/testdata/expected/geodesic_hub_no_split.qpj
  16. BIN python/plugins/processing/tests/testdata/expected/geodesic_hub_no_split.shp
  17. BIN python/plugins/processing/tests/testdata/expected/geodesic_hub_no_split.shx
  18. BIN python/plugins/processing/tests/testdata/expected/geodesic_hub_split.dbf
  19. +1 −0 python/plugins/processing/tests/testdata/expected/geodesic_hub_split.prj
  20. +1 −0 python/plugins/processing/tests/testdata/expected/geodesic_hub_split.qpj
  21. BIN python/plugins/processing/tests/testdata/expected/geodesic_hub_split.shp
  22. BIN python/plugins/processing/tests/testdata/expected/geodesic_hub_split.shx
  23. +42 −0 python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml
  24. +88 −4 src/analysis/processing/qgsalgorithmjoinwithlines.cpp
  25. +1 −0 src/analysis/processing/qgsalgorithmjoinwithlines.h
@@ -0,0 +1 @@
UTF-8
Binary file not shown.
@@ -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]]
@@ -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"]]
Binary file not shown.
Binary file not shown.
@@ -0,0 +1 @@
UTF-8
Binary file not shown.
@@ -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]]
@@ -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"]]
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -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]]
@@ -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"]]
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -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]]
@@ -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"]]
Binary file not shown.
Binary file not shown.
@@ -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:
@@ -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,14 +68,37 @@ 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 ) );
}

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;
@@ -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:

0 comments on commit f89d061

Please sign in to comment.
You can’t perform that action at this time.