The geo_interface (GeoJSON-like) protocol was proposed by Sean Gillies and can be used in Python with:
- Shapely
- Fiona
- descartes
- geojson
- osgeo.ogr
- pygeoif
- PySAL
- ArcPy
- Papyrus
- PyShp, implemented since version 1.1.7 by Christian Lederman (https://github.com/cleder/pyshp)
- Karim Bahgat has created a tempory fork of PyShp that can write new shapes based on geoJSON: Pyshp-fork-speedup-and-geojson-write
- QGIS with the new API with monkey patching by Nathan Woodrow or with generators in GeoJSON, nouveau lingua franca en géomatique ? by Martin Laloux
- mapnik implemented by Dane Springmeyer
- rasterstats
- SpatiaLite implemented here
- PostGIS implemented here
- GeoDjango
- pygml
- geodaisy
One big advantage of the protocol is its ability to quickly examine the contents of a shapefile as dictionaries:
>>> import fiona
>>> f = fiona.open('point.shp')
>>> f.next()
{'geometry': {'type': 'Point', 'coordinates': (161821.09375, 79076.0703125)}, 'id': '0', 'properties': {u'DIP_DIR': 120, u'STRATI_TYP': 1, u'DIP': 30}}
>>> f.next()['geometry']['coordinates']
(161485.09375, 79272.34375)
>>> f.next()['properties']['DIP']
55 def records(filename):
# generator
reader = shapefile.Reader(filename)
fields = reader.fields[1:]
field_names = [field[0] for field in fields]
for sr in reader.shapeRecords():
geom = sr.shape.__geo_interface__
atr = dict(zip(field_names, sr.record))
yield dict(geometry=geom,properties=atr)
>>> import shapefile
>>> a = records('point.shp')
>>> a.next()
{'geometry': {'type': 'Point', 'coordinates': (161821.09375, 79076.0703125)}, 'properties': {'DIP_DIR': 120, 'STRATI_TYP': 1, 'DIP': 30}}
>>> a.next()['geometry']['coordinates']
(161485.09375, 79272.34375)
>>> a.next()['properties']['DIP']
55
def records(file):
# generator
reader = ogr.Open(file)
layer = reader.GetLayer(0)
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
yield json.loads(feature.ExportToJson())
>>> from osgeo import ogr
>>> a = records('point.shp')
>>> a.next()
{'geometry': {'type': 'Point', 'coordinates': (161821.09375, 79076.0703125)}, 'properties': {'DIP_DIR': 120, 'STRATI_TYP': 1, 'DIP': 30}} layer = qgis.utils.iface.activeLayer()
def records(layer):
fields = layer.pendingFields()
field_names = [field.name() for field in fields]
for elem in layer.getFeatures():
geom= elem.geometry()
atr = dict(zip(field_names, elem.attributes()))
yield dict(geometry=geom.exportToGeoJSON(),properties=atr)
c = records(layer)
c.next()
{'geometry': {'type': 'Point', 'coordinates': (161821.09375, 79076.0703125)}, 'id': '0', 'properties': {u'DIP_DIR': 120, u'STRATI_TYP': 1, u'DIP': 30}}with shapely (Sean Gillies, as with pygeoif of Christian Lederman):
>>> from shapely.geometry import shape
>>> a = records('point.shp')
>>> print shape( a.next()['geometry'])
POINT (161821.09375 79076.0703125)