减去两个多边形之间的重叠时剩余点不正确

Question

我识别相交的多边形，然后计算它们的交集并将它们从第一个多边形中删除。我得到了一个奇怪的剩余点。如果我只考虑多边形之间的差异，这种情况就不会发生。

为什么？

import pandas as pd
import geopandas as gpd
import shapely

test = pd.DataFrame({"geo_uuid": ["a", "b"], "geometry": [shapely.from_wkt("POLYGON ((-85.489530540446 40.72197039620563, -85.54038813980434 40.686016211216575, -85.42963027177998 40.682309666676865, -85.489530540446 40.72197039620563))"),
                                                            shapely.from_wkt("POLYGON ((-85.52717584228726 40.732495186799525, -85.49542480146151 40.693215026642285, -85.44637022730953 40.73203119983377, -85.52717584228726 40.732495186799525))")]})

test = gpd.GeoDataFrame(test).set_geometry("geometry").set_crs(epsg=4326)

geo_1 = test.geometry.values[0]
geo_2 = test.geometry.values[1]

overlap = geo_1.intersection(geo_2)
geo_1_overlap_removed = geo_1.difference(overlap)  # this will have an extra point contained in geo_2

geo_1_minus_geo_2 = geo_1.difference(geo_2)  # this is good

为什么会出现这种情况？我做错了什么？

Answer 1

这是因为坐标没有无限精度。

为了说明这一点，这里绘制了不同多边形的图...

当您比较第一个图和第二个图时，您可以看到第二个图的最左边的点是计算为原始多边形的两条边界线的交点的点。

但是，由于计算机不使用无限精度，因此在相当多的情况下，这个计算出的点不会精确地位于交点上，而是在其旁边大约一纳米处。在这个例子中，它最终在它的右侧〜一纳米，所以当计算差异时，一个非常细的三角形不会被删除。

处理此问题的一个选项是使用 shapely.set_ precision 或 geopandas.set_ precision。这些函数会将所有坐标舍入到指定的精度+将删除多边形中比大致此精度窄的所有部分。结果见下面的第四张图。

下面的代码示例展示了如何使用它。

代码示例：

from matplotlib import pyplot as plt
import pandas as pd
import geopandas as gpd
import shapely
import shapely.plotting as plotter

test = pd.DataFrame({"geo_uuid": ["a", "b"], "geometry": [shapely.from_wkt("POLYGON ((-85.489530540446 40.72197039620563, -85.54038813980434 40.686016211216575, -85.42963027177998 40.682309666676865, -85.489530540446 40.72197039620563))"),
                                                            shapely.from_wkt("POLYGON ((-85.52717584228726 40.732495186799525, -85.49542480146151 40.693215026642285, -85.44637022730953 40.73203119983377, -85.52717584228726 40.732495186799525))")]})

test = gpd.GeoDataFrame(test).set_geometry("geometry").set_crs(epsg=4326)

geo_1 = test.geometry.values[0]
geo_2 = test.geometry.values[1]

overlap = geo_1.intersection(geo_2)
geo_1_overlap_removed = geo_1.difference(overlap)  # this will have an extra point contained in geo_2
geo_1_overlap_removed_prec = shapely.set_precision(geo_1_overlap_removed, grid_size=1e-9)

geo_1_minus_geo_2 = geo_1.difference(geo_2)  # this is good

print(geo_1_overlap_removed)
print(geo_1_minus_geo_2)

fig, ax = plt.subplots(ncols=5, sharex=True, sharey=True)
plotter.plot_polygon(geo_1, ax=ax[0])
plotter.plot_polygon(geo_2, ax=ax[0], color="red")
plotter.plot_polygon(overlap, ax=ax[1], color="red")
plotter.plot_polygon(geo_1_overlap_removed, ax=ax[2], color="orange")
plotter.plot_polygon(geo_1_overlap_removed_prec, ax=ax[3], color="orange")
plotter.plot_polygon(geo_1_minus_geo_2, ax=ax[4], color="green")

plt.show()