Python >> Tutoriel Python >  >> Python Tag >> JSON

Le chevauchement du polygone devient un trou lors de la conversion de geojson en shapefile

Le format shapefile ne connaît pas la géométrie "MultiPolygone", remplacée par une géométrie Polygone.

J'utilise ici Fiona et Shapely (beaucoup plus simple que ogr, nombreux exemples dans 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)))

Il y a 2 polygones dans le 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))

Enregistrez la géométrie dans un fichier de formes

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

Ouvrez maintenant le fichier de formes

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

Par conséquent le MultiPolygon composé de deux Polygons est remplacé dans le shapefile par un Polygon composé de deux parties =un Polygon avec trou

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

Modifiez votre géométrie ou n'utilisez pas de 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}})


Prochain article