[GIS] Overlapping Polygon becoming hole when converting geojson to shapefile

Question

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

Answer 1

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)

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": {}}]}
from shapely.geometry import shape, mapping
multi =  shape(data['features'][0]['geometry'])
print multi.wkt
MULTIPOLYGON (((-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)))

There are 2 polygons in the MultiPolygon

for poly in multi:
     print poly
POLYGON ((-98.7 49.6, -98.7 49.7, -98.8 49.7, -98.8 49.6, -98.7 49.6))
POLYGON ((-98.74 49.64, -98.74 49.66, -98.76 49.66, -98.76 49.64, -98.74 49.64))

Save the geometry to a shapefile

import fiona
# schema of the shapefile
schema = {'geometry': 'MultiPolygon','properties': {'id': 'int'}}
geom = data['features'][0]['geometry']
print geom
{'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]]]]}
# save to a shapefile
with fiona.open('multipol.shp', 'w', 'ESRI Shapefile', schema)  as output:
     output.write({'geometry':geom,'properties': {'id':1}})

Now open the shapefile

multi = fiona.open('multipol.shp')
# first feature
first = multi.next()
print first
{'geometry': {'type': 'Polygon', 'coordinates': [[(-98.7, 49.6), (-98.8, 49.6), (-98.8, 49.7), (-98.7, 49.7), (-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', 'id': '0', 'properties': OrderedDict([(u'id', 1)])}

Therefore the MultiPolygon consisting of two Polygons is replaced in the shapefile by a Polygon consisting of two parts = a Polygon with hole

coords =  first['geometry']['coordinates']
for coord in coord:
     print coord
 [(-98.7, 49.6), (-98.8, 49.6), (-98.8, 49.7), (-98.7, 49.7), (-98.7, 49.6)]
 [(-98.74, 49.64), (-98.74, 49.66), (-98.76, 49.66), (-98.76, 49.64), (-98.74, 49.64)]

Conclusion

Modify your geometry or don't use shapefiles

 schema = {'geometry': 'Polygon','properties': {'id': 'int'}}
 # save to a shapefile
  with fiona.open('polys.shp', 'w', 'ESRI Shapefile', schema)  as output:
      # split multipolygon
      for poly in multi:
           output.write({'geometry':mapping(poly),'properties': {'id':1}})