逆ジオコーディングをpythonで

逆ジオコーディング

大量に座標から住所変換したいのだが
APIは無料範囲では使い物にならんし、こんなことに金払う気にもならない。ひとつあったのが、Javascriptだったりして動かし方がわからない。
他にスタンドアロンで動くものがないから作った。

検索用のデータはGithubにおいた

import pickle
import os
import re
from math import cos, sin, sqrt, atan2
import json
import sys

datadir = os.path.join(os.path.dirname(__file__), "dat")
maxint = sys.maxsize

"""
    e-Stat 統計地理情報システム 境界データダウンロードより取得したデータに基づいて作成
    https://www.e-stat.go.jp/gis/statmap-search?page=1&type=2&aggregateUnitForBoundary=A&toukeiCode=00200521&toukeiYear=2020&serveyId=A002005212020&datum=2011
    この作品は、クリエイティブ・コモンズの 表示 4.0 国際 ライセンスで提供されています。ライセンスの写しをご覧になるには、 http://creativecommons.org/licenses/by/4.0/ をご覧頂くか、Creative Commons, PO Box 1866, Mountain View, CA 94042, USA までお手紙をお送りください。 
"""

def get_address(lat, lon):
    def contains(minmax):
        max_x, max_y, min_x, min_y = minmax
        if max_x < lat:
            return False
        if min_x > lat:
            return False
        if max_y < lon:
            return False
        if min_y > lon:
            return False
        return True

    def distance(dic):
        dlat = dic["X_CODE"] - lat
        dlon = dic["Y_CODE"] - lon

        if dlat ** 2 > 1 or dlon ** 2 > 1:
            return maxint

        a = sin(dlat / 2) ** 2 + cos(dic["X_CODE"]) * cos(lat) * sin(dlon / 2) ** 2
        return 12742.02 * atan2(sqrt(a), sqrt(1 - a))

    dist = maxint
    ret = {}

    for i in [i for i, (minmax, _) in enumerate(areas) if contains(minmax)]:
        for j in [j for j, (minmax, _) in enumerate(areas[i][1]) if contains(minmax)]:
            for dic in areas[i][1][j][1]:
                d = distance(dic)
                if d < dist:
                    dist = d
                    ret = dic

    return ret


if __name__ == "__main__":
    areas = pickle.load(open(datadir + "/revgeo.pkl", "rb"))

    for x in sys.argv[1:]:
        print(get_address(*x.split(",")))

但し書き

辞書も軽く、軽快に検索できるが精度を求める場合はhttps://y-mattu.hatenablog.com/entry/2017/09/18/185014
などを参考にした方が良い。

なぜかというと、上記のプログラムはポリゴンではなく町字レベルの住所毎に国土交通省が設定した重心点の座標があるので、直線距離が最も近い町字レベルの住所を表示しているのに過ぎない。
川を隔てて隣り町の中心の方が距離が近かったらそっちの住所になったりすることもあるので、今回の僕の場合みたいに多少違っててもいいから高速に住所を取りたい場合だけ使用できる。

今後のために何か役に立つかも？

def isin_poligon(x, y, polygon):
    """
    ポリゴンは反時計回りであること！など制約事項があるため
    今回使わなかったが、ポリゴンの内側にあるかどうかを判定するアルゴリズムもあった。
    今後何かの役に立つかもしれない。
    """
    cnt = 0
    L = len(polygon)
    for i in range(L):
        x0, y0 = polygon[i - 1]
        x1, y1 = polygon[i]
        x0 -= x
        y0 -= y
        x1 -= x
        y1 -= y

        cv = x0 * x1 + y0 * y1
        sv = x0 * y1 - x1 * y0
        if sv == 0 and cv <= 0:
            # a point is on a segment
            return True

        if not y0 < y1:
            x0, x1 = x1, x0
            y0, y1 = y1, y0

        if y0 <= 0 < y1 and x0 * (y1 - y0) > y0 * (x1 - x0):
            cnt += 1
    if (cnt % 2 == 1):
        return True