PyQgis: inaccurate points from geom.splitGeometry()

Question

I have some problems with inaccurate points from geom.splitGeometry() and wonder if this is a known bug (or even an intended mechanism), if it is possible to prevent it or work around it and if it is worth writing a bug report (or if it is a known and accepted bug e.g. cause of some rounding errors, that can't be prevented).

Summarizing the problem:

• find intersection point P of two lines, A and B
• distance from P to A is zero, and intersects
• distance from P to B is not zero, and doesn't intersect
• distance from P to backup of A is zero but does not intersect
• distance from A to B is also zero but does not intersect
• happens also with geom.intersects(), geom.closestSegmentWithContext(), geom.nearestPoint() and probably further
• is a big problem e.g. if using split Points of Polygons for creating new Lines and testing spatial operators with original Polygon

Some further remarks regarding splitGeometry():

Especially I don't exactly know what the 2nd statement in splitGeometry does (something with topological editing and if I turn it to False, I don't even get a point back - I am not quite sure if it is turned on standardly in a standalone script or how to test/change it and if it would change anything - maybe "true if topological editing is enabled" may also mean that with True it is enabled for this process? I am not quite sure).

To clarify some very simple code:

#creting points and lines and copies of lines
a = QgsPoint(0, 1)
b = QgsPoint(1, 10)
c = QgsPoint(-5, 6)
d = QgsPoint(4, 5)
line_a = QgsGeometry().fromPolyline([a, b])
line_a_backup = QgsGeometry(line_a)
line_b = QgsGeometry().fromPolyline([c, d])

#splitting one line with the other
succes, geoms, points = line_a.splitGeometry([c, d], True)

#testing returned point on intersection with lines
point = points[0]
print point
>>>(0.487805,5.39024)

qgsp = QgsGeometry().fromPoint(point)
if qgsp.intersects(line_a):
    print 'intersects'
else:
    print qgsp.distance(line_a)
>>>intersects
#all good with the splitted line

if QgsGeometry().fromPoint(point).intersects(line_b):
    print 'intersects'
else:
    print qgsp.distance(line_b)
>>>3.92331593257e-16
#distance between the split line and the intersection point with the splitted line? 
#Though very small is a big problem with spatial operators

if qgsp.intersects(line_a_backup):
    print 'intersects'
else:
    print qgsp.distance(line_a)
>>>0.0
#Still no intersection even with Zero distance?

if line_a.intersects(line_b):
    print 'intersects'
else:
    print line_a.distance(line_b)
>>>0.0
#Still no intersection even with Zero distance?

This is very frustrating and expensive to deal with if working with polygons - e.g. I want to test if one of the split points from two seperate splitGeometrys are in the specific subpolygon or test if a line beginning at the split point crosses the original polygon - which might happen only through these deviations, while it actually doesn't or shouldn't at least.

(tested on Qgis 2.18.3 and Windows and Qgis 2.8.6 and Ubuntu and Qgis 2.8.1 and Scientific Linux)

Answer 1

I could reproduce this, both with your code, and using another approach (QGIS 2.18.10/Ubuntu Xenial 32 bit)

This is very strange behaviour. To summarise, your problem is

• find intersection point P of two lines, A and B
• distance from P to A is zero, and intersects
• distance from P to B is not zero, and doesn't intersect

I tried this using a different approach - to rule out a problem with the splitGeometry() implementation.

In this case I find the intersection of your two linestrings, then measure this point's distance to each of the two line strings. This gave the same results as your splitGeometry approach.

>> inter = QgsGeometry.intersection(line_a,line_b)
>> print inter.exportToWkt()
Point (0.48780487804878048 5.39024390243902474)
>> inter.intersects(line_a)
True
>> inter.intersects(line_b)
False
>> inter.distance(line_a)
0.0
>> inter.distance(line_b)
3.9233159325692157e-16

(In my case I declared the geometries using floats rather than integers; it doesn't seem to make a difference).

This might not be a QGIS issue; QGIS delegates these calculations to the GEOS library.

I wonder if this might be best asked on the GEOS-dev mailing list ?

I wonder if this might be an obscure floating point precision or rounding issue. I'm not a C developer though, so I struggle to follow C code, and I'm not able to test the GEOS library directly :-(

Answer 2

This is due to the limitations of finite-precision floating point arithmetic (similar to the well-known issues with round-off when doing arithmetic calculation).

It's usually not possible to represent the exact value of the intersection point of two line segments with double-precision endpoints, using double-precision numbers. Also, there is usually round-off error in the calculation of the intersection point.

You will probably need to reframe your problem. Sometimes you can use a distance computation with a tolerance value to decide if a point is close enough to a line to be considered collinear.