6

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

Initial Geojson

Resulting Shapefile

Final Shapefile

0

The primary reason for your problem is that the MultiPolygon in your source data is invalid because it has nested shells. In the OGC standand for Simple features MultiPolygon is defined as:

6.1.14 MultiPolygon

A MultiPolygon is a MultiSurface whose elements are Polygons. The assertions for MultiPolygons are as follows.

a) The interiors of 2 Polygons that are elements of a MultiPolygon may not intersect.

b) The boundaries of any 2 Polygons that are elements of a MultiPolygon may not “cross” and may touch at only a finite number of Points.

You can't define your geometry as MultiPolygon but you must wrap it inside GeometryCollection. As WKT a valid presentation would be

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

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

enter image description here

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

enter image description here

  • I do not believe that limitations of the shapefile format have any role in this play. MultiPolygon which has overlapping parts is simply an invalid geometry. – user30184 Jul 18 '17 at 10:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.