日本国土の正方メッシュ座標

はじめに

無線通信システムの国土カバー率を求めたい、日本の地域に依存した電波伝搬特性を調べたい、などの理由で、日本の陸地全域について、一定間隔距離の座標列を知りたくなることがあります。

地球中心から赤道からの南北方向の角度である緯度については、地球1周は約4万キロメートルなので、1度あたりの距離は約111キロメートルです。一方、イギリスのグリニッジ天文台からの東西方向の角度である経度については、1度あたりの距離は緯度に依存し、先の約111キロメートルに緯度の余弦（コサイン）をかけた値になります。

そこで、大きくて座標のみやすい地図、例えば区分航空路に一定の緯度と経度ごとの線を引き、陸地の座標範囲をすべて読み取ろうとしましたが、すぐにこれは大変な手間のかかることだと気づきました。

代案を考えていたところ、地域ごとの統計を得るのには、総務省統計局の地域メッシュコードを利用すれば良いことに気づきました。

地域メッシュコード

地域メッシュは、緯度と経度に基づいて、地域をほぼ同じ大きさの網の目にわけたものです。総務省統計局で公開されているメッシュコードは、47都道府県（ただし北海道は3つに分かれています）地域別にCSV（comma separated values）形式にて公開されています。

メッシュコードは、緯度経度を数値化し階層化され、一辺のサイズに応じた次数が定義されます。このメッシュコードの大きさは、次数で表されます。1次メッシュは一辺80キロメートルで4桁表示、2次メッシュは一辺10キロメートルで6桁表示、3次メッシュは一辺1キロメートルで8桁表示です。そして、3次メッシュコードの末尾2桁を除くと、その2次メッシュコードになるように、階層化されています。統計局の地域メッシュコードは、「都道府県市町村コード」（5桁）とともに、市町村名と3次メッシュコードが記載されています。このメッシュコードには、日本人であっても訪れることが極めて困難な地域のものも含まれています。

はじめに、この地域メッシュCSVファイルをすべてダウンロードします。地域メッシュコードから座標への変換公式は地域メッシュ統計の特質・沿革（PDF：1,249KB）の13ページに書かれていますが、ここではGIS奮闘記・Python で地域メッシュコードを緯度経度に変換する方法のコードを関数化して利用しました。

mesh2latlon.py

#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# https://www.gis-py.com/entry/py-mesh

import sys

def mesh2latlon (meshCode):
    code_first_two = meshCode[0:2] # １次メッシュ用計算
    code_last_two = meshCode[2:4]
    code_first_two = int(code_first_two)
    code_last_two = int(code_last_two)
    lat = code_first_two * 2 / 3
    lon = code_last_two + 100
    if len(meshCode) > 4: # ２次メッシュ用計算
        if len(meshCode) >= 6:
            code_fifth = meshCode[4:5]
            code_sixth = meshCode[5:6]
            code_fifth = int(code_fifth)
            code_sixth = int(code_sixth)
            lat += code_fifth * 2 / 3 / 8
            lon += code_sixth / 8
        if len(meshCode) == 8: # ３次メッシュ用計算
            code_seventh = meshCode[6:7]
            code_eighth = meshCode[7:8]
            code_seventh = int(code_seventh)
            code_eighth = int(code_eighth)
            lat += code_seventh * 2 / 3 / 8 / 10
            lon += code_eighth / 8 / 10
    return lat, lon

if __name__ == '__main__':
    if len (sys.argv) < 2:
        print ("Usage {} meshCode".format(sys.argv[0]))
        meshCode = '52396594'
        lat, lon = mesh2latlon (meshCode)
        print ("meshCode {} is converted to {:.5f} {:.5f}."
            .format(meshCode, lat, lon))
        sys.exit (1)
    lat, lon = mesh2latlon(sys.argv[1])
    print (lat, lon)

# EOF

今回必要なメッシュは、1辺10キロメートルのものなので、統計局CSVファイルを全部読み込みながら、座標を出力するPythonコードjapanloc.pyを作成します。これは、先のmesh2latlon。pyをインポートしています。

japanloc.py

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

import csv
import glob
import mesh2latlon
import os
import sys
import time
import unicodedata

fnames = glob.glob ("csv/*.csv")
fnames.sort()

