Pythonで緯度経度から地域メッシュコードを算出してみる

どうも、
前回はこんな記事を書いてました。

pythonで電波伝播モデルから発信機の位置を計算してみる[Wi-Fi, Beacon編]

なので、順当にいけばGPS編だったんですけど、昨今のコロナ事情で弊社も辛く、
時間がカツカツなので別テーマで一旦場を持たせようと思います。

日本の地域メッシュの概要について

概要や全体感を知りたいだけなら、
弊社メンバーの別記事で地域メッシュについて概要を書いているのでそちらをお読みください。

日本の地域メッシュを地図上に可視化して理解する

こちらはもう少し「規格」とかディープな方の話になります。
数式もちょっと出ます。

メッシュコードの要点

いくつか見た中ではこの記事が一番わかりやすかったです。
標準地域メッシュ・システム

上記の内容と合わせてWikipediaやJIS X 0410:2002を読むと、
だいたい以下のことが基本的な知識になります。

・メッシュは四角形
・日本の国土を埋めるような形で地域メッシュが網の目に張り巡らされている(海上は除外)
・基準地域メッシュと呼ばれる1km四方のメッシュが存在する
・基準地域メッシュは通称1次メッシュ(80km四方)と呼ばれるものの一辺を80分の1分割したもの
・基準地域メッシュを分割したり、組み合わせたものを「分割/統合」地域メッシュとして区別する
・各地域メッシュの番号/桁は分割の仕方と密接に関係している

以上です。

要は
いろいろな名前のメッシュがあるが、全ては1次メッシュの派生メッシュ
約80km四方の1次メッシュを分割したり、統合したりして他のメッシュは作れる
分割/統合の仕組みを理解していれば、メッシュコードからメッシュの位置(緯度経度)が計算できる
と言うことです。

例として書くと
1次メッシュ(四桁)の一辺を8分割 = 2次メッシュ(六桁)
2次メッシュ(六桁)の一辺を10分割 = 基準地域メッシュ(八桁)
基準地域メッシュ(八桁)の一辺を2分割 = 4次メッシュコード(九桁)
・・・

となっているため、基本的には簡単な四則演算(と除算した余りを取り出す計算)で全てが再現できます。
文字だけではイメージしづらいと思うので、
こちらの表4から1次メッシュの採番ルールと、1次メッシュを2次メッシュにする絵を引用します。

緯度経度→メッシュコード

以上から、

①メッシュの採番を自分で再現するには南西端の緯度経度が必要
②メッシュの番号がわかれば南西端の緯度経度は逆算できる

ことがわかります。

このエントリではまず①を1次〜4次メッシュまで行います。
いつか余裕があれば②を1次〜4次メッシュまで行います。
言語はPythonです。

虎ノ門ヒルズの1次〜4次メッシュコードを出してみる

弊社のある虎ノ門ヒルズの緯度経度と3次メッシュ

緯度 / 経度 = 35.666863 / 139.74954 # 虎ノ門ヒルズ

3次メッシュ上の位置関係

このサイトで緯度経度から3次メッシュまでの変換ができるので、先に検索したものを使ってみてみます。

上記の画像で「メッシュ番号53394509」の吹き出しが出て居るあたりが虎ノ門ヒルズです。
メッシュ採番のルールから

3次メッシュコード　= 5339 45 09
2次メッシュコード　= 5339 45
1次メッシュコード　= 5339

となることがわかります。
最後に、4次メッシュコードを調べます。
画像の位置関係から虎ノ門ヒルズは3次メッシュコードの東南隅にある事がわかります。
4次メッシュコードは3次メッシュコード1辺を2分割し、左下/右下/左上/右上と採番しているので、

3次メッシュコード(全体)＝5339 45 09
4次メッシュコード(南西)＝5339 45 09 1
4次メッシュコード(南東)＝5339 45 09 2 ⇦ これ
4次メッシュコード(北西)＝5339 45 09 3
4次メッシュコード(北東)＝5339 45 09 4

となります。念の為、分割メッシュコードの採番ルールを説明してるサイトを再掲しておきます。
標準地域メッシュ・システム

以上で前段は終わりです。それではコードを動かしてみましょう。

コード

# -*- coding: utf-8 -*-
import math as mt

