13

この問題をデバッグしようとしましたが、できません。Polygon2 つのオブジェクトの交点を見つけようとしています。ほとんどの場合は機能しますが、次の場合は次の例外が発生します。

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)

2 つのポリゴンの面積が出力されるので、ポリゴンは正しく形成されていると思います。また、(どういうわけか) 最初のポリゴンを印刷して、それが実際に単純なポリゴンであることを確認しました。

この例外が発生する理由を誰か説明してもらえますか?

編集: p1.is_valid を印刷したところ、False であることが判明しました。ここにいくつかの説明があります。文字列 を検索しますis_valid。それは言う

有効な Polygon は、重複する外部リングまたは内部リングを持つことはできません。

誰かがこれが何を意味するのか、そして可能な回避策があるかどうかを説明してもらえますか? ところで、最後の座標を から削除すると、すべてが機能することにも気付きましたpoly1。おそらく、座標のリスト全体が多角形を複雑にしています。

4

2 に答える 2

26

前述のとおり、p1無効です。プロットすると、右下に小さな「ボウタイ」があることに気付きました。ポリゴンではこれは必要ないと思います。そうでない場合は、Shapely のbuffer(0)トリック (Shapely マニュアルに記載) を試して修正できます。

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

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

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

buffer(0)次の効果があります。

前：

ここに画像の説明を入力

後：

ここに画像の説明を入力

そして、これを行うことができます：

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-30T19:46:42.157 に答える
16

p1は有効なポリゴンではないため、この例外が発生しています。

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

ドキュメントには次のように記載されています。

有効な Polygon は、重複する外部リングまたは内部リングを持つことはできません。

ポリゴンの最初と最後の点は形状が異なるため、最初の点をリストの最後に追加することに注意してください。

>>> 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': 'svstEzthyOzkAkrxNfecL_fsAznpBcgv@ftSjsZnaeA~yRzst@_~P~mb@vwFzeXfqCvlg@w~Tvj@ra|NfjI{r@ngPfmEf`b@_ti@bvl@_oFbmx@~h]{r@~lgDsurIjdPk|hQgugAaqIntLlgFoaDwfQvsH', 'numLevels': 18, 'zoomFactor': 2, 'levels': 'PPLMKMKGKPNIKLMNPLLKJP'}

これをGoogle の Polyline Encoderに接続すると、アラバマ州であることがわかります。アラバマ州の左上部分を拡大すると、2 本の線が交差していることがわかります。これを修正するには、最後のポイントを削除するか、最後のポイントを最後から 2 番目のポイントと交換する必要があります。

于 2012-10-25T06:05:49.820 に答える