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

## はじめに

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

日本では現在、19の直角座標系が用いられています。([「分かりやすい平面直角座標系」from 国土地理院](https://www.gsi.go.jp/sokuchikijun/jpc.html))

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

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

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

という話とは『異なり』ます。[ガウス・クリューゲル投影法](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%E3%83%BB%E3%82%AF%E3%83%AA%E3%83%A5%E3%83%BC%E3%82%B2%E3%83%AB%E5%9B%B3%E6%B3%95)という計算方法を使うようなのですが、私は数式が複雑で着いていけませんでした。

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

### 国土地理院のサポート

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

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

https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/bl2xyf.html

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

https://vldb.gsi.go.jp/sokuchi/surveycalc/main.html

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

```python: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](https://epsg.io/) で説明を読むことができます。例えば[WGS84](https://epsg.io/4326)とか。
- [EPSG Geodetic Parameter Registry](http://www.epsg-registry.org/)で検索できます
- [spatial reference](https://spatialreference.org/ref/epsg/) というページにもリストや検索機能があります。

# pyproj　を使う

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

```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を使う方法が紹介されていましたが、利用してみると、、、

```python: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)