8分の1地域メッシュを緯度経度から計算する

Share Button

SASじゃなくともRなりExcelなり何かしらの言語では絶対どっかに書いてあるだろうからそれをリファクタリングしようと思ったのにマジでどこにも書いてないし仕方ないから自分で作るわ。
やり方自体はパッケージ使うわけでもないしSASだろうとPythonだろうとExcelだろうと関数名が変わるぐらいなんだけど。

別段難しくはないんだけど、位置情報界隈はとりあえず考え方がめんどくさい。
とりあえず簡単に理解し直すために各メッシュを図解する。

まず1-3次メッシュのサイズ感は


80km四方の1次メッシュを縦横8等分したのが10km四方の2次メッシュ、それをさらに10等分したのが1km四方の3次メッシュとなるわけだ。もう少し見やすくしたかったけどサイズ的にこれが限界だった。まあ1次と3次で80:1ではしかたないな。

んで3次メッシュを基準に1/2、1/4、1/8メッシュが


二分の一メッシュ以下の細かいメッシュ(専門用語がわからないのでここでは分数メッシュと呼んでおく)の算出に関する基本的な考え方は、とにかく基準となるひとつ大きなサイズのメッシュ(親メッシュと呼んでおく)の1枠を縦横2等分して、そのうちのどこの象限に属するかを左下、右下、左上、右上の順でナンバリングする。だから二分の一メッシュなら親メッシュたる3次メッシュを4分割(画像でいう深緑の枠と番号)するわけだ。そもそもこの番号付の順はなんなんだろう。

SASでの3次メッシュの算出まではググれば有能な後輩ちゃんのプログラムが落ちてるのでそれを読んでいただくとして、その先の分数メッシュの算出に関しては所属する親メッシュの中心点さえわかればそこと自分とこの位置を比較することで分類できる。つまり親メッシュ算出の過程で出てきた剰余の大小を見ればいい。

たとえば日本橋駅の緯度経度は
緯度:35度40分56秒
経度:139度46分29秒
Googleマップで緯度・経度を求めるより)

度分秒表記がまず死ぬほど鬱陶しいので十進数に直す。分は60で、秒は60*60=3600で割ればいいので、
緯度lat
= 35 + 40/60 + 56/3600
= 35.682222222
経度long
= 139 + 46/60 + 29/3600
= 139.774722222
検算と位置の確認は『緯度・経度の、10進数と60進数(度分秒)の変換』でやればいいです。Google Maps上でもプロットされるからわかりやすいですね。まあ一旦度分秒表記に戻してからプロットされるから数メートルぐらいズレるけど。
そしてこの時点での10進数での単位は度になる。ここに60かければ全体が分になるし,3600かければ秒になる。実際他のページ見てると全部最初に十進法に変換しといてから単位の次元を落としつつ計算を進めてますね。

1次メッシュに変換

一応自分でやっとかないとミスっても気づかないので念のためここからやっておく。
計算方法の詳細は他に譲るとして、要は1次メッシュの4桁の数字は上二桁が緯度の1.5倍切り捨て、下二桁が経度の下二桁切り捨てになるので、順々に計算して
mesh1_1
= int(lat * 60 / 40)
= int(35.682222222*1.5)
= 53

mesh1_2
= int(long – 100)
= int(139.774722222 – 100)
= 39

mesh1
= mesh1_1 * 100 + mesh1_2
= 53 * 100 + 39
= 5339

度分秒表記もそうなんだけど、単位や意味合いの異なる数字の連結にまずイライラするよね。
確認すると『地図上で標準地域メッシュを確認するページ』で見る限り23区内は大体5339で合ってますね。

2次メッシュに変換

さっきの計算の剰余が必要なので、
緯度部分 lat_r1
= lat * 60 – mesh1_1 * 40
= 35.682222222 * 60 – 53 * 40
= 20.93333332(分)

経度部分 long_r1
= long – 100 – mesh1_2
= 139.774722222 – 100 – 39
= 0.774722222(度)

これら剰余からまた2次メッシュそれぞれの数字を順に計算して、
mesh2_1
= int(lat_r1 / 5)
= int(20.93333332 / 5)
= 4

mesh2_2
= int(long_r1 * 60 / 7.5)
= int(0.774722222 * 60 / 7.5)
= 6(分)

mesh2
= mesh2_1 * 10 + mesh2_2
= 4 * 10 + 6
= 46(分)

mesh1とmesh2で5339-46、2次メッシュコード533946の位置を見る限り正しいですね。

