緯度経度から住所を取得する

公開されているポリゴンデータを使って、ある座標を含むポリゴンを探すプログラムを公開します。緯度経度から住所を取得することはReverse geocodingと呼ばれています。

ポリゴンデータを持っていない場合は、Google Maps APIなどでもReverse geocodingが提供されていて、簡単に使えます。ポリゴンデータを持っていてReverse geocodingをしたい場合は、MySQLなどのデータベースのGIS機能を使うこともできます。

今回は、大量の座標に対して、日本の中のどの市区町村に属するものかを判定したいときのPythonプログラムです。毎回インデックス生成をしているような雑魚プログラムでも、座標が多いのでMySQL 5.6を使う場合より実行にかかる時間は短かったです。

MySQLを使うメリット：

• データさえあればデータを入れてSQLを書くだけ

デメリット：

• ほしい住所レベルに合わせたポリゴンデータが必要（番地レベルのポリゴンはなかなか手に入らないかも）

Pythonで実装した理由：

• 番地レベルまで細かくなくていい
• データベースにアクセスしてると量が多くて終わらない（もちろんAPIもつらい）
• 照合に使ったエリアデータの出どころが明らか

必要なライブラリ

• Rtree
• Shapely

依存ライブラリがあるのでpip install rtree shapelyだけじゃうまくいかないかもしれません。なのでインストールはそれぞれのドキュメントを読んでください。

緯度経度からの住所の照合

rtreeとshapelyを使うとこう書けます。

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

使い方

e-statからダウンロードできる境界データを使う場合はこんなかんじになります。

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

実行例：

$ 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 山口県岩国市装束町６丁目

その他類似ライブラリとの比較

• reverse_geocoder https://github.com/thampiman/reverse-geocoder
  • エリアデータはGeoNamesのものを使っている
  • 距離的に近い街を返しているっぽい？（やっていることが違う）
• reverse_geocode https://bitbucket.org/richardpenman/reverse_geocode
  • たぶん上と同じ