[FEATURE][API] Add method to QgsDistanceArea to split a (multi)line g…

…eometry at the antimeridian Whenever line segments in the input geometry cross the antimeridian, they will be split into two segments, with the latitude of the breakpoint being determined using a geodesic line connecting the points either side of this segment. If the geometry contains M or Z values, these will be linearly interpolated for the new vertices created at the antimeridian.
qgis · Jan 11, 2019 · a8bd12c · a8bd12c
1 parent 4d609ff
commit a8bd12c
Show file tree

Hide file tree

Showing 4 changed files with 254 additions and 0 deletions.
diff --git a/python/core/auto_generated/qgsdistancearea.sip.in b/python/core/auto_generated/qgsdistancearea.sip.in
@@ -407,6 +407,35 @@ will also be in this same CRS.
 :return: - the latitude at which the geodesic crosses the date line
          - fractionAlongLine:  will be set to the fraction along the geodesic line joining ``p1`` to ``p2`` at which the date line crossing occurs.
 
+.. seealso:: :py:func:`breakGeometryAtDateLine`
+
+
+.. versionadded:: 3.6
+%End
+
+    QgsGeometry splitGeometryAtDateLine( const QgsGeometry &geometry ) const;
+%Docstring
+Splits a (Multi)LineString ``geometry`` at the international date line (longitude +/- 180 degrees).
+The returned geometry will always be a multi-part geometry.
+
+Whenever line segments in the input geometry cross the international date line, they will be
+split into two segments, with the latitude of the breakpoint being determined using a geodesic
+line connecting the points either side of this segment.
+
+The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
+
+``geometry`` must be in the sourceCrs() of this QgsDistanceArea object. The returned geometry
+will also be in this same CRS.
+
+If ``geometry`` contains M or Z values, these will be linearly interpolated for the new vertices
+created at the date line.
+
+.. note::
+
+   Non-(Multi)LineString geometries will be returned unchanged.
+
+.. seealso:: :py:func:`latitudeGeodesicCrossesDateLine`
+
 .. versionadded:: 3.6
 %End
 

diff --git a/src/core/qgsdistancearea.cpp b/src/core/qgsdistancearea.cpp
@@ -33,6 +33,7 @@
 #include "qgssurface.h"
 #include "qgsunittypes.h"
 #include "qgsexception.h"
+#include "qgsmultilinestring.h"
 
 #include <geodesic.h>
 
@@ -532,6 +533,133 @@ double QgsDistanceArea::latitudeGeodesicCrossesDateLine( const QgsPointXY &pp1,
   return lat;
 }
 
