はじめに
みなさん。平面直角座標系はご存知ですか？
こんな記事をご覧の皆様は既にご存知かもしれませんが、日本の公共測量などに利用されている座標系です。
（わかりやすい平面直角座標系：https://www.gsi.go.jp/sokuchikijun/jpc.html ）
平面直角座標系
「平面直角座標系」なんて単語は一般の方は絶対知らないと思うので、簡単に説明すると「日本を19のゾーンに分割し、各ゾーンに座標原点を設けて平面として捉えた座標系」です。
は？って感じですよね。
もう少し詳しく説明すると、「地球は楕円体なので、平面（例えば画面や紙）の地図に起こすと何かしらの要素が必ず歪んでしまうため、歪みを最小限に抑えるべく狭い範囲に限定して平面に投影する手法」と思っていただければ良いと思います。
狭い範囲であれば歪んでしまっても、まぁ誤差みたいなもんだし問題ないよね！っていう考えで、日本で生み出されて、国内では広く利用されている一般的な手法です。
例えば国土地理院の位置は平面直角座標9系で
X座標：11543.6883m
Y座標：22916.2436m
の位置にあります。
「国土地理院は平面直角座標9系でX座標が11543.6883、Y座標が22916.2436の位置にあるよ！！！」って元気に言われてもキレますよね？
工事や測量の図面は歪みが最小限でなければいけないので平面直角座標系が多用されますが、この図面に記載の建物は地球上のどこにあるんだろう？緯度経度でいったらどの辺？と疑問に思うことはよくあります。
現場に向かう際にも緯度経度でわからないとGoogleマップで検索できないですよね。
（まぁ実際にはカーナビなどに住所を入力して検索してしまうと思いますが…）
なので、緯度経度を平面直角座標系に変換するツールは国土地理院が用意してくれており、Web上から簡単に変換することができます。
(https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/xy2blf.html)
また、逆に緯度経度を平面直角座標系に変換するツールというのも親切に用意してくれています。
（https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xyf.html）
実際、国土地理院は緯度経度で表すと
緯度：36°6'13.58925"
経度：140°5'16.27815"
に存在していますが、これを平面直角座標9系で示すと
X座標：11543.6883m
Y座標：22916.2436m
とXY座標に変換できます。（単位はm）
度分秒でも十進法でも入力できるし、平面直角座標から緯度経度も割り出せたりしてとても便利なんですが、このツールには一つ難点がありまして、「そもそも指定した緯度経度がどこの座標系に属しているのか把握していなければならない」という点です。
緯度経度から平面直角座標系を割り出してみる
例えば先程の国土地理院の緯度経度はあらかじめ9系であることがわかっていたので、9系を想定して変換を行ったため、正確に座標を割り出すことができました。
- 9系/360613.58925/1400516.27815
しかし系番号が不明な場合はどうでしょうか？
当然、計算することは不可能ですし、適当な座標系を入力すると歪みが大きすぎるため計算不能になります。
- 16系/360613.58925/1400516.27815
長年、GIS業界に関わられている方であれば緯度経度から座標系を想定したり、前述の国土地理院のホームページなどから判別したりすることが可能ですが、普段あまりGISに関わらない方であればそもそも何が間違っているかもわからずGoogleで検索することすら困難なのではないでしょうか？
簡単に調べてみたところ、座標系を識別するツールみたいなものはなさそうな感じでした…
と、いうことで
ないなら作っちゃえばいいじゃない！！！ということで作成してみました！
https://github.com/MIERUNE/search_zone_number
こちら、Pythonパッケージとしても利用可能なCLIツールで、「市区町村名・緯度経度」から平面直角座標系を取得できます！
詳しくはリポジトリの方をご覧ください！というわけにもいかないので簡単に説明していきますね。
平面直角座標の情報を持つデータを作成
今回、地名と緯度経度から座標系を割り出すツールを作りましたが、割り出すために必要な元データは、100MBを超えてくるファイルサイズになっているため、リポジトリの中に格納するのは躊躇われます。
なので、生成用のツールも作成しました！
https://github.com/MIERUNE/create_gpkg_for_city_boundaries
まずはこちらを使ってデータを生成していきましょう。
とはいえリポジトリをcloneして、以下のコマンドを入力するだけで生成が可能です。
% cd scripts/
% ./run.sh
生成にはとても時間がかかるので気長にお待ちください…
すべての処理が終了して、merge_city_boundary.featherが生成されていればOKです！
インストール
こちら（https://github.com/MIERUNE/search_zone_number）のリポジトリに戻ります。
まずはコードをcloneしておきましょう！
ただ、ツールを利用する際に、自身の環境にPipenvがインストールされていないと利用できないので、されていない方はまずPipenvの構築からお願いします！
こちらの記事などがよくまとまっていてわかりやすいかなと思います。
リポジトリのルートに移動したら、パッケージをインストールしましょう
pipenv sync
その後、先程生成したmerge_city_boundary.featherを.../search_zone_number/search_zone_number/assets/に格納してください。
これで準備完了です！
使い方
緯度経度で検索
緯度経度から市区町村名と系番号を取得するには、lnglatサブコマンドを利用します。
% pipenv run python -m search_zone_number lnglat --help
Usage: __main__.py lnglat [OPTIONS]
緯度経度から対象の行を取得
Options:
-n, --lng FLOAT 経度を入力 [required]
-a, --lat FLOAT 緯度を入力 [required]
--help
早速使ってみましょう！
以下の様に緯度経度を指定することで、市区町村名と平面直角座標系の番号を取得できます！
% pipenv run python -m search_zone_number lnglat --lng 141.34694 --lat 43.06417
GST_CSS_NAME system_number
1042 札幌市中央区 12
市区町村名で検索
市区町村名を入力し、系番号を取得することも可能です。
% pipenv run python -m search_zone_number name -n 札幌市
GST_CSS_NAME system_number
1042 札幌市中央区 12
1043 札幌市北区 12
1044 札幌市南区 12
1045 札幌市厚別区 12
1046 札幌市手稲区 12
1047 札幌市東区 12
1048 札幌市清田区 12
1049 札幌市白石区 12
1050 札幌市西区 12
1051 札幌市豊平区 12
終わりに
ということで、市区町村名や緯度経度を入力するだけで平面直角座標系の番号をさっと割り出すことのできるツールを作成していきました！
これで、投影変換などの作業をスムーズに行うことができますね！
ただ、現状２つのツールを行き来したり、重たいデータを作成する必要があったりと色々手間なのは承知しております…！！！
なので、Webブラウザーから地図をクリックするだけで簡単に系番号を閲覧できるツールも作りました！
こちらも近いうちに公開しようと思っているのでお待ちください！