如何处理Shapely中的舍入误差

Question

我有一个案例，基于将一个点投影到一条线上，然后在这条线上将其分离。我的用例稍微复杂一些，但我的问题可以通过以下代码进行复现：

from shapely import *
line1 = LineString([(1,1.2), (2,2), (3, 2.), (4,1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))

按照构造，"pr" 应该在第一行，它们的交点也应在第一行：

line1.contains(pr)
line1.intersects(LineString([pt, pr]))

打印两次“True”。但是稍微改变输入的坐标会导致工作流程中断：

from shapely import *
line1 = LineString([(1,1.2), (2,2), (3, 2.3), (4,1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))
line1.contains(pr)
line1.intersects(LineString([pt, pr]))

打印出“False”。

我明白这背后的浮点精度问题，但这是否意味着我永远不能测试点是否在直线上？当我基于点列表构建一条直线时，我能够确定至少所有“构建”点都在该直线上吗？

- Fabzi

---

你能否选择更细粒度的单位，比如毫米而不是米？ - Paulo Scardine

@PauloScardine：谢谢。是的，如果我能获得稳定性，我可以轻松放弃精度。将我的值乘以10就可以解决问题。但这在所有情况下都有效吗？Shapely在内部继续使用浮点数。 - Fabzi

Answer 1

基本上，需要一个精度模型，并且有各种计划在将来将其实现到GEOS中（不要抱太大希望，因为这已经讨论了好几年）。

否则，选项是基于距离的测试（推荐）或更昂贵的缓冲区技术进行微小调整（请参见机器epsilon）：

from shapely.geometry import LineString, Point

line1 = LineString([(1, 1.2), (2, 2), (3, 2.3), (4, 1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))

# Distance based
print(line1.distance(pr) == 0.0)  # True

# Buffer based
EPS = 1.2e-16
print(line1.buffer(EPS).contains(pr))  # True
print(line1.buffer(EPS).intersects(LineString([pt, pr])))  # True

您还可以使用or运算符来链接较便宜和昂贵的测试，例如：

print(line1.contains(pr) or line1.buffer(EPS).contains(pr))

如果第一个测试返回False，则仅运行第二个更昂贵的测试。

- Mike T

---

嗨@MikeT，感谢您的回答，这就是我担心的。我真的希望我们有一个精度模型。或者使用整数而不是浮点数(见上面保罗的评论)的技巧能确保稳定性吗？我的用例更复杂，基于联合: shapely.ops.polygonize(LinearRing(line1).union(LineString([pt, pr])))如果两条线相交或不相交，则返回不同的多边形......(我希望我可以给您的回答点赞，但我没有足够的声望;)) - Fabzi

这些距离测试能否用于检查两个线串是否接触 - 即其中一个的最后一个点是否是另一个的第一个点？我遇到过线串，其中点之间存在那些特征为0.0000000001（由于浮点数不准确性而导致）的差异。 - Mr_and_Mrs_D

Answer 2

从 shapely 2.0 开始（请参见发布说明），现在提供了一个精度模型：set_precision 函数和叠加函数的 grid_size 参数。

使用您之前的示例并将相同的任意精度设置为几何图形。

from shapely import set_precision
from shapely.geometry import LineString, Point

precision = 1e-10
line1 = set_precision(LineString([(1, 1.2), (2, 2), (3, 2.3), (4, 1.2)]), precision)
pt = set_precision(Point(2.5, 1.2), precision)
pr = set_precision(line1.interpolate(line1.project(pt)), precision)

print(line1.intersects(LineString([pt, pr])))  # True

注意：我本来以为对于pr对象，使用set_precision是多余的，因为所有涉及到的几何图形都是用任意精度创建的，但事实并非如此...

contains问题

从我的实验中可以看出，无论精度和精度组合如何，contains谓词都不会成立。我在GEOS/shapely方面的资格还不够，无法给出适当的解释。

尽管如此，直觉告诉我这与所涉及的几何类型有关：使用@MikeT解释的通常的缓冲技术，我们实际上将线/点几何图形更改为具有边界和内部维度影响的多边形。

- jmon12