readshapefile returns a multipolygon as several polygons

I wanted to draw a choropleth map using basemap. My plan was to read the polygons from the shapefile vis readshapefile and then further process them using shapely, descartes, etc.

The problem that I encountered was that readshapefile returned multiple polygons for each MultiPolygon. This causes confusion, as for the enclosed polygon we will have more than one record.

I think the fix will be to return a MultiPolygon instead of multiple polygons.

Here is an example:

from collections import OrderedDict
import fiona
from mpl_toolkits.basemap import Basemap


rec1 = dict()
rec2 = dict()

rec1['id'] = 0
rec1['type'] = 'Feature'
rec1['properties'] = OrderedDict([(u'NAME', 'outer'),(u'VALUE', 100)])
rec1['geometry'] = dict()
rec1['geometry']['type'] = 'MultiPolygon'
rec1['geometry']['coordinates'] = [[[(-50,40),(50,40),
                                    (50,-40),(-50,-40),
                                    (-50,40)
                                   ]],
                                  [[(-45,35),(45,35),
                                    (45,-35),(-45,-35),
                                    (-45,35)
                                   ]]]

rec2['id'] = 1
rec2['type'] = 'Feature'
rec2['properties'] = OrderedDict([(u'NAME', 'inner'),(u'VALUE', 1)])
rec2['geometry'] = dict()
rec2['geometry']['type'] = 'Polygon'
rec2['geometry']['coordinates'] = [[(-45,35),(45,35),
                                    (45,-35),(-45,-35),
                                    (-45,35)
                                   ]]
driver = 'ESRI Shapefile'
crs = {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'}
schema = {'geometry': 'Polygon',
         'properties': {'NAME': 'str',
                'VALUE': 'float:15.2'}}

with fiona.open(
     '/tmp/foo.shp',
     'w',
     driver=driver,
     crs=crs,
     schema=schema) as c:
    c.write(rec1)
    c.write(rec2)

m = Basemap()
m.readshapefile('/tmp/foo','foo')
for info,region  in zip(m.foo_info, m.foo):
    print info
    print region

for which the output is:

{'RINGNUM': 1, 'NAME': 'outer', 'VALUE': 100.0, 'SHAPENUM': 1}
[(-50.0, 40.0), (50.0, 40.0), (50.0, -40.0), (-50.0, -40.0), (-50.0, 40.0)]
{'RINGNUM': 2, 'NAME': 'outer', 'VALUE': 100.0, 'SHAPENUM': 1}
[(-45.0, 35.0), (-45.0, -35.0), (45.0, -35.0), (45.0, 35.0), (-45.0, 35.0)]
{'RINGNUM': 1, 'NAME': 'inner', 'VALUE': 1.0, 'SHAPENUM': 2}
[(-45.0, 35.0), (45.0, 35.0), (45.0, -35.0), (-45.0, -35.0), (-45.0, 35.0)]

Commencez au point d’entrée Basemap.readshapefile et reproduisez l’exemple Fiona fourni, en comparant les enregistrements renvoyés aux géométries sources. C’est terminé lorsque le cas MultiPolygon conserve la géométrie et les métadonnées d’une feature source sans casser le cas Polygon séparé.