3次メッシュに変換

再び剰余がいるので
緯度部分 lat_r2
= lat_r1 – mesh2_1 * 5
= 20.93333332 – 4 * 5
= 0.93333332(分)

経度部分 long_r2
= long_r1 * 60 – mesh2_2 * 7.5
= 0.774722222 * 60 – 6 * 7.5
= 1.48333332(分)

mesh3_1
= int(lat_r2 * 60 / 30)
= int(0.93333332 * 60 / 30)
= 1(秒)

mesh3_2
= int(long_r2 * 60 / 45)
= 1.48333332 * 60 / 45
= 1(秒)

mesh3
= mesh3_1 * 10 + mesh3_2
= 11
よって3次メッシュは53394611ですね。同じく53394611を確認すれば…あってますね。

さて、ここまでの計算の流れが理解できていればその先は難しいことはないのです。ここまでどんな計算をしてきたかというと、結局60かけて度→分、分→秒に変換しつつ、たとえばmesh2_1の計算次には5で割っていて、mesh2_2では7.5で割っているけど、これはWikipediaに言う所の2次メッシュの「緯度差は5分、経度差は7分30秒」に沿っているだけ。つまり親メッシュの算出時に出た剰余の単位を合わせた上で、子メッシュのメッシュの大きさに分割しているに過ぎない。ということはその先は半分に分割しているだけなのでもっと簡単ですわな。

以下に簡単にメッシュのサイズを貼っておくと、
メッシュ種別 緯度幅 経度幅
1次メッシュ 40分 1度
2次メッシュ 5分 7分30秒
3次メッシュ 30秒 45秒
1/2メッシュ 15秒 22.5秒
1/4メッシュ 7.5秒 11.25秒
1/8メッシュ 3.75秒 5.625秒
縮尺に従って割るだけです。

ただまあ、分数メッシュの面倒なところは、単なる大小関係では位置が決まらないことですね。緯度経度それぞれの位置に基づいてナンバリングしないといけないのが最高にとろくさいです。

1/2メッシュに変換

ということで、まずは3次メッシュ算出の剰余の計算から。これ以降単位は全部秒にします。
lat_r3
= lat_r2 * 60 – mesh3_1 * 30
= 0.93333332 * 60 – 1 * 30
= 25.9999992(秒)

単純な話、親メッシュの緯度差が30秒なので剰余が30超えてたらそれはもうなんかおかしいわけです。同じく経度も、

long_r3
= long_r2 * 60 – mesh3_2 * 45
= 1.48333332 * 60 – 1 * 45
= 43.9999992(秒)

こっちも幅は45秒なんで大丈夫ですね。
そして次にどの象限に入るのか確認します。

1/2メッシュは緯度差が15秒なので、
1_2mesh_1
= int(lat_r3 / 15)
= int(25.9999992 / 15)
= 1
2等分する以上は1_2mesh_1は0から2まで(2は含まない。2を超える[子メッシュの幅2つ分を超える]ということは親メッシュでひとつ隣の枠に行ってしまう)の値を切り下げるので、簡単にいえば0, 1の離散値になります。緯度は値が大きい方が北なので、親メッシュたる3次メッシュを南北で2等分した上で、0なら南半分、1なら北半分になるわけです。

同じく経度差は22.5秒なので、
1_2mesh_2
= int(long_r3 / 22.5)
= int(43.9999992 / 22.5)
= 1
こちらも0/1の離散値になります。経度は値の大きい方が東になるので、親メッシュたる3次メッシュを東西で2等分した上で、0なら西半分,1なら東半分ですね。

1_2meshは1,2の結果として北東になりました。
同じく3次メッシュ53394611を下の画像のように4分割すると、確かに日本橋駅は北東です。


従って1/2メッシュの値は4ですね。ここからはナンバリングが必要になります。
{1_2mesh_1, 1_2mesh_2} → 1_2meshとするとその関係は
{0, 0} → 1
{0, 1} → 2
{1, 0} → 3
{1, 1} → 4
となっています。ということでその関係を定式化するとすれば
1_2mesh = 1 + 2 * 1_2mesh_1 + 1_2mesh_2
ですね。緯度メッシュの増加で1/2メッシュの値は2増え、経度メッシュに対しては1増えています。あとはどちらも0の場合に切片が1だと考えれば正しいですね。

1/4メッシュに変換

ここからは基本的なやり方は同じですね。ひたすら剰余を出しては子メッシュのスケールで割るだけです。単位も秒から下には下がりません。
lat_r4
= lat_r3 – 1_2mesh_1 * 15
= 25.9999992 – 1 * 15
= 10.9999992(秒)

