reverse_geo.py

#!/usr/bin/env python
# -*- coding:utf-8 -*-

import sys
import sqlite3
import codecs
import Geohash
import pyproj

from collections import defaultdict

geodict = defaultdict(list)

with codecs.open(sys.argv[1], 'r', 'shift-jis') as f:
  f.next()
  for line in f:
    line = line[:-1].split(',')

    lat = float(line[6].strip('"'))
    lon = float(line[7].strip('"'))

    data = {
      'pref_code': line[1].strip('"'),
      'city_code': line[3].strip('"'),
      'block_code': line[5].strip('"'),
      'lat': lat,
      'lon': lon
    }

    geodict[Geohash.encode(lat, lon, precision=3)].append(data)

from pyproj import Geod
geod = Geod(ellps='WGS84')
def distance(p1, p2):
  fa, ba, d = geod.inv(p1[0], p1[1], p2[0], p2[1])        
  return d

lat = float(sys.argv[2])
lon = float(sys.argv[3])

geohash = Geohash.encode(lat, lon, precision=3)

if geohash in geodict:
  data = geodict[geohash]
  if len(data) == 1: 
    print data[0]['pref_code'], data[0]['city_code'], data[0]['block_code']
  else:
    best = None
    best_d = sys.maxint
    for item in data:
      d = distance((lon, lat), (item['lon'], item['lat']))
      if d < best_d:
        best = item
        best_d = d
    print best['pref_code'], best['city_code'], best['block_code']
else:
  print 'Not found'