日々適当

日本の全域の地図。
国土地理院のここ [国土地理院] に置かれている地図は正距方位図法だそうです。

というわけで、その地図上に点をプロットしたい時にどうするか探ってみました。緯度経度を求めて、それをこの図法の座標にマップし直すという作業になりますね。ネットをあさった結果、やはりpyproj [code.google.com]が役に立つということが何となく分かりました。

書いてみたのが以下の物。

以下のurlを設定する場所。url = u"http://maps.google... の行が正しく書かれていないので、流用元のページを見にいったほうが吉。

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

from pyproj import Proj
import urllib
import xml.dom.minidom

#ほぼそのまま流用 http://d.hatena.ne.jp/kiwikiwi789/20100912/1284310744
ENCODING = 'utf-8'
def geocode(address):
	url = u"http://maps.google.com/maps/api/geocode/xml?&language=ja&sensor=false®ion=ja&address="
	url = url + urllib.quote(address.encode(ENCODING))
	buffer = urllib.urlopen(url).read()
	return buffer

def getLngLatFromXML( buffer ):
	dom = xml.dom.minidom.parseString(buffer)
	location = dom.getElementsByTagName('location')
	if location.length > 0:
		lat = location[0].getElementsByTagName('lat')[0].firstChild.data
		lng = location[0].getElementsByTagName('lng')[0].firstChild.data
	else:
		lat = lng = -10000.0, -10000.0
	
	return [lng, lat]
#ここまでほぼそのまま流用

addresses = {u'札幌駅':u'北海道札幌市北区北６条西4丁目', u'東京駅':u'千代田区丸の内1-9-1', u'鹿児島中央駅':u'鹿児島県鹿児島市中央町',}

jpLng0 = 139.0 + 44.0/60 + 29.0/3600
jpLat0 = 35.0 + 39.0/60 + 29.0/3600

print("zero: %f, %f"%(jpLng0, jpLat0))

convertor = Proj( proj='aeqd',  lon_0=jpLng0, lat_0=jpLat0 )

for addressKey in addresses.keys():
	_tmp = geocode( addresses[ addressKey ] )
	lngLat = getLngLatFromXML( _tmp )
	x, y = convertor( lngLat[0],  lngLat[1] )
	print(u"%s: 東経:%s, 北緯:%s, x:%f, y:%f"%( addressKey , lngLat[0],  lngLat[1] , x, y))

鹿児島中央駅、東京駅、札幌駅の位置を取得しようとするスクリプトとなります。住所から緯度経度を取得する部分はコード内のURLにある物をほぼそのまま使わせてもらっています。
実行すると、

zero: 139.741389, 35.658056
鹿児島中央駅: 東経:130.5424311, 北緯:31.5805856, x:-872753.568892, y:-412673.754915
東京駅: 東経:139.7667614, 北緯:35.6817879, x:2296.834276, y:2633.468804
札幌駅: 東経:141.3493514, 北緯:43.0688226, x:131327.426495, y:823902.392672

という値が返ってきます。xとyの値はメートルですね。日本経緯度原点から見て、どれだけの位置にあるか、という事がこれで示されます。
ざっくりと東京駅から札幌と鹿児島中央のそれぞれの地図上の位置と、帰ってきた値の位置を見比べてみると、だいたいあっているようです。
ということで、敵が正距方位図法ならなんとかなりそうな予感です。