kamoland - ページ 5 - KMIソフトウェアとしてソフトウェア開発を行っています

世界測地系から日本測地系への変換の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) 三次元直交座標系に変換して平行移動

ダウンロード版TKY2JDG (パラメータファイル無し)
proj4を使った変換(*1)

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

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

１. 緯度，経度から三次元直交座標系(x,y,z)に変換(Bessel楕円体使用)

２. x,y,zを平行移動する(*2)

３. 三次元直交座標系から緯度，経度に変換(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.414ｍ、
地球重心から緯度ゼロ、経度90度の方向に、+507.337ｍ、
地球重心から北極方向に+680.507ｍ

とのことです．これが先ほどのproj4の設定値と対応しています．

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

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

Web版 TKY2JGD

ダウンロード版TKY2JD (パラメータファイルあり)

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

実装例

(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;

としたのですが，この平行移動量については，以下のサイトによると，

[GIS] 日本測地系からWGS84への変換パラメータの謎を追え – ここのことはなかったことにするかも

-146.336, 506.832, 680.254

という話もあるようです．

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

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