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.
Since celticflute did not post his code, I had to figure it out again.
Finally got a working sample. This seems to work pretty well.
The key thing is that you do not write polygon parts and holes differently!
Instead, you construct a polygon from nested arrays which represent parts made of up rings. When you write the geometry by writing it to a shape field or use it in a tool, behind the scenes arc objects planarizes the rings within in each polygon part and determines what's a hole and what isn't. And sorts the coordinates clockwise, and burns the xys into the coordinate system and domain of the geodatabase feature class.
#
# Write a polygon feature class
#
import os
import arcpy
from arcpy import env
def makepoly(coord_list, SR=None):
"""Convert a Python list of coordinates to an ArcPy polygon feature
Author: Curtis Price, USGS, cprice@usgs.gov
Examples, from Desktop Help 10.x: Reading Geometries
Feat0 = [
[[3.0, 8.0],
[1.0, 8.0],
[2.0, 10.0],
[3.0, 8.0]]
]
Feat1 = [
[[5.0, 3.0],
[3.0, 3.0],
[3.0, 5.0],
[5.0, 3.0]],
[[7.0, 5.0],
[5.0, 5.0],
[5.0, 7.0],
[7.0, 5.0]],
]
# this feature has an interior ring (donut)
Feat2 = [
[[9.0, 11.0],
[9.0, 8.0],
[6.0, 8.0],
[6.0, 11.0],
[9.0, 11.0],
None,
[7.0, 10.0],
[7.0, 9.0],
[8.0, 9.0],
[8.0, 10.0],
[7.0, 10.0]]
]
"""
parts = arcpy.Array()
rings = arcpy.Array()
ring = arcpy.Array()
for part in coord_list:
for pnt in part:
if pnt:
ring.add(arcpy.Point(pnt[0], pnt[1]))
else:
# null point - we are at the start of a new ring
rings.add(ring)
ring.removeAll()
# we have our last ring, add it
rings.add(ring)
ring.removeAll()
# if we only have one ring: remove nesting
if len(rings) == 1:
rings = rings.getObject(0)
parts.add(rings)
rings.removeAll()
# if single-part (only one part) remove nesting
if len(parts) == 1:
parts = parts.getObject(0)
return arcpy.Polygon(parts, SR)
# test data from:
# Desktop Help 10.0: Reading Geometries
Feat0 = [
[[3.0, 8.0],
[1.0, 8.0],
[2.0, 10.0],
[3.0, 8.0]]
]
Feat1 = [
[[5.0, 3.0],
[3.0, 3.0],
[3.0, 5.0],
[5.0, 3.0]],
[[7.0, 5.0],
[5.0, 5.0],
[5.0, 7.0],
[7.0, 5.0]],
]
# this last feature has an interior ring (donut)
Feat2 = [
[[9.0, 11.0],
[9.0, 8.0],
[6.0, 8.0],
[6.0, 11.0],
[9.0, 11.0],
None,
[7.0, 10.0],
[7.0, 9.0],
[8.0, 9.0],
[8.0, 10.0],
[7.0, 10.0]]
]
# test code
# create the empty feature class
# with real data, provide a SR code, name or dataset for SR
# SR = arcpy.SpatialReference(4326)
SR = None
env.workspace = env.scratchGDB
Data = arcpy.CreateScratchName("","","featureclass",env.workspace)
print "writing: " + Data
print
arcpy.CreateFeatureclass_management(os.path.dirname(Data),
os.path.basename(Data),"Polygon",
spatial_reference=SR)
# create the polygons and write them
Rows = arcpy.da.InsertCursor(Data, "SHAPE@")
for f in [Feat0, Feat1, Feat2]:
print "coords: " + repr(f)
p = makepoly(f)
print "feature: " + repr(p)
Rows.insertRow([p])
del Rows
Best Answer
Quoting the shapefile specification:
So, by definition, a polygon must have at least one exterior ring. Interior rings (aka "holes") can only exist within a parent exterior ring. If there are multiple exterior rings, they can touch each other at a single point, but they cannot share a boundary or overlap. Interior rings must be wholly contained within their parent rings (they cannot touch the exterior ring at all), and can only touch other interior rings at a single point (again, no shared segment or overlap). Interior rings are distinguished from exterior rings by vertex order (clockwise for exterior and counter-clockwise for interior, as per "right hand rule"). There is no upper limit to the number of exterior rings, or to the number of interior rings within any given exterior ring (though file size limitations would impose an upper limit if you work it backwards, four points per ring, sixteen bytes per point -- roughly 33.5 million rings).