多邊形交點錯誤| Python Shapely

我一直在試圖調試這個問題，但無法做到這一點。我正在嘗試找到兩個對象Polygon的交集。它大部分時間工作，但對於下列情況，它會產生以下例外情況：多邊形交點錯誤| Python Shapely

P1 area: 13.125721955 
P2 area: 1.0 
Traceback (most recent call last): 
File "geom2d.py", line 235, in <module> 
print p1.intersection(p2) 
File "/usr/local/lib/python2.7/dist-packages/shapely/geometry/base.py", line 334, in  intersection 
return geom_factory(self.impl['intersection'](self, other)) 
    File "/usr/local/lib/python2.7/dist-packages/shapely/topology.py", line 47, in __call__ 
    "The operation '%s' produced a null geometry. Likely cause is invalidity of the geometry %s" % (self.fn.__name__, repr(this))) 
shapely.geos.TopologicalError: The operation 'GEOSIntersection_r' produced a null  geometry. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon  object at 0x8e5ad6c>

代碼如下。

from shapely.geometry import Point,Polygon,MultiPolygon 

poly1 = [(35.0041000000000011, -88.1954999999999956), (34.9917999999999978,   -85.6068000000000069), (32.8404000000000025, -85.1756000000000029), (32.2593000000000032, -84.8927000000000049), (32.1535000000000011, -85.0341999999999985), (31.7946999999999989, -85.1358000000000033), (31.5199999999999996, -85.0438000000000045), (31.3384000000000000, -85.0836000000000041), (31.2092999999999989, -85.1069999999999993), (31.0023000000000017, -84.9943999999999988), (30.9953000000000003, -87.6008999999999958), (30.9422999999999995, -87.5926000000000045), (30.8538999999999994, -87.6256000000000057), (30.6744999999999983, -87.4072000000000031), (30.4404000000000003, -87.3687999999999931), (30.1463000000000001, -87.5240000000000009), (30.1545999999999985, -88.3863999999999947), (31.8938999999999986, -88.4742999999999995), (34.8937999999999988, -88.1020999999999930), (34.9478999999999971, -88.1721000000000004), (34.9106999999999985, -88.1461000000000041)] 
poly2 = [(34.7998910000000024, -88.2202139999999986), (34.7998910000000024, -87.2202139999999986), (35.7998910000000024, -87.2202139999999986), (35.7998910000000024, -88.2202139999999986)] 

p1 = Polygon(poly1) 
p2 = Polygon(poly2) 
print 'P1 area:',p1.area 
print 'P2 area:',p2.area 
print p1.intersection(p2)

由於它打印兩個多邊形的面積，我假設多邊形形成正確。我也（不知何故）打印第一個多邊形以確保它確實是一個簡單的多邊形。

任何人都可以請解釋爲什麼我得到這個異常？

編輯：我打印p1.is_valid，結果是False。有一些解釋here。搜索字符串is_valid。它說

有效的多邊形可能不具有任何重疊的外部或內部環。

有人能解釋一下這是什麼意思，如果有可能的解決方法嗎？順便說一句，我也注意到，如果我從poly1刪除最後一個座標，整個事情就起作用了。也許整個座標列表使多邊形複雜。

來源

2012-10-25 Nik

你會得到這個例外，因爲p1不是有效的多邊形。

>>> p1.is_valid 
False 
>>> p2.is_valid 
True

的documentation說：

有效的多邊形可能不具有任何重疊的外部或內部環。

請記住，由於多邊形的第一個和最後一個點是不同的，所以會將第一個點追加到列表的末尾。

>>> list(p1.exterior.coords) 
[(35.004100000000001, -88.195499999999996), (34.991799999999998, -85.606800000000007), (32.840400000000002, -85.175600000000003), (32.259300000000003, -84.892700000000005), (32.153500000000001, -85.034199999999998), (31.794699999999999, -85.135800000000003), (31.52, -85.043800000000005), (31.3384, -85.083600000000004), (31.209299999999999, -85.106999999999999), (31.002300000000002, -84.994399999999999), (30.9953, -87.600899999999996), (30.942299999999999, -87.592600000000004), (30.853899999999999, -87.625600000000006), (30.674499999999998, -87.407200000000003), (30.4404, -87.368799999999993), (30.1463, -87.524000000000001), (30.154599999999999, -88.386399999999995), (31.893899999999999, -88.474299999999999), (34.893799999999999, -88.102099999999993), (34.947899999999997, -88.1721), (34.910699999999999, -88.146100000000004), (35.004100000000001, -88.195499999999996)]

您的多邊形的外觀是一個線性環，這也似乎是無效的：

>>> p1.exterior.type 
'LinearRing' 
>>> p1.exterior.is_valid 
False

你也可以看到，如果你把多邊形的外觀成線串它將是複雜的：

>>> l1 = LineString(p1.exterior.coords) 
>>> l1.is_simple 
False

不知何故多邊形的外部交叉或觸及自身。

挖掘到更多的數據位，我們可以想像它在地圖上：

>>> import cgpolyencode 
>>> encoder = cgpolyencode.GPolyEncoder() 
>>> encoder.encode((y, x) for x, y in p1.exterior.coords) 
{'points': '[email protected][email protected][email protected]@[email protected]|NfjI{[email protected]`[email protected][email protected]@[email protected]~h]{[email protected]~lgDsurIjdPk|hQgugAaqIntLlgFoaDwfQvsH', 'numLevels': 18, 'zoomFactor': 2, 'levels': 'PPLMKMKGKPNIKLMNPLLKJP'}

如果你把它插入到Google's Polyline Encoder，你可以看到它是美國阿拉巴馬州。如果你放大阿拉巴馬州的左上角，你會發現兩條線相互交叉。要解決這個問題，您需要刪除最後一個點或將第二個點與第二個點交換到最後一個點。

來源

2012-10-25 06:05:49

有沒有辦法處理_complex_多邊形？ – Nik

理想情況下，您可以將所有多邊形定義爲正確的開頭。既然這是對付美國的一個國家，而不是某種線索，我會認爲修正這將是更好的路線。但是，如果這是不可能的，你總是可以使用'p1.convex_hull'而不是'p1'。 'convex_hull'永遠是有效的，永不復雜。 –

如前所述，p1無效。在繪製它時，我注意到右下方有一個「領結」。我假設你不需要你的多邊形;如果沒有，你可以嘗試身材勻稱的buffer(0)絕招（勻稱手冊中記錄）來解決這個問題：

In [382]: p1.is_valid 
Out[382]: False 

In [383]: p1 = p1.buffer(0) 

In [384]: p1.is_valid 
Out[384]: True

buffer(0)具有以下作用：

前：

enter image description here

後：

enter image description here

你現在可以這樣做：

print p1.intersection(p2) 
POLYGON ((34.9396324323625151 -88.1614025927056559, 34.8937999999999988 -88.1020999999999930, 34.7998910000000024 -88.1137513649788247, 34.7998910000000024 -87.2202139999999986, 34.9994660069532983 -87.2202139999999986, 35.0041000000000011 -88.1954999999999956, 34.9396324323625151 -88.1614025927056559))

請注意，我有一些情況下（與看上去更像是「鳥巢」不是簡單的領結區），其中並沒有工作;檢查確保你找回一個單一的Polygon對象而不是MultiPolygon。

來源

2012-12-30 19:46:42 subnivean

多邊形交點錯誤| Python Shapely

回答

相關問題