【python】緯度経度から標準地域コードの取得

標準地域コードとは
wikipedia地域メッシュ
日本地図の上に碁盤の目を張り巡らせて、1つ1つの枠にユニークなコードを割り振ました。ということですね。

以降の記述は、世界測地系、10進数で緯度経度を扱うものとします。

1次メッシュ 定義&取得

• 1辺の長さが約80kmの四角形となります。
• 4桁の数字で構成され、枠の西南端の緯度を上2桁、経度を下2桁で表します。
• 緯度を1.5倍した値の上2桁 + 経度から100引いた値となります。

#encoding:utf-8
import math

class MeshCodeUtility(object):
	#1次メッシュコードの取得
	@staticmethod	def get_1st_mesh(lat, lon):
		left_operator  = int( math.floor(lat * 15 / 10 ) )
		right_operator = int( math.floor(lon - 100 ) )
		#南西端のlat,lonを取得
		dest_lat = left_operator / 15.0 * 10
		dest_lon = right_operator + 100.0
		src = {"mesh_code":str(left_operator)+str(right_operator), "lat":dest_lat, "lon":dest_lon}
		return src

1次メッシュコードの取得は簡単ですね。
同時に変換したコードが指すメッシュの西南端の緯度、経度を取得しておきます。

2次メッシュ 定義&取得

• 1次メッシュを縦、横ともに8等分した四角形(1次メッシュを64分割)。
• 1辺の長さは約4kmで、緯度5分(5/60)、経度7分30秒(7/60+30/60/60)となります。
• 1次メッシュの西南端を基点(0,0)とし、上1桁が北方向へのメッシュ位置、下1桁が東方向へのメッシュ位置を表します。
• 1次メッシュと組み合わせてnnnn-nnと記述します。

まずは、上1桁の値を取得してみましょう。

#2次メッシュコードの取得
@staticmethod
def get_2nd_mesh(lat, lon):
	#所属する1次メッシュの西南端のlat,lonを取得します。
	base_data = MeshCodeUtility.get_1st_mesh(lat, lon)
	#2次メッシュは緯度方向5分(5/60=0.08333)区切りとなります。
	left_operator = int(math.floor((lat - base_data["lat"]) * 100000 / 8333))

引数の緯度から1次メッシュの西南端の緯度を引くことによって、北方向に何度進んだ地点かがわかります。
あとは、緯度5分区切りの四角形なので、0.08333…で割りますと所属するメッシュ地点がわかります。
同様に下2桁も

#経度方向7分30秒(7/60+30/60/60=0.11666+0.008333=0.1249))区切りとなります。
right_operator = int(math.floor((lon - base_data["lon"]) * 1000 / 125))

これで、2次メッシュが取得出来ました。1次メッシュの計算時と同様に2次メッシュ西南端の緯度、経度を取得しておきます。

#南西端のlat,lon
dest_lat = left_operator * 8333.0 / 100000 + base_data["lat"]
dest_lon = right_operator * 125.0 / 1000   + base_data["lon"]
src = {"mesh_code":base_data["mesh_code"]+"-"+str(left_operator)+str(right_operator), "lat":dest_lat, "lon":dest_lon}
return src

3次メッシュ 定義&取得

• 2次メッシュを縦、横ともに10等分した四角形(2次メッシュを100分割、1次メッシュを6400分割)。
• 1辺の長さは約1kmで、緯度30秒(30/60/60)、経度45秒(45/60/60)となります。
• 2次メッシュの西南端を基点(0,0)とし、上1桁が北方向へのメッシュ位置、下1桁が東方向へのメッシュ位置を表します。
• 1次メッシュ、2次メッシュと組み合わせてnnnn-nn-nnのように記述します。

3次メッシュも分割距離が変わっただけすので、2次メッシュの計算を応用します。

#3次メッシュコードの取得
@staticmethod
def get_3rd_mesh(lat, lon):
	base_data = MeshCodeUtility.get_2nd_mesh(lat, lon)
	#3次メッシュは緯度方向 30秒(30/60/60=0.008333)区切りとなります
	left_operator = int(math.floor((lat - base_data["lat"]) * 1000000 / 8333))
	#経度方向 45秒(45/60/60=0.0125)区切りとなります。
	right_operator = int(math.floor((lon - base_data["lon"]) * 10000 / 125))
	#南西端のlat,lon
	dest_lat = left_operator * 8333 / 1000000.0 + base_data["lat"]
	dest_lon = right_operator * 125 / 10000.0 + base_data["lon"]

	src = {"mesh_code":base_data["mesh_code"]+"-"+str(left_operator)+str(right_operator), "lat":dest_lat, "lon":dest_lon}
	return src

では早速、試してみましょう。

#渋谷区役所の緯度、経度となります。
print MeshCodeUtility.get_3rd_mesh(35.6640352,139.6982122)
#{'mesh_code': '5339-35-95', 'lat': 35.658320333333329, 'lon': 139.6875}
#札幌市役所
print MeshCodeUtility.get_3rd_mesh(43.0620958,141.3543763)
#{'mesh_code': '6441-42-78', 'lat': 43.058317666666667, 'lon': 141.34999999999999}
#沖縄市役所
print MeshCodeUtility.get_3rd_mesh(26.3344266,127.8055832)
#{'mesh_code': '3927-46-04', 'lat': 26.333320000000001, 'lon': 127.8}

正しく取得出来ました。
完成したコードは、githubに置きました。
pythonのコード
phpのコード