Shapely Polygons | Calculate Intersection of Union

areafor looppolygonpython 3shapely

Goal: Calculate 8 different Intersection of Union area numbers, each concerning the intersection of 3 MultiPolygons.

There are 3 sources, each representing the same 8 groups of shapes.

Mathematically, my instinct is to refer to the Jaccard Index.

Data

I've 3 MultiPolygon lists:

  1. extracted_multipoly
  2. original_multipoly
  3. wkt_multipoly

They each contain e.g.:

[<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>]

Extracting 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]

They each contain e.g.:

[17431020.0,
 40348778.0,
 5453911.5,
 5982124.5,
 8941145.5,
 11854195.5,
 10304965.0,
 31896495.0]

Procedure Attempts

Using 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>

Using 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'

Best Answer

Solution

IoU values should be between 0 to 1.

intersection_of_union = []
for e, o in zip(extracted_multipoly, original_multipoly):
    e = e.buffer(0)
    o = o.buffer(0)
    intersection_area = e.intersection(o).area
    intersection_of_union.append(intersection_area / (e.area + o.area - intersection_area))

[0.8970148657684971,
 0.9377700784370339,
 0.8136220015019057,
 0.8980586930524846,
 0.8496839666124079,
 0.8428598403182237,
 0.8599616483904042,
 0.9550894396247209]

Adapted from tutorial.

Related Question