GeoJSON and the geo interface for Python

Sean Gillies

2013-06-27 00:00

It's pretty clear that GeoJSON is good for the web. I think it's also good for Python programmers in two different ways.

Fiona employs Python mappings modeled on GeoJSON to represent GIS features and their geometries. By reusing these existing concepts – the GeoJSON data model and the Python mapping type – Fiona keeps its mental footprint small. There's less to memorize, less need to check documentation, fewer gotchas. There's also less code and less room for bugs. For example, given a Fiona feature, there is no need to look up the classes and methods for accessing the feature's non-geometry fields, it's just feature['properties'].items().

>>> import fiona
>>> with fiona.open(
...         "/Users/seang/data/ne_50m_admin_0_countries/"
...         "ne_50m_admin_0_countries.shp") as c:
...     first = next(c)
...
>>> items = sorted(first['properties'].items())
>>> import pprint
>>> pprint.pprint(items)
[('abbrev', 'Aruba'),
 ('abbrev_len', 5.0),
 ('adm0_a3', 'ABW'),
 ('adm0_a3_is', 'ABW'),
 ('adm0_a3_un', -99.0),
 ('adm0_a3_us', 'ABW'),
 ('adm0_a3_wb', -99.0),
 ('adm0_dif', 1.0),
 ('admin', 'Aruba'),
 ('brk_a3', 'ABW'),
 ('brk_diff', 0.0),
 ('brk_group', None),
 ('brk_name', 'Aruba'),
 ('continent', 'North America'),
 ('economy', '6. Developing region'),
 ('featurecla', 'Admin-0 country'),
 ('fips_10', None),
 ('formal_en', 'Aruba'),
 ('formal_fr', None),
 ('gdp_md_est', 2258.0),
 ('gdp_year', -99.0),
 ('geou_dif', 0.0),
 ('geounit', 'Aruba'),
 ('gu_a3', 'ABW'),
 ('homepart', -99.0),
 ('income_grp', '2. High income: nonOECD'),
 ('iso_a2', 'AW'),
 ('iso_a3', 'ABW'),
 ('iso_n3', '533'),
 ('labelrank', 5.0),
 ('lastcensus', 2010.0),
 ('level', 2.0),
 ('long_len', 5.0),
 ('mapcolor13', 9.0),
 ('mapcolor7', 4.0),
 ('mapcolor8', 2.0),
 ('mapcolor9', 2.0),
 ('name', 'Aruba'),
 ('name_alt', None),
 ('name_len', 5.0),
 ('name_long', 'Aruba'),
 ('name_sort', 'Aruba'),
 ('note_adm0', 'Neth.'),
 ('note_brk', None),
 ('pop_est', 103065.0),
 ('pop_year', -99.0),
 ('postal', 'AW'),
 ('region_un', 'Americas'),
 ('region_wb', 'Latin America & Caribbean'),
 ('scalerank', 3),
 ('sov_a3', 'NL1'),
 ('sovereignt', 'Netherlands'),
 ('su_a3', 'ABW'),
 ('su_dif', 0.0),
 ('subregion', 'Caribbean'),
 ('subunit', 'Aruba'),
 ('tiny', 4.0),
 ('type', 'Country'),
 ('un_a3', '533'),
 ('wb_a2', 'AW'),
 ('wb_a3', 'ABW'),
 ('wikipedia', -99.0),
 ('woe_id', -99.0)]

If we parse the GeoJSON version of that shapefile using Python's own standard json library, we get the same data.

>>> import requests
>>> r = requests.get(
...     "https://raw.github.com/sgillies/sgillies.github.com/master/"
...     "ne-sample/ne_50m_admin_0_countries.json")
>>>
>>> r.ok
True
>>> r.json()['type'] == 'FeatureCollection'
True
>>> first = r.json()['features'][0]
>>> sorted(first['properties'].items()) == items
True

The same data, accessed in the very same way. A familiar uniform interface for both Shapefile and JSON features.

For a while I've been suggesting that GeoJSON, or rather Python mappings modeled on GeoJSON, might be as useful a protocol for passing geospatial data between Python packages as it is for passing data around the web. A few Python packages have implemented the Python Geo Interface: ArcPy, descartes, Fiona, PySAL, Shapely. Descartes can, for example, help you draw ArcPy, Fiona, PySAL, or Shapely objects; or any other object that satisfies the protocol. It's dirt-simple Python GIS interoperability based on mappings with agreed-upon items.

Yesterday, I saw Nathan Woodrow's post about monkey patching the protocol into QGIS classes. I'm happy to see this because adoption of the protocol by QGIS will be a big deal and also because it demonstrates a manner in which Esri Python users might fix ArcPy's slightly broken implementation of the protocol. In the same discussion, Calvin Metcalf confirmed for me that ArcPy feature objects also satisfy the protocol, which means that writing GeoJSON from ArcGIS using Python could, in theory, be very simple.

import json
json.dumps(
    dict(
        type='FeatureCollection',
        features=[row.__geo_interface__ for row in cursor]))

Esri's support for GeoJSON is (potentially) better than most people realize.

Comments

Re: GeoJSON and the geo interface for Python

Author: Nathan

I found a easier, and better way to attach the __geo_interface__ to QgsFeature.

QgsFeature.__geo_interface__ = property(mapping_feature)

Much nicer then the MethodType stuff and fits in nicer with what you have with Fiona.

Re: GeoJSON and the geo interface for Python

Author: Martin

The geo_interface was also implemented in Pyshp by Christian Lederman (https://github.com/cleder/pyshp)

import shapefile

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)

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

and you can use the same method to implement it in QGIS 2

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)

>>> layer = qgis.utils.iface.activeLayer()

>>> c = records(layer)

>>> c.next()

{'geometry': u'{ "type": "LineString", "coordinates": [ [164917.66073716, 115230.16694565], [170806.15565476, 116401.17445767] ] }', 'properties': {u'id': 1,u'test':u'ligne'}}

My problem is to implement it in writing as in Fiona