有趣的问题。
您有机会拥有许多多边形(总共 34 个),这使得调试更容易:
multi_poly_1 = RGeo::Geos.factory(srid: 4326).parse_wkt(wkt_1)
multi_poly_2 = RGeo::Geos.factory(srid: 4326).parse_wkt(wkt_2)
polygons = multi_poly_1.each.to_a + multi_poly_2.each.to_a
可能存在有问题的多边形或多边形对,所以让我们尝试找到它们。此代码遍历索引的每个组合并检查 union 是否为 nil :
(0...polygons.size).to_a.combination(2).each do |i1, i2|
if polygons[i1] + polygons[i2] == nil then
p [i1, i2]
end
end
它返回
# [1, 11]
# [1, 12]
1 是非空多边形,11 和 12 看起来像线。
保留 1 并删除 11 和 12 是不够的:所有多边形的并集仍然为零。
可能仍然有线条或非常平坦的多边形:
polygons.reject!{|poly| poly.area < 1E-20}
p polygons.size
# 25
现在 9 个多边形消失了,可以计算并集:
sum = polygons.inject(&:union)
p sum.area - multi_poly_1.area - multi_poly_2.area
#=> -5.800481622797449e-18
面积差异不大,可能来自相交的多边形或已删除的(非常)小多边形。