+QgsGeometry QgsDistanceArea::splitGeometryAtDateLine( const QgsGeometry &geometry ) const
+{
+  if ( QgsWkbTypes::geometryType( geometry.wkbType() ) != QgsWkbTypes::LineGeometry )
+    return geometry;
+
+  QgsGeometry g = geometry;
+  // TODO - avoid segmentization of curved geometries (if this is even possible!)
+  if ( QgsWkbTypes::isCurvedType( g.wkbType() ) )
+    g.convertToStraightSegment();
+
+  std::unique_ptr< QgsMultiLineString > res = qgis::make_unique< QgsMultiLineString >();
+  for ( auto part = g.const_parts_begin(); part != g.const_parts_end(); ++part )
+  {
+    const QgsLineString *line = qgsgeometry_cast< const QgsLineString * >( *part );
+    if ( !line )
+      continue;
+    if ( line->isEmpty() )
+    {
+      continue;
+    }
+
+    std::unique_ptr< QgsLineString > l = qgis::make_unique< QgsLineString >();
+    try
+    {
+      double x = 0;
+      double y = 0;
+      double z = 0;
+      double m = 0;
+      QVector< QgsPoint > newPoints;
+      newPoints.reserve( line->numPoints() );
+      double prevLon = 0;
+      double prevLat = 0;
+      double lon = 0;
+      double lat = 0;
+      double prevZ = 0;
+      double prevM = 0;
+      for ( int i = 0; i < line->numPoints(); i++ )
+      {
+        QgsPoint p = line->pointN( i );
+        x = p.x();
+        if ( mCoordTransform.sourceCrs().isGeographic() )
+        {
+          x = std::fmod( x, 360.0 );
+          if ( x > 180 )
+            x -= 360;
+          p.setX( x );
+        }
+        y = p.y();
+        lon = x;
+        lat = y;
+        mCoordTransform.transformInPlace( lon, lat, z );
+
+        //test if we crossed the dateline in this segment
+        if ( i > 0 && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon  < -120 ) ) )
+        {
+          // we did!
+          // when crossing the antimeridian, we need to calculate the latitude
+          // at which the geodesic intersects the antimeridian
+          double fract = 0;
+          double lat180 = latitudeGeodesicCrossesDateLine( QgsPointXY( prevLon, prevLat ), QgsPointXY( lon, lat ), fract );
+          if ( line->is3D() )
+          {
+            z = prevZ + ( p.z() - prevZ ) * fract;
+          }
+          if ( line->isMeasure() )
+          {
+            m = prevM + ( p.m() - prevM ) * fract;
+          }
+
+          QgsPointXY antiMeridianPoint;
+          if ( prevLon < -120 )
+            antiMeridianPoint = mCoordTransform.transform( QgsPointXY( -180, lat180 ), QgsCoordinateTransform::ReverseTransform );
+          else
+            antiMeridianPoint = mCoordTransform.transform( QgsPointXY( 180, lat180 ), QgsCoordinateTransform::ReverseTransform );
+
+          QgsPoint newPoint( antiMeridianPoint );
+          if ( line->is3D() )
+            newPoint.addZValue( z );
+          if ( line->isMeasure() )
+            newPoint.addMValue( m );
+
+          if ( std::isfinite( newPoint.x() ) && std::isfinite( newPoint.y() ) )
+          {
+            newPoints << newPoint;
+          }
+          res->addGeometry( new QgsLineString( newPoints ) );
+
+          newPoints.clear();
+          newPoints.reserve( line->numPoints() - i + 1 );
+
+          if ( lon < -120 )
+            antiMeridianPoint = mCoordTransform.transform( QgsPointXY( -180, lat180 ), QgsCoordinateTransform::ReverseTransform );
+          else
+            antiMeridianPoint = mCoordTransform.transform( QgsPointXY( 180, lat180 ), QgsCoordinateTransform::ReverseTransform );
+
+          if ( std::isfinite( antiMeridianPoint.x() ) && std::isfinite( antiMeridianPoint.y() ) )
+          {
+            // we want to keep the previously calculated z/m value for newPoint, if present. They're the same each
+            // of the antimeridian split
+            newPoint.setX( antiMeridianPoint.x() );
+            newPoint.setY( antiMeridianPoint.y() );
+            newPoints << newPoint;
+          }
+        }
+        newPoints << p;
+
+        prevLon = lon;
+        prevLat = lat;
+        if ( line->is3D() )
+          prevZ = p.z();
+        if ( line->isMeasure() )
+          prevM = p.m();
+      }
+      res->addGeometry( new QgsLineString( newPoints ) );
+    }
+    catch ( QgsCsException & )
+    {
+      QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform linestring. Unable to calculate break point." ) );
+      res->addGeometry( line->clone() );
+      break;
+    }
+  }
+
+  return QgsGeometry( std::move( res ) );
+}
+
+
 QVector< QVector<QgsPointXY> > QgsDistanceArea::geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, const double interval, const bool breakLine ) const
 {
   if ( !willUseEllipsoid() )

diff --git a/src/core/qgsdistancearea.h b/src/core/qgsdistancearea.h
@@ -333,10 +333,35 @@ class CORE_EXPORT QgsDistanceArea
      * \param fractionAlongLine will be set to the fraction along the geodesic line joining \a p1 to \a p2 at which the date line crossing occurs.
      *
      * \returns the latitude at which the geodesic crosses the date line
+     * \see breakGeometryAtDateLine()
+     *
      * \since QGIS 3.6
      */
     double latitudeGeodesicCrossesDateLine( const QgsPointXY &p1, const QgsPointXY &p2, double &fractionAlongLine SIP_OUT ) const;
 
+    /**
+     * Splits a (Multi)LineString \a geometry at the international date line (longitude +/- 180 degrees).
+     * The returned geometry will always be a multi-part geometry.
+     *
+     * Whenever line segments in the input geometry cross the international date line, they will be
+     * split into two segments, with the latitude of the breakpoint being determined using a geodesic
+     * line connecting the points either side of this segment.
+     *
+     * The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
+     *
+     * \a geometry must be in the sourceCrs() of this QgsDistanceArea object. The returned geometry
+     * will also be in this same CRS.
+     *
+     * If \a geometry contains M or Z values, these will be linearly interpolated for the new vertices
+     * created at the date line.
+     *
+     * \note Non-(Multi)LineString geometries will be returned unchanged.
+     *
+     * \see latitudeGeodesicCrossesDateLine()
+     * \since QGIS 3.6
+     */
+    QgsGeometry splitGeometryAtDateLine( const QgsGeometry &geometry ) const;
+
   private:
 
     /**

diff --git a/tests/src/python/test_qgsdistancearea.py b/tests/src/python/test_qgsdistancearea.py
@@ -694,6 +694,11 @@ def testGeodesicIntersectionAtDateLine(self):
                                                         QgsPointXY(179.92498611649580198, 7.24703528617311754))
         self.assertAlmostEqual(lat, 7.6112109902580265, 5)
         self.assertAlmostEqual(fract, 0.04111771567489498, 5)
+        lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(360 - 178.20070563806575592, 16.09649962419504732),
+                                                        QgsPointXY(179.92498611649580198, 7.24703528617311754))
+        self.assertAlmostEqual(lat, 7.6112109902580265, 5)
+        self.assertAlmostEqual(fract, 0.04111771567489498, 5)
+
         lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(175.76717768974583578, 8.93749416467257873),
                                                         QgsPointXY(-175.15030911497356669, 8.59851183021221033))
         self.assertAlmostEqual(lat, 8.80683758146703966, 5)
@@ -763,6 +768,73 @@ def testGeodesicLine(self):
         self.assertEqual(g.asWkt(0),
                          'MultiLineString ((-13536427 14138932, -16514348 11691516, -17948849 9406595, -18744235 7552985, -19255354 6014890, -19622372 4688888, -19909239 3505045, -20037508 2933522),(20037508 2933522, 19925702 2415579, 19712755 1385803, 19513769 388441, 19318507 -600065, 19117459 -1602293, 18899973 -2642347, 18651869 -3748726, 18351356 -4958346, 17960498 -6322823, 17404561 -7918366, 16514601 -9855937, 14851845 -12232940, 13760912 -13248201))')
 
+    def testSplitGeometryAtDateline(self):
+        da = QgsDistanceArea()
+        crs = QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId)
+        da.setSourceCrs(crs, QgsProject.instance().transformContext())
+        da.setEllipsoid("WGS84")
+
+        # noops
+        g = da.splitGeometryAtDateLine(QgsGeometry())
+        self.assertTrue(g.isNull())
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('Point(1 2)'))
+        self.assertEqual(g.asWkt(), 'Point (1 2)')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('MultiPoint(1 2, 3 4)'))
+        self.assertEqual(g.asWkt(), 'MultiPoint ((1 2),(3 4))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('PointZ(1 2 3)'))
+        self.assertEqual(g.asWkt(), 'PointZ (1 2 3)')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('PointM(1 2 3)'))
+        self.assertEqual(g.asWkt(), 'PointM (1 2 3)')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString()'))
+        self.assertEqual(g.asWkt(), 'MultiLineString ()')
+
+        # lines
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(0 0, -170 0)'))
+        self.assertEqual(g.asWkt(), 'MultiLineString ((0 0, -170 0))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(-170 0, 0 0)'))
+        self.assertEqual(g.asWkt(), 'MultiLineString ((-170 0, 0 0))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(179 0, -179 0)'))
+        self.assertEqual(g.asWkt(), 'MultiLineString ((179 0, 180 0),(-180 0, -179 0))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(179 0, 181 0)'))
+        self.assertEqual(g.asWkt(), 'MultiLineString ((179 0, 180 0),(-180 0, -179 0))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(-179 0, 179 0)'))
+        self.assertEqual(g.asWkt(), 'MultiLineString ((-179 0, -180 0),(180 0, 179 0))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(181 0, 179 0)'))
+        self.assertEqual(g.asWkt(), 'MultiLineString ((-179 0, -180 0),(180 0, 179 0))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(179 10, -179 -20)'))
+        self.assertEqual(g.asWkt(3), 'MultiLineString ((179 10, 180 -5.362),(-180 -5.362, -179 -20))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(179 -80, -179 70)'))
+        self.assertEqual(g.asWkt(3), 'MultiLineString ((179 -80, 180 -55.685),(-180 -55.685, -179 70))')
+
+        # multiline input
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('MultiLineString((1 10, 50 30),(179 -80, -179 70))'))
+        self.assertEqual(g.asWkt(3), 'MultiLineString ((1 10, 50 30),(179 -80, 180 -55.685),(-180 -55.685, -179 70))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('MultiLineString((1 10, 50 30),(179 -80, 179.99 70))'))
+        self.assertEqual(g.asWkt(3), 'MultiLineString ((1 10, 50 30),(179 -80, 179.99 70))')
+
+        # with z/m
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineStringZ(179 -80 1, -179 70 10)'))
+        self.assertEqual(g.asWkt(3), 'MultiLineStringZ ((179 -80 1, 180 -55.685 2.466),(-180 -55.685 2.466, -179 70 10))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineStringM(179 -80 1, -179 70 10)'))
+        self.assertEqual(g.asWkt(3), 'MultiLineStringM ((179 -80 1, 180 -55.685 2.466),(-180 -55.685 2.466, -179 70 10))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineStringZM(179 -80 1 -4, -179 70 10 -30)'))
+        self.assertEqual(g.asWkt(3), 'MultiLineStringZM ((179 -80 1 -4, 180 -55.685 2.466 -8.234),(-180 -55.685 2.466 -8.234, -179 70 10 -30))')
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('MultiLineStringZ((179 -80 1, -179 70 10),(-170 -5 1, -181 10 5))'))
+        self.assertEqual(g.asWkt(3), 'MultiLineStringZ ((179 -80 1, 180 -55.685 2.466),(-180 -55.685 2.466, -179 70 10),(-170 -5 1, -181 10 5))')
+
+        # different ellipsoid - should change intersection latitude
+        da.setEllipsoid("Phobos2000")
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString(179 10, -179 -20)'))
+        self.assertEqual(g.asWkt(3), 'MultiLineString ((179 10, 180 -5.459),(-180 -5.459, -179 -20))')
+
+        # with reprojection
+        da.setEllipsoid("WGS84")
+        # with reprojection
+        da.setSourceCrs(QgsCoordinateReferenceSystem('EPSG:3857'), QgsProject.instance().transformContext())
+
+        g = da.splitGeometryAtDateLine(QgsGeometry.fromWkt('LineString( -13536427 14138932, 13760912 -13248201)'))
+        self.assertEqual(g.asWkt(1), 'MultiLineString ((-13536427 14138932, -20037508.3 2933521.7),(20037508.3 2933521.7, 13760912 -13248201))')
+
 
 if __name__ == '__main__':
     unittest.main()