python - using shapely: translating between Polygons and MultiPolygons

[EDIT: the solution to this was simply to use OGR to read shapefiles. See geographika's example.]

In an ESRI shapefile, there is no distinction between Polygons and MultiPolygons. Furthermore, there is no explicit distinction between interior holes and exterior rings (besides the "handedness" of a given polygon).

So after reading a shapefile, I have a list of coordinate sequences describing rings, but without some more intensive processing, I cannot distinguish which of these rings are exterior rings, interior holes, or additional polygons.

It appears that for shapely's Polygon and MultiPolygon constructors, there must be a clear distinction between exterior and interior rings, so how should I move from an unclear list of rings to an ordered set of separated polygons, with clearly designated interior and exterior rings?

To summarize: if I have a list of polygon rings, but I don't know which rings are holes in the interior or are separate polygons, how should I best sort them into separate polygons with designated interior holes?

I'm looking for a simple algorithmic solution that I can implement in python, can be used to process hundreds of polygons in ~a minute or less, and I'm doing this in order to perform a large number of intersections.

Answer

Further to relet's answer on how to get individual polygons, you can then run an intersection on all the polygons to create the holes. If your dataset contains overlapping polygons though you're out of luck.

Explain again what is wrong with existing shapefile readers?

Would it not be easier to export feature IDs and M values from the shapefile and then join them back to the polygons after using an existing shapefile reader?

For multipatches you can use the same technique of assigning polygon IDs to a "patch ID" and then adding this attribute back to the features.

Edit: Whilst you say you don't want to use OGR, just in case you change your mind..

import ogr
# Get the driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# Open a shapefile
shapefileName = "D:/temp/myshapefile.shp"
dataset = driver.Open(shapefileName, 0)

layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):

    feature = layer.GetFeature(index)
    geometry = feature.GetGeometryRef()
    #geometry for polygon as WKT, inner rings, outer rings etc. 
    print geometry

The geometry should be output as follows:

POLYGON ((79285 57742,78741 54273...),(76087 55694,78511 55088,..))

The first bracket contains the coords of the exterior ring, subsequent brackets the coords of interior rings. If you have Z values points should be in the format 79285 57742 10 (where the last coord is a height).

Otherwise you could use the Shapely Contains and Within functions to assess every polygon with each other and apply a spatial index beforehand - http://pypi.python.org/pypi/Rtree/ to speed up processing.