@@ -1854,12 +1854,71 @@ QgsGeometry QgsGeos::shortestLine( const QgsGeometry& other, QString* errorMsg )
1854
1854
return QgsGeometry ( line );
1855
1855
}
1856
1856
1857
+
1858
+ /* * Extract coordinates of linestring's endpoints. Returns false on error. */
1859
+ static bool _linestringEndpoints ( const GEOSGeometry* linestring, double & x1, double & y1, double & x2, double & y2 )
1860
+ {
1861
+ const GEOSCoordSequence* coordSeq = GEOSGeom_getCoordSeq_r ( geosinit.ctxt , linestring );
1862
+ if ( !coordSeq )
1863
+ return false ;
1864
+
1865
+ unsigned int coordSeqSize;
1866
+ if ( GEOSCoordSeq_getSize_r ( geosinit.ctxt , coordSeq, &coordSeqSize ) == 0 )
1867
+ return false ;
1868
+
1869
+ if ( coordSeqSize < 2 )
1870
+ return false ;
1871
+
1872
+ GEOSCoordSeq_getX_r ( geosinit.ctxt , coordSeq, 0 , &x1 );
1873
+ GEOSCoordSeq_getY_r ( geosinit.ctxt , coordSeq, 0 , &y1 );
1874
+ GEOSCoordSeq_getX_r ( geosinit.ctxt , coordSeq, coordSeqSize - 1 , &x2 );
1875
+ GEOSCoordSeq_getY_r ( geosinit.ctxt , coordSeq, coordSeqSize - 1 , &y2 );
1876
+ return true ;
1877
+ }
1878
+
1879
+
1880
+ /* * Merge two linestrings if they meet at the given intersection point, return new geometry or null on error. */
1881
+ static GEOSGeometry* _mergeLinestrings ( const GEOSGeometry* line1, const GEOSGeometry* line2, const QgsPoint& interesectionPoint )
1882
+ {
1883
+ double x1, y1, x2, y2;
1884
+ if ( !_linestringEndpoints ( line1, x1, y1, x2, y2 ) )
1885
+ return nullptr ;
1886
+
1887
+ double rx1, ry1, rx2, ry2;
1888
+ if ( !_linestringEndpoints ( line2, rx1, ry1, rx2, ry2 ) )
1889
+ return nullptr ;
1890
+
1891
+ bool interesectionAtOrigLineEndpoint =
1892
+ ( interesectionPoint.x () == x1 && interesectionPoint.y () == y1 ) ||
1893
+ ( interesectionPoint.x () == x2 && interesectionPoint.y () == y2 );
1894
+ bool interesectionAtReshapeLineEndpoint =
1895
+ ( interesectionPoint.x () == rx1 && interesectionPoint.y () == ry1 ) ||
1896
+ ( interesectionPoint.x () == rx2 && interesectionPoint.y () == ry2 );
1897
+
1898
+ // the intersection must be at the begin/end of both lines
1899
+ if ( interesectionAtOrigLineEndpoint && interesectionAtReshapeLineEndpoint )
1900
+ {
1901
+ GEOSGeometry* g1 = GEOSGeom_clone_r ( geosinit.ctxt , line1 );
1902
+ GEOSGeometry* g2 = GEOSGeom_clone_r ( geosinit.ctxt , line2 );
1903
+ GEOSGeometry* geoms[2 ] = { g1, g2 };
1904
+ GEOSGeometry* multiGeom = GEOSGeom_createCollection_r ( geosinit.ctxt , GEOS_MULTILINESTRING, geoms, 2 );
1905
+ GEOSGeometry* res = GEOSLineMerge_r ( geosinit.ctxt , multiGeom );
1906
+ GEOSGeom_destroy_r ( geosinit.ctxt , multiGeom );
1907
+ return res;
1908
+ }
1909
+ else
1910
+ return nullptr ;
1911
+ }
1912
+
1913
+
1857
1914
GEOSGeometry* QgsGeos::reshapeLine ( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos , double precision )
1858
1915
{
1859
1916
if ( !line || !reshapeLineGeos )
1860
1917
return nullptr ;
1861
1918
1862
1919
bool atLeastTwoIntersections = false ;
1920
+ bool oneIntersection = false ;
1921
+ QgsPoint oneIntersectionPoint;
1863
1922
1864
1923
try
1865
1924
{
@@ -1869,6 +1928,16 @@ GEOSGeometry* QgsGeos::reshapeLine( const GEOSGeometry* line, const GEOSGeometry
1869
1928
{
1870
1929
atLeastTwoIntersections = ( GEOSGeomTypeId_r ( geosinit.ctxt , intersectGeom ) == GEOS_MULTIPOINT
1871
1930
&& GEOSGetNumGeometries_r ( geosinit.ctxt , intersectGeom ) > 1 );
1931
+ // one point is enough when extending line at its endpoint
1932
+ if ( GEOSGeomTypeId_r ( geosinit.ctxt , intersectGeom ) == GEOS_POINT )
1933
+ {
1934
+ const GEOSCoordSequence* intersectionCoordSeq = GEOSGeom_getCoordSeq_r ( geosinit.ctxt , intersectGeom );
1935
+ double xi, yi;
1936
+ GEOSCoordSeq_getX_r ( geosinit.ctxt , intersectionCoordSeq, 0 , &xi );
1937
+ GEOSCoordSeq_getY_r ( geosinit.ctxt , intersectionCoordSeq, 0 , &yi );
1938
+ oneIntersection = true ;
1939
+ oneIntersectionPoint = QgsPoint ( xi, yi );
1940
+ }
1872
1941
GEOSGeom_destroy_r ( geosinit.ctxt , intersectGeom );
1873
1942
}
1874
1943
}
@@ -1878,27 +1947,18 @@ GEOSGeometry* QgsGeos::reshapeLine( const GEOSGeometry* line, const GEOSGeometry
1878
1947
atLeastTwoIntersections = false ;
1879
1948
}
1880
1949
1950
+ // special case when extending line at its endpoint
1951
+ if ( oneIntersection )
1952
+ return _mergeLinestrings ( line, reshapeLineGeos, oneIntersectionPoint );
1953
+
1881
1954
if ( !atLeastTwoIntersections )
1882
1955
return nullptr ;
1883
1956
1884
1957
// begin and end point of original line
1885
- const GEOSCoordSequence* lineCoordSeq = GEOSGeom_getCoordSeq_r ( geosinit.ctxt , line );
1886
- if ( !lineCoordSeq )
1887
- return nullptr ;
1888
-
1889
- unsigned int lineCoordSeqSize;
1890
- if ( GEOSCoordSeq_getSize_r ( geosinit.ctxt , lineCoordSeq, &lineCoordSeqSize ) == 0 )
1891
- return nullptr ;
1892
-
1893
- if ( lineCoordSeqSize < 2 )
1958
+ double x1, y1, x2, y2;
1959
+ if ( !_linestringEndpoints ( line, x1, y1, x2, y2 ) )
1894
1960
return nullptr ;
1895
1961
1896
- // first and last vertex of line
1897
- double x1, y1, x2, y2;
1898
- GEOSCoordSeq_getX_r ( geosinit.ctxt , lineCoordSeq, 0 , &x1 );
1899
- GEOSCoordSeq_getY_r ( geosinit.ctxt , lineCoordSeq, 0 , &y1 );
1900
- GEOSCoordSeq_getX_r ( geosinit.ctxt , lineCoordSeq, lineCoordSeqSize - 1 , &x2 );
1901
- GEOSCoordSeq_getY_r ( geosinit.ctxt , lineCoordSeq, lineCoordSeqSize - 1 , &y2 );
1902
1962
QgsPointV2 beginPoint ( x1, y1 );
1903
1963
GEOSGeometry* beginLineVertex = createGeosPoint ( &beginPoint, 2 , precision );
1904
1964
QgsPointV2 endPoint ( x2, y2 );
0 commit comments