POSTED BY
2023-07-15

Google Mapsをはじめ現在座標表現はほとんど世界測地系(WGS84)が使われていると思われるが、まだ日本測地系(TOKYO)を採用しているシステムもあるので、測地系の相互変換が必要になるが、Cとかで使える計算式が検索しても出てこなかったのでメモ。

以下ソースです。

/*
GPS 座標操作関係
日本測地系<->世界測地系 相互変換(近似値)
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* 世界測地系->日本測地系変換 */
/* double lat,lon=度単位10進数 例=35.910722,139.459500 */
extern void com_wgs84_to_tokyo( double lat_w,double lon_w,double *lat_t,double *lon_t )
{
  *lat_t = (lat_w + 0.00010696 * lat_w - 0.000017467 * lon_w - 0.0046020);
  *lon_t = (lon_w + 0.000046047 * lat_w + 0.000083049 * lon_w - 0.010041);
}

/* 日本測地系->世界測地系変換 */
extern void com_tokyo_to_wgs84( double lat_t,double lon_t,double *lat_w,double *lon_w )
{
  *lat_w = (lat_t - 0.00010695 * lat_t + 0.000017464 * lon_t + 0.0046017);
  *lon_w = (lon_t - 0.000046038 * lat_t - 0.000083043 * lon_t + 0.010040);
}

/* TEST MAIN
gcc -D_COM_GEO_MAIN geo.c
*/
#ifdef _COM_GEO_MAIN
extern int main(int argc, char **argv) {
  double lat_t_kurume = 35.742798; //イオンモール東久留米(バス停)(日本測地)
  double lon_t_kurume = 139.529953; //イオンモール東久留米(バス停)(日本測地)
  double lat_w_kurume;
  double lon_w_kurume;
  double lat_t_musasi = 35.742694; //イオンモールむさし村山(バス停)(日本測地)
  double lon_t_musasi = 139.386639; //イオンモールむさし村山(バス停)(日本測地)
  double lat_w_musasi;
  double lon_w_musasi;

  com_tokyo_to_wgs84(lat_t_kurume, lon_t_kurume, &lat_w_kurume, &lon_w_kurume);
  fprintf(stdout,"イオンモール東久留米(バス停)日本測地: %.6lf,%.6lf 世界測地: %.6lf,%.6lf\n",
    lat_t_kurume, lon_t_kurume, lat_w_kurume, lon_w_kurume);

  com_wgs84_to_tokyo(lat_w_kurume, lon_w_kurume, &lat_t_kurume, &lon_t_kurume);
  fprintf(stdout,"イオンモール東久留米(バス停)世界測地: %.6lf,%.6lf 日本測地: %.6lf,%.6lf\n",
    lat_w_kurume, lon_w_kurume, lat_t_kurume, lon_t_kurume );

  com_tokyo_to_wgs84(lat_t_musasi, lon_t_musasi, &lat_w_musasi, &lon_w_musasi);
  fprintf(stdout,"イオンモールむさし村山(バス停)日本測地: %.6lf,%.6lf 世界測地: %.6lf,%.6lf\n",
    lat_t_musasi, lon_t_musasi, lat_w_musasi, lon_w_musasi);

  com_wgs84_to_tokyo(lat_w_musasi, lon_w_musasi, &lat_t_musasi, &lon_t_musasi);
  fprintf(stdout,"イオンモールむさし村山(バス停)世界測地: %.6lf,%.6lf 日本測地: %.6lf,%.6lf\n",
    lat_w_musasi, lon_w_musasi, lat_t_musasi, lon_t_musasi );

  return 0;
}
#endif /* _COM_GEO_MAIN */

wgs84_to_tokyo関数で世界測地系→日本測地系変換、tokyo_to_wgs84関数で日本測地系→世界測地系変換。
サンプルは以下のようにコンパイルします。

gcc -D_COM_GEO_MAIN geo.c

テスト実行

./a.out

イオンモール東久留米(バス停)日本測地: 35.742798,139.529953 世界測地: 35.746014,139.526760
イオンモール東久留米(バス停)世界測地: 35.746014,139.526760 日本測地: 35.742798,139.529953
イオンモールむさし村山(バス停)日本測地: 35.742694,139.386639 世界測地: 35.745907,139.383458
イオンモールむさし村山(バス停)世界測地: 35.745907,139.383458 日本測地: 35.742694,139.386639

Google Mapsでも使われている度単位10進数形式(35.742798,139.529953)のdouble型座標値に、地球を球体と見なした丸みの近似値を掛けることで相互変換を行います。近似なので多少はズレるものの、かなり高い精度で合いました。

世界測地:の出力(35.745907,139.383458)はそのままの形式でGoogle Mapsに貼り付け、確認可能です。
変換後の値をそのまま逆の変換にかけた結果、元の値に戻っているのが、正常動作のポイントです。