I'm trying to convert a geojson MultiPolygon which contains two overlapping polygons to a shapefile, however I noticed that the resulting shape has a smaller area. Upon inspection it's because the overlapping polygon is becoming a hole in the shapefile. Do I need to change a setting or do something different to preserve the overlapping polygon? I'm using GDAL 1.10.1 and the following code:
def convert():
data = {
"type": "FeatureCollection",
"features": [
{
"geometry": {
"type": "MultiPolygon", "coordinates": [
[[[-98.7, 49.6], [-98.7, 49.7], [-98.8, 49.7], [-98.8, 49.6], [-98.7, 49.6]]],
[[[-98.74, 49.64], [-98.74, 49.66], [-98.76, 49.66], [-98.76, 49.64], [-98.74, 49.64]]]
]
},
"type": "Feature",
"properties": {}
}
]
}
geometry = []
if (data['type'] == 'FeatureCollection'):
for feature in data['features']:
geom = ogr.CreateGeometryFromJson(json.dumps(feature['geometry']))
print 'Geom Area', geom.GetArea()
geometry.append({
'geometry': geom})
driver = ogr.GetDriverByName("ESRI Shapefile")
srs = osr.SpatialReference()
res = srs.ImportFromEPSG(4326)
print 'Import Res', res
data_source = driver.CreateDataSource('out.shp')
layer = data_source.CreateLayer('shape-layer', srs, ogr.wkbMultiPolygon)
featureDefn = layer.GetLayerDefn()
layerDefn = layer.GetLayerDefn()
for data in geometry:
feature = ogr.Feature(featureDefn)
feature.SetGeometry(data['geometry'])
layer.CreateFeature(feature)
data_source.Destroy()
Viewing the geojson and resulting shapefile in QGIS, this is what I see:
Initial Geojson
Resulting Shapefile
Best Answer
The shapefile format does not know the "MultiPolygon" geometry, replaced by a Polygon geometry.
I use here Fiona and Shapely (much easier than ogr, many examples in GIS.SE)
There are 2 polygons in the MultiPolygon
Save the geometry to a shapefile
Now open the shapefile
Therefore the MultiPolygon consisting of two Polygons is replaced in the shapefile by a Polygon consisting of two parts = a Polygon with hole
Conclusion
Modify your geometry or don't use shapefiles