How to convert NumPy arrays obtained from cv2.findContours to Shapely polygons?

Question

I am using CV2 to find contours from an image and then converting them into polygons using Shapely. I am currently stuck because when I try putting one of the contour arrays into Polygon() from Shapely it throws an unspecified error.

I have double-checked that I imported everything I needed, and that creating a Shapely polygon works when I manually enter the array coordinate points.

Here is the problematic section of the code:

contours, hierarchy = cv2.findContours(thresh, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
testcontour = contours[1]

ply = Polygon(testcontour)

Where the list of contours looks like this:

contours = [np.array([[[700, 700]],
                      [[700, 899]],
                      [[899, 899]],
                      [[899, 700]]]), 
            np.array([[[774, 775]], 
                      [[775, 774]],
                      [[824, 774]],
                      [[825, 775]],
                      [[825, 824]],
                      [[824, 825]],
                      [[775, 825]],
                      [[774, 824]]]), 
            np.array([[[200, 200]],
                      [[200, 399]],
                      [[399, 399]],
                      [[399, 200]]]), 
            np.array([[[274, 275]],
                      [[275, 274]],
                      [[324, 274]],
                      [[325, 275]],
                      [[325, 324]],
                      [[324, 325]],
                      [[275, 325]],
                      [[274, 324]]])]

The error I get is:

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-65-4124f49b42e1> in <module>
----> 1 ply = Polygon(testcontour)

~\AppData\Local\Continuum\anaconda3\envs\geocomp\lib\site-packages\shapely\geometry\polygon.py in __init__(self, shell, holes)
    238 
    239         if shell is not None:
--> 240             ret = geos_polygon_from_py(shell, holes)
    241             if ret is not None:
    242                 self._geom, self._ndim = ret

~\AppData\Local\Continuum\anaconda3\envs\geocomp\lib\site-packages\shapely\geometry\polygon.py in geos_polygon_from_py(shell, holes)
    492 
    493     if shell is not None:
--> 494         ret = geos_linearring_from_py(shell)
    495         if ret is None:
    496             return None

~\AppData\Local\Continuum\anaconda3\envs\geocomp\lib\site-packages\shapely\speedups\_speedups.pyx in shapely.speedups._speedups.geos_linearring_from_py()

AssertionError: 

Answer 1

The problem is that for some reason cv2.findContours returns each contour as a 3D NumPy array with one redundant dimension:

>>> contours[1]
array([[[774, 775]],
       [[775, 774]],
       [[824, 774]],
       [[825, 775]],
       [[825, 824]],
       [[824, 825]],
       [[775, 825]],
       [[774, 824]]])

but Shapely expects a 2D array in this form (see the docs):

array([[774, 775],
       [775, 774],
       [824, 774],
       [825, 775],
       [825, 824],
       [824, 825],
       [775, 825],
       [774, 824]])

So, what we can do is to use np.squeeze to remove that redundant dimension, and use the result to obtain our polygon:

import numpy as np
from shapely.geometry import Polygon

contour = np.squeeze(contours[1])
polygon = Polygon(contour)
print(polygon.wkt)
# POLYGON ((774 775, 775 774, 824 774, 825 775, 825 824, 824 825, 775 825, 774 824, 774 775))

---

In case if you want to convert all of the contours at once, I would do it like this:

contours = map(np.squeeze, contours)  # removing redundant dimensions
polygons = map(Polygon, contours)  # converting to Polygons
multipolygon = MultiPolygon(polygons)  # putting it all together in a MultiPolygon

The resulting multipolygon will look like this:

And to get the second polygon from here you would just write:

my_polygon = multipolygon[1]
print(my_polygon.wkt)
# POLYGON ((774 775, 775 774, 824 774, 825 775, 825 824, 824 825, 775 825, 774 824, 774 775))