long_r4
= long_r3 – 2mesh_2 * 22.5
= 43.9999992 – 1 * 22.5
= 21.4999992(秒)

1_4mesh_1
= int(lat_r4 / 7.5)
= int(10.9999992 / 7.5)
= 1

1_4mesh_2
= int(long_r4 / 11.25)
= int(21.4999992 / 11.25)
= 1

定式化した式は変わらないので
1_4mesh
= 1 + 2 * 1_4mesh_1 + 1_4mesh_2
= 1 + 2 + 1
= 4
結果1_4mesh = 4でした。画像で確認しても


確かに北東4ですね。

1/8メッシュに変換

最後です。サクサクやっていきます。
lat_r5
= lat_r4 – 1_4mesh_1 * 7.5
= 10.9999992 – 1 * 7.5
= 3.4999992(秒)

long_r5
= long_r4 – 1_4mesh_2 * 11.25
= 21.4999992 – 1 * 11.25
= 10.2499992(秒)

1_8mesh_1
= int(lat_r5 / 3.75)
= int(3.4999992 / 3.75)
= 0

1_8mesh_2
= int(long_r5 / 5.625)
= int(10.2499992 / 5.625)
= 1

1_8mesh
= 1 + 2 * 1_8mesh_1 + 1_8mesh_2
= 1 + 2 * 0 + 1
= 2
ですね。


ずっと4ばっかりだったので場所選びミスったなと思ってたんですが、2で正しいですね。

ということでこれをまとめます。(このまま順に走らせてももちろんいける)
lat_deg = lat1 + lat2/60 + lat3/3600
long_deg = long1 + long2/60 + long3/3600 
mesh1_1 = int(lat_deg * 60 / 40)
mesh1_2 = int(long_deg - 100)
mesh1 = mesh1_1 * 100 + mesh1_2
_lat_r1 = lat_deg * 60 - mesh1_1 * 40
_long_r1 = long_deg - 100 - mesh1_2
mesh2_1 = int(_lat_r1 / 5)
mesh2_2 = int(_long_r1 * 60 / 7.5)
mesh2 = mesh2_1 * 10 + mesh2_2
_lat_r2 = _lat_r1 - mesh2_1 * 5
_long_r2 = _long_r1 * 60 - mesh2_2 * 7.5
mesh3_1 = int(_lat_r2 * 60 / 30)
mesh3_2 = int(_long_r2 * 60 / 45)
mesh3 = mesh3_1 * 10 + mesh3_2
_lat_r3 = _lat_r2 * 60 - mesh3_1 * 30
_long_r3 = _long_r2 * 60 - mesh3_2 * 45
1_2mesh_1 = int(_lat_r3 / 15)
1_2mesh_2 = int(_long_r3 / 22.5)
_lat_r4 = _lat_r3 - 1_2mesh_1 * 15
_long_r4 = _long_r3 - 2mesh_2 * 22.5
1_4mesh_1 = int(_lat_r4 / 7.5)
1_4mesh_2 = int(_long_r4 / 11.25)
1_4mesh = 1 + 2 * 1_4mesh_1 + 1_4mesh_2
_lat_r5 = _lat_r4 - 1_4mesh_1 * 7.5
_long_r5 = _long_r4 - 1_4mesh_2 * 11.25
1_8mesh_1 = int(_lat_r5 / 3.75)
1_8mesh_2 = int(_long_r5 / 5.625)
1_8mesh = 1 + 2 * 1_8mesh_1 + 1_8mesh_2
よく見たら変数名が競合してますね。まあ上から順にやって行く分には問題はないんだけど。
投入時のlat1/long1, lat2/long2, lat3/long3にはそれぞれ度、分,秒をいれてもらえれば、10進数(単位:度)に直した上で処理を走らせます。初めから度数のみになっているのであれば頭2行の処理は不要ですね。
途中の剰余の変数にはすべて頭にアンダーバーがつけてあるので、SASならdrop _:;とか、python+pandasならdf[[v for v in df.columns if not v.startswith(‘_’)]]とかで一括で落としてもらえれば消えます。Rは知りません。

style="display:block; text-align:center;"
data-ad-layout="in-article"
data-ad-format="fluid"
data-ad-client="ca-pub-3546003055292762"
data-ad-slot="5749192034">


Share Button

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

This site uses Akismet to reduce spam. Learn how your comment data is processed.