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.
I don't really understand your problem. Some explanations: for example, we want to
- create a new shapefile
modify the original schema: for that, it's easier to just copy things to a new shapefile and make the changes as you copy
- Create a new Polyline shapefile:
import fiona
# schema: it is a simple dictionary with geometry and properties as keys
schema = {'geometry': 'LineString','properties': {'test': 'int'}}
# for defining the geometry, you need Shapely
from shapely.geometry import LineString, mapping
# two simples geometries
lines = [LineString([(272830.63,155125.73),(273770.32,155467.75)]),LineString([(273536.47,155914.07),(272033.12,152265.71)])]
with fiona.open('myshp.shp', 'w', 'ESRI Shapefile', schema) as layer:
for line in lines:
# filling schema
elem = {}
# geometry with mapping function of shapely
elem['geometry'] = mapping(line)
# attribute value (the same here)
elem['properties'] = {'test': 145}
# writing element in the file
layer.write(elem)
- Now we want to modify the schema of the original shapefile in a new shapefile:
1) Open the original shapefile:
shapefile =fiona.open('myshp.shp')
#read the schema
schema2 = shapefile.schema
print schema2
{'geometry': 'LineString', 'properties': OrderedDict([(u'test', 'int:10')])}
2) As it is a dictionary, it is easy to add new fields/keys in the properties:
schema2['properties']['string']='str'
3) Now we create a new shapefile copying myshp.shp with the new schema :
with fiona.open('myshp.shp', 'r') as input:
schema = schema2
# writing the new shapefile
with fiona.open('myshp_copy.shp', 'w', 'ESRI Shapefile', schema) as output:
for elem in input:
# add the new attribute value
elem['properties']['string']="hello"
output.write({'properties': elem['properties'],'geometry': mapping(shape(elem['geometry']))})
c = fiona.open('myshp_copy.shp')
c.schema
{'geometry': 'LineString', 'properties': OrderedDict([(u'test', 'int:10'), (u'string', 'str')])}
# first element of the shapefile
c.next()
{'geometry': {'type': 'LineString', 'coordinates': [(272830.63, 155125.73000000001), (273770.32000000001, 155467.75)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'test', 145), (u'string', u'hello')])}
If you don't want to modify the shapefile, it is easier with schema.copy() that is only used to get a copy of the original schema (no definition here)
with fiona.open('myshp.shp', 'r') as input:
schema = input.schema.copy()
with fiona.open('myshp_copy2.shp', 'w', 'ESRI Shapefile', schema) as output:
for elem in input:
output.write({'properties': elem['properties'],'geometry': mapping(shape(elem['geometry']))})
Best Answer
The Shapely method object.within(other) looks tailored for your purpose.
You may also be interested in object.touches(other).