【python】WGS84から平面直角座標系への変換

はじめに

測量では成果を緯度経度でなく、平面直角座標系という２次元の平面の座標値で表現します。その相互変換を行う方法のメモです。

日本では現在、19の直角座標系が用いられています。(「分かりやすい平面直角座標系」from 国土地理院)

よくある(?)誤解: 慣性航法で使う座標系とは違う！

座標系の変換というと簡単に聞こえますが、これは

• WGS84の楕円体の極座標(緯度,経度,楕円体高)を地球中心の座標系(ECEF)の直交座標系の座標値XYZに変換
• ECEFの直交座標系の値を、該当する楕円体における接平面に対応するように回転させる

という話とは『異なり』ます。ガウス・クリューゲル投影法という計算方法を使うようなのですが、私は数式が複雑で着いていけませんでした。

当初、自分は上記のような単純な回転計算と誤解していて、周りに「そんなの簡単な計算だよ」などと話していて、思い返すと恥ずかしいです。

国土地理院のサポート

変換するために、国土地理院から提供されているものを利用できます。（さすが公共機関！）

まず、ブラウザに直接座標値を入力して変換することができます。

また、WEB APIからも利用できます。こちらに、提供されている機能と説明があります。アクセスが集中しないように同一アドレスから10秒間で10回の制限がある等の注意も書かれています。

平面直角座標系の番号を指定して取り出す例です。

use_gsi_webapi.py

import requests
from numpy import round
lat, lon = 35.71, 139.74

url = "http://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xy.pl?"
params={'latitude':round(lat, 8), 'longitude': round(lon, 8), "refFrame":2, "zone": 9, 'outputType':'json'}

res =　requests.get(url, params=params)
if res.status_code == requests.codes.ok:
    print(res.json())

ここでは次の辞書が返されます。

{'OutputData': {
  'publicX': '-32170.0990',
  'publicY': '-8445.1361',
  'gridConv': '0.054477778',
  'scaleFactor': '0.99990088'}}

EPSG

日本の平面直角座標系の番号は日本が定めるものですが、国際的にはEPSGとして定めらています。国際的には平面直角座標系の番号ではなくEPSGで参照します。(例えば写真測量のソフトであるmetashape などを使う場合とか。)

EPSG	測地系
EPSG:4326	WGS 84
EPSG:6668	JGD2011
EPSG:6669	JGD2011 / Japan Plane Rectangular CS I
EPSG:6670	JGD2011 / Japan Plane Rectangular CS II
EPSG:6671	JGD2011 / Japan Plane Rectangular CS III
EPSG:6672	JGD2011 / Japan Plane Rectangular CS IV
EPSG:6673	JGD2011 / Japan Plane Rectangular CS V
EPSG:6674	JGD2011 / Japan Plane Rectangular CS VI
EPSG:6675	JGD2011 / Japan Plane Rectangular CS VII
EPSG:6676	JGD2011 / Japan Plane Rectangular CS VIII
EPSG:6677	JGD2011 / Japan Plane Rectangular CS IX
EPSG:6678	JGD2011 / Japan Plane Rectangular CS X
EPSG:6679	JGD2011 / Japan Plane Rectangular CS XI
EPSG:6680	JGD2011 / Japan Plane Rectangular CS XII
EPSG:6681	JGD2011 / Japan Plane Rectangular CS XIII
EPSG:6682	JGD2011 / Japan Plane Rectangular CS XIV
EPSG:6683	JGD2011 / Japan Plane Rectangular CS XV
EPSG:6684	JGD2011 / Japan Plane Rectangular CS XVI
EPSG:6685	JGD2011 / Japan Plane Rectangular CS XVII
EPSG:6686	JGD2011 / Japan Plane Rectangular CS XVIII
EPSG:6687	JGD2011 / Japan Plane Rectangular CS XIX

情報源ですが、本家本元や管理団体などよく理解していません。以下がわりとauthority のような気がします。

• epsg.io で説明を読むことができます。例えばWGS84とか。
• EPSG Geodetic Parameter Registryで検索できます
• spatial reference というページにもリストや検索機能があります。

pyproj　を使う

さきほどはWEB API をりようしましたが、これはpython ライブラリを用いてもできます。
インストールで苦労した記憶はありません。

test_webapi.py

from pyproj import Transformer

lat, lon = 35.71, 139.74
wgs84_epsg, rect_epsg = 4326, 6677
tr = Transformer.from_proj(wgs84_epsg, rect_epsg)
x, y = tr.transform(lat, lon)
print(x, y)

実行すると以下を得ます。

-32170.09900022997 -8445.136058616452

有効数字は分かりませんが、地理院の結果とは一致しています。（小数点以下４桁、0.1ミリメートル！）

CRSの利用

巷では、CRSを使う方法が紹介されていましたが、利用してみると、、、

test_use_pyproj.py

import pyproj
wgs84=pyproj.CRS("EPSG:4326")
jgd2011_9 = pyproj.CRS("EPSG:6677")
lat, lon = 35.71, 139.74
print( pyproj.transform(wgs84, jgd2011_9, lat, lon) )

実行すると警告が出ますが、同じ結果を得ます。

/usr/bin/ipython3:1: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
(-32170.09900022997, -8445.136058616452)

これは多分、古い方法なのだと思います。しっかりチュートリアルを読まなければ！
https://pyproj4.github.io/pyproj/stable/examples.html

まとめ・雑感

• WGS84と平面直角座標系の変換は地理院提供のＡＰＩでできるが、pyproj でも同じことはできそう
• pyproj は多機能で4Dtansform など気になる機能もあり、余力があれば調べてみたい

等々(2020/08/08)