世界測地系から日本測地系への変換のJava実装

世界測地系と日本測地系の変換

拙作のAndroidアプリ「地図ロイド」では,世界測地系の緯度経度で位置を扱っているのですが, 日本測地系の対応を入れるかどうか検討しています.

それで,日本測地系に対応する方法について調べた結果をまとめます.

変換の仕様

だいたい,以下の3種類の方法があるようです.

(1) 1次式による近似

ネットを検索すると出てくる,ポピュラーな方式です. 具体的には,例えば日本測地系から世界測地系への変換の場合,単位は度で
  • lat(世界) = lat(日本) - lat(日本) * 0.00010695 + lon(日本) * 0.000017464 + 0.0046017;
  • lon(世界) = lon(日本) - lat(日本) * 0.000046038 - lon(日本) * 0.000083043 + 0.010040;

という式で変換を行うという話です.

(2) 三次元直交座標系に変換して平行移動

などで使われている方式.

具体的には,例えば日本測地系から世界測地系への変換の場合,

  1. 緯度,経度から三次元直交座標系(x,y,z)に変換(Bessel楕円体使用)
  2. x,y,zを平行移動する(*2)
  3. 三次元直交座標系から緯度,経度に変換(GRS80楕円体使用)

で求める方法です.

(*1)proj4というのは,proj4jsで以下のような設定で使われているような場合の話です.

proj4.defs([
    ["EPSG:4301",
    "+proj=longlat +ellps=bessel +towgs84=-146.414,507.337,680.507,0,0,0,0 +no_defs"
    ]
]);

(*2)平行移動の内容は, 我が国の測地基準系の改訂について 測地学研究連絡委員会報告 によると,

地球重心から緯度、経度ゼロの方向に-146.414m、
地球重心から緯度ゼロ、経度90度の方向に、+507.337m、
地球重心から北極方向に+680.507m
とのことです.これが先ほどのproj4の設定値と対応しています.

(3) 地域ごとの変換パラメータで変換する

これが正式な変換なんですかね.

の方式です.ただ,自分では実装できなさそうです.

実装例

(1) 1次式による近似

「日本測地系 世界測地系 javascript」あたりでぐぐると出てきますので,割愛.

(2) 三次元直交座標系に変換して平行移動

変換の計算内容としては,

になります.これに平行移動を組み合わせたもののJava実装です.

// 楕円体定数.A=長半径(m), F=扁平率の逆数(1/f)
// https://www.gsi.go.jp/sokuchikijun/datum-main.html#p1
private static final double A_BESSEL = 6377397.155;
private static final double F_BESSEL = 299.1528128;
private static final double A_GRS80 = 6378137;
private static final double F_GRS80 = 298.257222101;

// 平行移動量(m)
private static final double X_SHIFT = -146.414;
private static final double Y_SHIFT = 507.337;
private static final double Z_SHIFT = 680.507;

// 世界測地系から日本測地系に変換する
public static double[] worldDatumToTokyoDatum(double latD, double lonD, double height) {
    double[] xyz = latLonHeightToXYZ(latD, lonD, height);

    double[] latLonD = xyzToLatLon(xyz[0] - X_SHIFT, xyz[1] - Y_SHIFT, xyz[2] - Z_SHIFT);

    return new double[]{latLonD[0], latLonD[1]};
}
// (緯度,経度,ジオイド高+標高)から,地心直交座標へ変換する
private static double[] latLonHeightToXYZ(double latD, double lonD, double height) {
    double lat = Math.toRadians(latD);
    double lon = Math.toRadians(lonD);

    double A = A_GRS80;
    double f = F_GRS80;
    // 離心率の2乗
    double e2 = (2 * f - 1) /f / f;
    // 卯酉線曲率半径
    double N = A / Math.sqrt(1 - e2 * Math.pow(Math.sin(lat), 2));

    double x = (N + height) * Math.cos(lat) * Math.cos(lon);
    double y = (N + height) * Math.cos(lat) * Math.sin(lon);
    double z = (N * (1 - e2) + height) * Math.sin(lat);
    return new double[]{x, y, z};
}
// 地心直交座標から,緯度経度へ変換する
private static double[] xyzToLatLon(double x, double y, double z) {
    double A = A_BESSEL;
    double f = F_BESSEL;
    double e2 = (2 * f - 1) /f / f;
    double P = Math.sqrt(x*x + y*y);

    double lat0 = Math.atan2(z, P);
    // 反復計算.無限ループ防止のため上限回数を設けておく
    for (int i = 0; i < 10; i++) {
        // 卯酉線曲率半径
        double N0 = A / Math.sqrt(1 - e2 * Math.pow(Math.sin(lat0), 2));
        double latTmp = Math.atan2(z, P - e2 * N0 * Math.cos(lat0));
        if (Math.abs(latTmp - lat0) < 0.000000000001) {
            break;
        }
        lat0 = latTmp;
    }
    double lon0 = Math.atan2(y, x);

    return new double[]{Math.toDegrees(lat0), Math.toDegrees(lon0)};
}
もし日本測地系から世界測地系へ変換する場合は,楕円体の定数(ベッセル,GRS80)を逆にして,平行移動量の符号を逆にする感じです.

引数のheightとして,その位置のジオイド高+標高 を渡す必要があるのが,使い方によっては厄介かもしれません.

ちなみに ダウンロード版TKY2JDG のソース(modTky2jgd.bas の Sub ITRF94toTokyo97())では, このheightがわからない前提での実装がされています.

平行移動量の補足

上記の実装では,
private static final double X_SHIFT = -146.414;
private static final double Y_SHIFT = 507.337;
private static final double Z_SHIFT = 680.507;
としたのですが,この平行移動量については,以下のサイトによると,
-146.336, 506.832, 680.254
という話もあるようです.

この値で試したところ,わずかにWeb版TKY2JGDの結果に近い値になりました. (わずかというのは,緯度の10進数で小数点以下7桁目あたりの違い)

ですので,Web版TKY2JGDに近づけるという点では,こちらの-146.336〜 の方がよいのかなと思います.

kamolandをフォローしましょう


© 2021 KMIソフトウェア