class latlon2grid:
    def grid1st(lon, lat):  # 1次メッシュ(4桁) 分割なし
        return int(mt.floor(lat*1.5)) * 100 + int(mt.floor(lon-100))

    def grid2nd(lon, lat):  # 2次メッシュ(6桁) 8分割
       return (int(mt.floor(lat*12       / 8))   * 10000 + int(mt.floor((lon-100)*8         / 8))  * 100   +   
               int(mt.floor(lat*12 %  8     ))   * 10    + int(mt.floor((lon-100)*8))  %  8               )  

    def grid3rd(lon, lat):  # 3次メッシュ(8桁) 8分割x10分割=80分割
        return (int(mt.floor(lat*120      / 80)) * 1000000 + int(mt.floor((lon-100))             ) * 10000 +  # 1次メッシュ
                int(mt.floor(lat*120 % 80 / 10)) * 1000    + int(mt.floor((lon-100)*80 % 80 / 10)) * 100 +    # 2次メッシュ
                int(mt.floor(lat*120 % 10))      * 10      + int(mt.floor((lon-100)*80)) % 10               ) 

    def grid4th(lon, lat):  # 4次メッシュ(9桁) 8分割x10分割x2分割=160分割
        return (int(mt.floor(lat*240       / 160)) * 10000000 + int(mt.floor((lon-100)*160       / 160)) * 100000 +    # 1次メッシュ
                int(mt.floor(lat*240 % 160 / 20))  * 10000    + int(mt.floor((lon-100)*160 % 160 / 20))  * 1000   +    # 2次メッシュ
                int(mt.floor(lat*240 % 20  / 2))   * 100      + int(mt.floor((lon-100)*160 % 20  / 2))   * 10     +    # 3次メッシュ
                int(mt.floor(lat*240)) % 2         * 2        + int(mt.floor((lon-100)*160)) % 2                  + 1) # 4次メッシュ

if __name__ == "__main__":
    lon = 139.74954 # 虎ノ門ヒルズ
    lat = 35.666863 # 虎ノ門ヒルズ
    print('1次メッシュ： ', latlon2grid.grid1st(lon, lat))  # 5339
    print('2次メッシュ： ', latlon2grid.grid2nd(lon, lat))  # 533945
    print('基準地域メッシュ： ', latlon2grid.grid3rd(lon, lat))  # 53394509
    print('4次メッシュ： ', latlon2grid.grid4th(lon, lat))  # 533945092

どうでしょうか？一定のパターンが見えるかと思います。
複雑な数式に見えますが、やって居ることは単純で
南西の緯度経度をメッシュの次数ごとにメッシュコードに変換。その後桁数を揃えているだけです。
メッシュコードはIDなので、整数データではなく文字データだと思うのですが、
パターンが定義されていると数式で導きだせる事がわかります。

上記コードの出力結果は以下になります。

1次メッシュ： 5339
2次メッシュ： 533945
基準地域メッシュ： 53394509
4次メッシュ： 533945092

問題なさそうですね。
同じ要領で5次メッシュ以降もできるので、興味のある方はご自身でお試しください。

余談

日本の地域メッシュを地図上に可視化して理解する
にあるように、日本の地域メッシュは海の上では定義されていません。

しかしながら、メッシュコードの採番ルール自体は緯度経度からなされるため、
海の上においても番号自体は振られます。

じゃあ、これで世界中を覆えるんじゃないか？という話ですと・・・

赤道方向の直径は、約1万2756kmで、北極と南極の方向の直径は、約1万2714km。

緯度は 12756km / 80km = 159.45
経度は 12756km / 80km = 158.925

と言った感じで帯に短し襷に長しでどうにも割り切れない結果になりました。
概ね覆う事はできるけども・・・緯度経度のどちらかが0になる付近特に本初子午線と赤道がぶつかるアフリカ西の沖？あたりが重複してしまいそうです。 まぁ、海の上なんでぶっちゃけ大した問題じゃありませんし、本当は厳密に80km四方になっていない感じなので思い過ごしかもしれません

実際問題として、地上は球面なので球面幾何学的に考える事は大局的には必要になってきます。2点間の距離は局地的にはユークリッドの距離で近似できますが、2点が離れれば離れるほど近似誤差は大きくなります。このへんについては常に注意が必要で、個人的に深掘りしてみたい分野です。

四角形で埋めるのがいいのか？
三角形か、六角形か？
はたまたサッカーボールの様に五角形と六角形を組み合わせるべきか？

もしみなさんが位置情報サービスを検討される際には、
どの様なメッシュの形/サイズが理想的なのか、考えてみてはいかがでしょうか？