測地系変換の精度比較：通常変換・pyproj・TKY2JGDの違い

# はじめに
地理座標には、地球全体を基準とした「世界測地系」と、日本独自の基準点や楕円体を用いる「日本測地系」があります。両者は相互に変換することができます。
本記事では、その代表的な変換手法（通常変換、pyproj、高精度変換(TKY2JGD)）をPythonで実装し、各手法で生じる誤差を比較してみます。

# 測地系の基礎知識
測地系とは、地球上の位置を表すための基準となる座標系です。
日本では、従来の日本測地系（Tokyo Datum）と、国際標準の世界測地系（WGS84/JGD2000）が使われています。
両者は基準とする楕円体や原点が異なるため、同じ地点でも座標値が異なります。
そのため、測地系をまたいだデータ利用や現場作業では、正確な変換が不可欠です。

# 測地系変換の手法
測地系変換にはいくつかの方法があります。
本記事では、以下の3つの手法を取り上げます。

- **通常変換**：全国一律の簡易補正式を用いた手法
- **pyproj**：座標変換ライブラリ（PROJ）を利用した手法
- **高精度変換（TKY2JGD）**：国土地理院が提供する地域ごとの補正値を用いた高精度な手法

それぞれの手法で、変換後の座標精度にどのような違いが生じるかを比較します。

# 測地系変換
## 検証地点
- 変換誤差が大きく現れる離島の代表例として「石垣島」を検証地点とします
- 日本測地系座標から世界測地系座標に変換します

| 測地系 | 緯度 | 経度 |
|:-----------:|:------------:|:------------:|
| 日本測地系 | 24.3344444444 | 124.1561111111 |
| 世界測地系 | 24.339932361 | 124.156508419 |

:::note
離島や本土から遠く離れた地域（小笠原諸島、南北大東島、沖縄の一部、北海道北端の島々など）では、簡易式変換の誤差が数十メートルから200メートルを超えることがあります
:::

## 通常変換
日本測地系→世界測地系の簡易式（全国一律補正）を用いて変換します

```python
def tky2wgs(lat, lon):
    lat_w = lat - lat * 0.00010695 + lon * 0.000017464 + 0.0046017
    lon_w = lon - lat * 0.000046038 - lon * 0.000083043 + 0.010040
    return lat_w, lon_w

lat_tky, lon_tky = 24.3344444444, 124.1561111111

lat_wgs, lon_wgs = tky2wgs(lat_tky, lon_tky)
print(lat_wgs, lon_wgs)
# 出力：24.338611837891115 124.15472050601167
```

## pyproj
pyprojライブラリを用いた測地系変換の例です

```python
from pyproj import Transformer

lat_tky, lon_tky = 24.3344444444, 124.1561111111

# EPSG:4301 = 日本測地系, EPSG:4326 = WGS84
transformer = Transformer.from_crs("EPSG:4301", "EPSG:4326", always_xy=True)
lon_wgs, lat_wgs = transformer.transform(lon_tky, lat_tky)
print(lat_wgs, lon_wgs)
# 出力：24.338624537501115 124.15449819512936
```

## 高精度変換（TKY2JGD）
国土地理院が提供するTKY2JGDパラメータファイルを利用し、地域ごとの補正値を適用した高精度変換を行います

**座標変換パラメータ・ファイル**
[ファイル配布サイト](https://www.gsi.go.jp/sokuchikijun/tky2jgd_download.html)

```python
import tky2jgd

lat_tky, lon_tky = 24.3344444444, 124.1561111111

tky2jgd.load_parameter("TKY2JGD.par")
dB, dL = tky2jgd.bilinear(lat_tky, lon_tky)
lat_jgd = lat_tky + dB / 3600
lon_jgd = lon_tky + dL / 3600
print(lat_jgd, lon_jgd)
# 出力：24.339932361951444 124.1565084197749
```

# 誤差の比較
各変換手法で得られた世界測地系座標と、国土地理院の基準点成果値（世界測地系）との距離（地表誤差）を比較します。

```python
from math import radians, sin, cos, sqrt, atan2

def haversine(lat1, lon1, lat2, lon2):
    R = 6378137  # 地球半径[m]
    dlat = radians(lat2 - lat1)
    dlon = radians(lon2 - lon1)
    a = sin(dlat/2)**2 + cos(radians(lat1)) * cos(radians(lat2)) * sin(dlon/2)**2
    c = 2 * atan2(sqrt(a), sqrt(1-a))
    return R * c

# 例
lat_ref, lon_ref = 24.339932361, 124.156508419  # 基準点（世界測地系）
# 各手法で得た座標をlat_calc, lon_calcとして
distance = haversine(lat_ref, lon_ref, lat_calc, lon_calc)
print(distance)
```

| 変換手法               | 緯度           | 経度           | 誤差（m） |
|------------------------|---------------|---------------|-----------|
| 通常変換               | 24.338611837891115 | 124.15472050601167 | 233.43776301573956 |
| pyproj                 | 24.338624537501115 | 124.15449819512936 | 250.53078177987675 |
| 高精度変換 (TKY2JGD)   | 24.339932361951444 | 124.1565084197749 | 0.00013188913222721108 |


# まとめ
- 測地系変換は、手法によって誤差が大きく異なります
- 通常変換やpyprojでは、離島など本土から離れた地域で数百メートル単位の誤差が発生する場合があります
- 国土地理院の高精度変換（TKY2JGD）を用いることで、誤差をほぼゼロに抑えることが可能です
- 測地系変換を行う際は、用途や地域に応じて適切な手法を選択することが重要です