Monday, 23 March 2015

PyQgis: inaccurate points from geom.splitGeometry()


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)




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...