Shapely MultiPolygon Construction will not accept entire set of Polygons

Question

I am new to Shapely, so I am using a work request to learn a bit more about it. I need to execute a number of unions to build larger geometries from assessment neighborhoods in DC. I am curious about why a MultiPolygon does not seem to be able to be constructed from a list unless it's sliced (or why what I am doing is incorrect). More importantly, I cannot seem to get all of my polygons in there. Here is the set up.

#Open assessment neighborhood shapefile
asmt=fiona.open(data_dir+'ASNBHDPLY.shp')

#Initiate iterator
nbhd=asmt.next()

#Create container for features
nbhd_list=[]

#For each feature...
for nbhd in asmt:
    #...throw the feature in the list
    nbhd_list.append(shape(nbhd['geometry']))

From here, I attempted to generate a MultiPolygon object with both of the following choices.

nbhds=MultiPolygon(nbhd_list)

and

nbhds = MultiPolygon([shape(nbhd['geometry']) for nbhd in asmt])

The second version was tried because of the success of the technique here.

In both cases, an exception was thrown:

NotImplementedError: A polygon does not itself provide the array interface. Its rings do.

I did note, however, that both the following approaches worked:

nbhds=MultiPolygon([nbhd_list[0],nbhd_list[1]])

and

nbhds=MultiPolygon(nbhd_list[:62])

The kicker, however, is that there are 71 assessment neighborhoods and I get the above error as soon as I exceed 62.

Answer 1

look at Creating Shapely MultiPolygons from shape file MultiPolygons and I detail below the process:

1) browse through the shapefile :

• with next()

import fiona
c = fiona.open('polygons.shp')
# the result is a Python iterator, so to read the first feature of the shapefile
print c.next()
# the result is a Python dictionary (GeoJSON)format
{'geometry': {'type': 'Polygon', 'coordinates': [[(249744.23153029341, 142798.16434689672), (250113.79108725351, 142132.95714436853), (250062.62130244367, 141973.76225829343), (249607.77877080048, 141757.71205576291), (249367.77424759799, 142304.68402918623), (249367.77424759799, 142304.68402918623), (249744.23153029341, 142798.16434689672)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1), (u'couleur', 1), (u'id_class', None)])}
# if you want the geometry only
print c.next()['geometry']
{'type': 'Polygon', 'coordinates': [[(249175.78991730965, 142292.53526406409), (249367.77424759799, 142304.68402918623), (249607.77877080048, 141757.71205576291), (249014.45396077307, 141876.13484290778), (249175.78991730965, 142292.53526406409)]]}
....

but as all the iterators, this method raises a built-in StopIteration exception at the end of the iterator or end of the file:

c.next()
Traceback (most recent call last): .....   
StopIteration

• so the solution is a list (iterator without the StopIteration exception)

c = fiona.open('polygons.shp')
features = list(c)
# and the result is a list of all the features in the shapefile:
# or
for features in fiona.open('polygons.shp'):
     ....

2) Therefore to convert the polygons to a multipolygon:

from shapely.geometry import MultiPolygo
multi = []
# append the geometries to the list
for pol in fiona.open('polygons.shp'):
    multi.append(shape(pol['geometry'])
# create the MultiPolygon from the list of Polygons
Multi = MultiPolygon(multi)
print Multi.wkt
'MULTIPOLYGON (((249744.2315302934148349 142798.1643468967231456, 250113.7910872535139788 142132.9571443685272243, ..., 249870.8182051893090829 142570.3083320840960369)))'

Or in one line with a list comprehension:

Multi = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('polygons.shp')]

Note:

If there are different types of geometries in the shapefile (MultiPolygon embedded in polygons) use unary_union: "the unary union function can operate on different geometry types, not only polygons as is the case for the older cascaded unions"

from shapely.ops import unary_union
polygon1 = Polygon([(0, 0), (1, 1), (1, 0)])
polygon2 = Polygon([(1, 1), (2,2), (2, 0)])
poly  = MultiPolygon([polygon1,polygon2])
list = [polygon1,polygon2,poly]
result = unary_union(list)
print result
MULTIPOLYGON (((0.00 0.00, 1.00 1.00, 1.00 0.00, 0.00 0.00)), ((1.00 1.00, 2.00 2.00, 2.00 0.00, 1.0 1.00)))

So:

geoms =[shape(feature['geometry']) for feature in fiona.open("polygons.shp")]
result = unary_union(geoms) 

Answer 2

Thank to Marvin's answer. I had the same problem. So based on the hint, I rewrite the code to check if there is any case a shape(...) call lead to the creation of a MultiPolygon. Then this works perfectly.

polys = []
for pol in fiona.open(shapefile):
    poly = shape(pol['geometry'])
    if isinstance(poly, MultiPolygon):
        polys += poly
    else:
        polys.append(poly)
mp = MultiPolygon(polys)

Answer 3

It turns out the issue was not with the read in of the data. For some reason, two of the Polygon objects were read in as a single MultiPolygon object. Whether or not this is an issue with Shapely or the input file is not altogether clear. The resolution was to just extract and bust up the MP. The Polygon object components were just appended to the original list. So, just in case anyone else encounters this issue, validation is your friend.