for fname in fnames:
    with open (fname, 'r', encoding='CP932') as fp:
        line = csv.reader (fp)
        header = next (line)
        mcode_prev = ''
        for row in line:
            mcode = row[2]
            if mcode[0:6] == mcode_prev[0:6]: continue # 2次メッシュまで利用
            mcode_prev = mcode
            loc = row[1]
            try:
                if row[3]: loc += row[3]
            except IndexError:
                pass
            lat, lon = mesh2latlon.mesh2latlon (mcode)
            # 桁揃え https://mulberrytassel.com/python-practice-zenhan/
            # cnt = 0
            # for c in loc:
            #     if unicodedata.east_asian_width (c) in "FWA": cnt += 1
            # col = 28 - cnt
            # print ("{:{wd}s} {} {:.2f} {:.2f}"
            #   .format (loc, mcode, lat, lon, wd = col))
            print ("{} {:.6f} {:.6f}".format (mcode, lat, lon))

# EOF

コメントされているコードは、市町村名を表示するためのもので、デバッグに使いました。このコードは、2次メッシュコード、緯度、経度を1行ごとに出力します。市町村の境界には、その複数市町村名に対して同一メッシュ番号が与えられます。この重複メッシュを除くために、市町村名を表示しないようにした上で、メッシュコードの昇順にて並べ替えを行い、uniqコマンドでこの重複を取り除きます。

./japanloc.py | sort -k 1n | uniq > japanloc.txt

メッシュコードで並べ替えを行ったこの座標は、行ごとに大きく離れた場所を指しますので、第1順位を緯度に、第2順位を経度として再び並べ替えを行い、2次元平面に対して1次元走査できるようにします。私は、一旦vimに読み込んだ上で、:%! sort -k 2n -k 3nとして、対話的に結果を確認し、ZZにて結果を保存しました。

このようにして得た、日本国土の10キロメートル正方の座標ファイルjapanloc.txtの行数は11,986行です。統計局CSVファイルの先頭行にあるヘッダを除いた3次メッシュ全行数は461,373行であり、2次メッシュ抽出後の行数は12,879行でした。

メッシュコードの応用

実は、このようなことを考えたり、急に日本の国土座標に興味が湧いてきたきっかけは、東京海洋大学の高須知二先生が解析されたALESセンチメートル級測位サービスの基準局座標地図と、どこかでこの別バージョンの特定座標からの同心円にて表した地図を見たことです。距離に対するRTK（reltime kinematic）のTTFF（time to first fix）だけでなく、その方位角依存性も解析できるかもしれないと思ったからです。ALESセンチメートル級測位サービスは、本来の高精度測位ツールとしてだけでなく、研究ツールとしても、有用であることに気づきました。

そこで、基準局座標を知るために、上述の座標列をNMEA GGAとしてこのNTRIP Casterに送りながら、RTCM（Radio Technical Commission for Maritime Services）メッセージタイプ1005（座標）と1033（局番号、アンテナ名、アンテナ設定番号、受信機名やシリアル番号）を1行にまとめて表示するコードを作成しました。この表示は、RTCMメッセージに変更があったときのみ出力されるようにして、この出力を継続的に記録します。このRTCMメッセージ読み込みコードは、QZS L6 Toolsのshowrtcm.pyをインポートして作成しました。申し訳ございませんが、このコードは非公開とします。

幸いなことに、ALESセンチメートル級測位サービスでは、1秒ごとの観測データだけでなく、同一タイミングにて基準局座標も伝送されていました。ここでは、ひとつの座標送信後に1.5秒間待機して、次の座標を送信する手順をくり返しましたので、この座標列をすべて送信し終えるまでに5時間弱かかりました。この結果には、重複した基準局がリストアップされます。そこで、局番号の昇順にて並べ替えた後にuniqコマンドで重複基準局のエントリを取り除きました。

この方法により、3,086基準局を見つけました。高須先生は3,194基準局を発見されていますので、残念ながら私の結果には100局以上のスキャン漏れがあるようです。例えば、局番号172は離れた2つの基準局に重複して割り当てられていて、また、一部の基準局には受信機バージョン番号4.4.0が付加されていました。すみませんが、この生データも公表できません。ごめんなさい、もう二度とこんなことしません。

まとめ

日本国土の10キロメートル正方座標を総務省統計局のCSVファイルを用いて得ました。稠密（ちゅうみつ）に設置されたRTK基準局は、測位研究の有用なツールになり得ます。日本には、みちびきや優れた民間サービスがあるので、日本で測位研究を行い、また測位技術を学ぶ価値は高いといえます。