目标:计算8
不同交点的联合 area
数,每一个与3
MultiPolygons
的交点有关。
有3
源,每个源代表相同的8
形状组。
从数学上讲,我的直觉是指Jaccard指数。
数据
我有3个MultiPolygon
列表:
extracted_multipoly
original_multipoly
wkt_multipoly
它们每一项都包括:
[<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e5a8cbb0>,
<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e319fb50>,
<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e303fe20>,
<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e30805e0>,
<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e302d7f0>,
<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e5a2aaf0>,
<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e5a2a160>,
<shapely.geometry.multipolygon.MultiPolygon at 0x7f18e5a2ae20>]
提取area
extracted_multipoly_area = [mp.area for mp in extracted_multipoly]
original_multipoly_area = [mp.area for mp in original_multipoly]
wkt_multipoly_area = [mp.area for mp in wkt_multipoly]
它们每一项都包括:
[17431020.0,
40348778.0,
5453911.5,
5982124.5,
8941145.5,
11854195.5,
10304965.0,
31896495.0]
程序尝试
使用MultiPolygon
for i, e in enumerate(extracted_multipoly):
for j, o in enumerate(original_multipoly):
for k, w in enumerate(wkt_multipoly):
if e.intersects(o) and e.intersects(w):
print(i, j, k, (e.intersection(o, w).area/e.area)*100)
[2022-11-18 10:06:40,387][ERROR] TopologyException: side location conflict at 8730 14707. This can occur if the input geometry is invalid.
---------------------------------------------------------------------------
PredicateError Traceback (most recent call last)
File ~/miniconda3/envs/pdl1lung/lib/python3.9/site-packages/shapely/predicates.py:15, in BinaryPredicate.__call__(self, this, other, *args)
14 try:
---> 15 return self.fn(this._geom, other._geom, *args)
16 except PredicateError as err:
17 # Dig deeper into causes of errors.
File ~/miniconda3/envs/pdl1lung/lib/python3.9/site-packages/shapely/geos.py:609, in errcheck_predicate(result, func, argtuple)
608 if result == 2:
--> 609 raise PredicateError("Failed to evaluate %s" % repr(func))
610 return result
PredicateError: Failed to evaluate <_FuncPtr object at 0x7f193af77280>
During handling of the above exception, another exception occurred:
TopologicalError Traceback (most recent call last)
Cell In [38], line 4
2 for j, o in enumerate(original_multipoly):
3 for k, w in enumerate(wkt_multipoly):
----> 4 if e.intersects(o) and e.intersects(w):
5 print(i, j, k, (e.intersection(o, w).area/e.area)*100)
File ~/miniconda3/envs/pdl1lung/lib/python3.9/site-packages/shapely/geometry/base.py:799, in BaseGeometry.intersects(self, other)
797 def intersects(self, other):
798 """Returns True if geometries intersect, else False"""
--> 799 return bool(self.impl['intersects'](self, other))
File ~/miniconda3/envs/pdl1lung/lib/python3.9/site-packages/shapely/predicates.py:18, in BinaryPredicate.__call__(self, this, other, *args)
15 return self.fn(this._geom, other._geom, *args)
16 except PredicateError as err:
17 # Dig deeper into causes of errors.
---> 18 self._check_topology(err, this, other)
File ~/miniconda3/envs/pdl1lung/lib/python3.9/site-packages/shapely/topology.py:37, in Delegating._check_topology(self, err, *geoms)
35 for geom in geoms:
36 if not geom.is_valid:
---> 37 raise TopologicalError(
38 "The operation '%s' could not be performed. "
39 "Likely cause is invalidity of the geometry %s" % (
40 self.fn.__name__, repr(geom)))
41 raise err
TopologicalError: The operation 'GEOSIntersects_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.multipolygon.MultiPolygon object at 0x7f18e5be2f70>
使用area
for e, o, w in zip(extracted_multipoly_area, original_multipoly_area, wkt_multipoly_area):
print(e, o, w)
print(e.intersection(o, w))
22347776.0 22544384.0 17431020.0
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In [31], line 3
1 for e, o, w in zip(extracted_multipoly_area, original_multipoly_area, wkt_multipoly_area):
2 print(e, o, w)
----> 3 print(e.intersection(o, w))
AttributeError: 'float' object has no attribute 'intersection'
https://stackoverflow.com/questions/74488188
复制相似问题