日本測地系からWGS84へ変換 #Python

こんにちは。
日本測地系からWGS84へ変換するプログラムを作りました。

測地系 WGS84¹ から Tokyo Datum へ3パラメータ変換² ³で変換された日本測地系⁴ 経緯度座標値（多くのGPS受信機の日本測地系出力）を、WGS84へ戻す変換に相当します。
なお高度値も用いるならば、同じく3パラメータ変換で戻すのは容易ですが、ここでは高度値不要の変換を求めることとし、その代わりに内部でその地点の地表面高度を算出しています（楕円体高はパラメータファイルを入手して計算）。
変換式は近似多項式の形で求めました。この2次式の変換の正確度は下記のように島嶼を含む陸域で概ね 80cm 以下です（これに対し 1次式では 10m を超える地点が生じます）。

tky2wgs_approx.py

# angle units in degrees
def tky2wgs_approx(lon_tky, lat_tky, degree=2):
    L = lon_tky - 135
    B = lat_tky - 35
    if degree == 1:
        dL = -2.80641517e-3 - 4.98891541e-5*B - 8.33469263e-5*L
        dB =  3.20321021e-3 - 1.15452990e-4*B + 2.47032059e-5*L
    elif degree == 2:
        a00,a01,a02,a03,a04,a05,a10,a11,a12,a13,a14,a15 = -2.79648156e-03, -3.64571151e-05, -1.00958714e-06, -8.83091873e-05, -8.33002662e-07, 3.54561248e-07, 3.19774649e-03, -1.13997082e-04, -7.51313530e-07, 2.34640477e-05, 6.06731757e-07, 3.28340222e-07
        dL = a00+a01*B+a02*B**2+a03*L+a04*B*L+a05*L**2
        dB = a10+a11*B+a12*B**2+a13*L+a14*B*L+a15*L**2
    return lon_tky + dL, lat_tky + dB

WGS84 とは米国国防省策定の World Geodetic System （世界測地系）（1984年版）を指します。 ↩
「日本測地系と JGD2000 との変換について」の中の (3) です。 ↩
「日本測地系における離島位置の補正量」（測地学会誌）の後半にも説明があります。 ↩
なお 3パラメータ変換による変換だけでは正しくは日本測地系と呼ぶべきではなく、日本測地系の内部歪みも付加する必要があります。すなわち正しい日本測地系とは「TKY2JGD」（国土地理院）で取り扱われているものです。 ↩