緯度経度から住所を取得する ref: http://qiita.com/elnikkis/items/def3f325446ce7a54eff

file2.txt

$ python findcity.py
Shapefiles: ['data/japan-shape/A002005212010DDSWC35208/h22ka35208.shp', 'data/japan-shape/A002005212010DDSWC35207/h22ka35207.shp', 'data/japan-shape/A002005212010DDSWC35201/h22ka35201.shp', 'data/japan-shape/A002005212010DDSWC35203/h22ka35203.shp', 'data/japan-shape/A002005212010DDSWC35206/h22ka35206.shp', 'data/japan-shape/A002005212010DDSWC35202/h22ka35202.shp', 'data/japan-shape/A002005212010DDSWC35204/h22ka35204.shp']
35208602 山口県岩国市装束町６丁目

findcity.py

# -*- coding: utf-8 -*-

import glob
import fiona
from shapely.geometry import Polygon, shape
from revgeocoder import ReverseGeocoder


def parse_shapefile(shapefile):
    '''(area_id, Polygon, name) を返すジェネレータ'''
    with fiona.open(shapefile, 'r') as source:
        for obj in source:
            # ポリゴンデータ
            polygon = shape(obj['geometry'])

            # そのエリアの住所
            names = []
            for prop in ['KEN_NAME', 'GST_NAME', 'CSS_NAME', 'MOJI']:
                if obj['properties'][prop] is not None:
                    names.append(obj['properties'][prop])
            name = ''.join(names)

            # ユニークなidを生成（KEN + CITY + SEQ_NO2）
            area_id = int(''.join(map(str, (obj['properties']['KEN'], obj['properties']['CITY'], obj['properties']['SEQ_NO2']))))

            yield area_id, polygon, name


def main():
    shapefiles = glob.glob('data/japan-shape/A002005212010DDSWC3520*/*.shp')
    print('Shapefiles:', shapefiles)

    # インデックスを作る
    rgeocoder = ReverseGeocoder()
    id2name = {}
    def gen(shapefile):
        for area_id, polygon, name in parse_shapefile(shapefile):
            id2name[area_id] = name
            yield area_id, polygon

    for shapefile in shapefiles:
        rgeocoder.insert_from_iterator(gen(shapefile))

    # test
    area_id = rgeocoder.contains(132.257269, 34.108815)[0]
    print(area_id, id2name[area_id])

    return rgeocoder


if __name__ == '__main__':
    rgeocoder = main()

revgeocoder.py

# -*- coding: utf-8 -*-

import collections
from shapely.geometry import Polygon, Point
from rtree import index

Area = collections.namedtuple('Area', ['area_id', 'polygon'])


class ReverseGeocoder():
    def __init__(self):
        self.idx = index.Index()

    def insert_from_iterator(self, itr):
        '''(id, Polygon)を返すイテレータからRtreeを作る

        Polygon.boundをもとに作成したRtreeは、idとPolygonを保持する。

        Args:
            itr: (id, Polygon)を返すイテレータ
        '''
        for i, (area_id, polygon) in enumerate(itr):
            obj = Area(area_id=area_id, polygon=polygon)
            self.idx.insert(i, polygon.bounds, obj)

    def contains(self, lat, lon):
        '''Point(lat, lon)を含むPolygonのarea_idを返す。

        2つ以上のPolygonとマッチした場合は面積が小さい順にソートして返す。
        （包含関係になっているエリアがあるとき面積が小さいほうを選ぶとよい。
        しかし実際は、PolygonとPolygonのあいだが重なっていたり……。）
        '''
        result = []
        point = Point(lat, lon)
        for hit in self.idx.intersection(point.bounds, objects=True):
            if hit.object.polygon.contains(point):
                result.append(hit.object)
        if len(result) > 1:
            result.sort(key=lambda x: (x.polygon.area, x.area_id))
        return [r.area_id for r in result]

    def __repr__(self):
        return '<ReverseGeocoder contains {} polygons>'.format(self.idx.count(self.idx.bounds))