Polygon se solapa cuando se convierte geojson en shapefile

Question

Estoy tratando de convertir un geojson Multipolígono que contiene dos superposición de polígonos a un shapefile, sin embargo me di cuenta de que la forma resultante tiene un área más pequeña. Luego de la inspección es debido a la superposición de polígonos se está convirtiendo en un agujero en el shapefile. Necesito cambiar una configuración o hacer algo diferente para preservar la superposición de polígonos? Estoy usando GDAL 1.10.1 y el siguiente código:

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()

Visualización de la geojson y que resulta shapefile en QGIS, esto es lo que yo veo:

Inicial Geojson

Resultando Shapefile

Answer 1

El formato shapefile no conoce el "Multipolígono" de la geometría, reemplazado por un Polígono, geometría.

Yo uso aquí Fiona y bien formada (mucho más fácil que la ogr, muchos ejemplos en los SIG.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)))

Hay 2 polígonos en la Multipolígono

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))

Guardar la geometría a un 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}})

Ahora abra el archivo de forma

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)])}

Por lo tanto, la Multipolígono que consta de dos Polígonos es reemplazado en el shapefile por un Polígono, que consta de dos partes = un Polígono con agujero

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)]

Conclusión

Modificar su geometría o no utilizar los 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}})

Answer 2

La razón principal para su problema es que el Multipolígono en su fuente de datos no es válido porque ha anidado conchas. En el OGC standand de características Simples Multipolígono se define como:

6.1.14 Multipolígono

Un Multipolígono es un Multisuperficie cuyos elementos son Polígonos. El afirmaciones para Multipolígonos son como sigue.

a) Los interiores de 2 Polígonos que son elementos de un Multipolígono puede no se cruzan.

b) los límites de Los 2 Polígonos que son elementos de un Multipolígono no puede "cruz" y puede tocar en sólo un número finito de Puntos.

No puede definir su geometría como un Multipolígono pero usted debe ajustar dentro de GeometryCollection. Como WKT válido de presentación sería

GEOMETRYCOLLECTION (
    POLYGON ((
            -98.7 49.6, 
            -98.8 49.6, 
            -98.8 49.7, 
            -98.7 49.7, 
            -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
        )))