7

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.

3 Answers 3

5

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) 
5
  • I appreciate you taking the time to respond, and the mapping business will be helpful as well. However, I am afraid I already linked to your earlier response in my question. The proposed technique, while effective generally, was not applicable to my situation. The problem was a single MultiPolygon was embedded amongst the Polygon objects. Jan 15, 2014 at 19:06
  • see my new note in the answer.
    – gene
    Jan 15, 2014 at 19:28
  • Now that is good to know! A quick follow up: I have been working on creating unions across specific neighborhoods so that we can identify sub-District regions. I was going to try and do this by matching indices across lists (one list of neighborhood names, and one list of neighborhood Polygon objects). Is there a more efficient way of going about such a task? The groups are basically chosen by visual reference. Jan 15, 2014 at 20:26
  • For this, try rtree with shapely.
    – gene
    Jan 15, 2014 at 21:11
  • Well this does look useful. I'll have to dive in on this material. Thanks for the heads up. Jan 15, 2014 at 22:13
3

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)
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.