# 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.

• This question is missing the crucial information as to what you use to read the Shapefile. – inc42 Nov 13 '19 at 19:32
• @inc42 I was using python to read the file directly. – BenjaminGolder Nov 14 '19 at 15:53
• Ah, then the Shapely bit is misleading. Your actual problem was finding out how to determine the "kind" of ring in the Shapefile format. :) – inc42 Nov 16 '19 at 17:51

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.

• thanks, this isn't what I'm looking for, but I saw some more clarifications I could make to my question. To answer your questions: 1, because I'm sorting them to do tons of intersections in the first place. 2, nothing, I just need one that is lightweight and that I can use easily from python, and I haven't learned ogr. 3, no, it wouldn't be easier in this case. – BenjaminGolder Jan 24 '11 at 10:54
• wow! Thanks geographika! I'm now 90% convinced that ogr is the solution here. I ran into some installation issues with ogr, but it looks like they're worth working through. So ogr organizes the rings in its wkt output? and is it fine with 3d shapefiles (I use a lot of 3d)? – BenjaminGolder Jan 24 '11 at 11:26
• It looks like the answers to my questions are yes and yes. And thank you for introducing me to rtree. OGR appears to solve my problem entirely - as long as I can configure the installation correctly this time. – BenjaminGolder Jan 24 '11 at 11:44
• Re installation of OGR - remember you'll need the binaries for Windows and the matching Python bindings. – geographika Jan 24 '11 at 12:05

First, use ogr to open shapefile:

``````from osgeo import ogr
source = ogr.Open("mpolys.shp")
layers =  source.GetLayerByName("mpoly")
len(layers)
1
``````

convert shapefile geometries into shapely geometries

``````from shapely.wkb import loads
element=layers[0] #(because lenght of layer =1, else you need "for element in layers: ...")
geom.geom_type
'MultiPolygon'
print geom
MULTIPOLYGON ((..... # the geometry in shapely wkt format
``````

For the polygons in the multipolygon:

``````poly=[]
for pol in geom:
poly.append(pol)
poly[0]
<shapely.geometry.polygon.Polygon object at 0x00B82CB0>
poly[0].geom_type
'Polygon'
print(poly[1])
POLYGON ((.... # the geometry in shapely wkt format
``````

And now, you can use all the functions of shapely (shapely)

I'm not too familiar with how polygons are actually stored in shape files, but - shouldn't a polygon ring be a closed loop if and only if the start coordinate is repeated? So if you compare each subsequent coordinate with the start coordinate you will find the first point where a polygon is closed. If that is the last coordinate of the polygon, it is a simple polygon, if not, it is a multipolygon and requires processing the other loops.

That might be the 'more intensive processing' that you want to avoid, but it is really just an iteration through the coordinates that comes for free when you have to read them in anyway.

• I apologize - my question was a bit misleading and I've edited it. I have the polygon rings, and I know when I have a list of more than one polygon ring, but aside from the clockwise or counterclockwise order of the points, I do not know if they are interior or exterior rings. If they are interior rings, I have no way of being sure which exterior ring they belong to, apart from measuring their location. – BenjaminGolder Jan 24 '11 at 10:44
• I see. I'm also happy to see that you found your way to ogr. ;) – relet Jan 24 '11 at 11:48

As indicated by @inc42, on page 8 of the ESRI Shapefile Technical Description: An ESRI White Paper--July 1998 (page 12 of 34 in the PDF) there is a discussion of polygon record contents which may be what you are after:

A polygon may contain multiple outer rings. The order of vertices or orientation for a ring indicates which side of the ring is the interior of the polygon.

• p.10 has it "A polygon may contain multiple outer rings. The order of vertices or orientation for a ring indicates which side of the ring is the interior of the polygon." – inc42 Nov 13 '19 at 19:32
• @inc42 updated accordingly although I think the page number should be 8 (or 12 of 34). – PolyGeo Nov 13 '19 at 22:18