satoh

知っているようで知らないデジタル図形処理 ~任意形状図形の選択から交点・接点・接線計算まで、方程式では解けないノウハウを紹介

   

緯度経度と平面直角座標系

測地系

日本の測地系は2001年以前はベッセル楕円体という、1841年頃に提案されたものを使っていました。いわゆる日本測地系です。
それ以降はGRS楕円体を使った世界測地系です。とは言うものの、国土地理院がそう言ってるだけで古い測地系の地図がたくさんありますし、地図関係のコンテンツは高価なので古い地図そのまま使っている場合が多いと思います。
測地系の違いにより緯度経度がずれてくるので、まずは、どの測地系の地図データを使うのかを調べないとなりません。
ゼンリンの地図も今は「日本測地系」か「世界測地系」か選択できると思います。
過去と違ってGPSで宇宙から現在位置を確定できますし、全世界同一の測地系用いたほうが利便性は高いです。
必要であれば国土地理院からも日本測地系→世界測地系変換のソフトウェアがダウンロードできます。

平面直角座標系

理想的な緯度経度は世界測地系を用いてGPSからの信号で世界共通で確定できます。
物理的には理想的ですが、実際面を考えると緯度経度では成り立ちません。
「市役所から自宅までの距離は緯度方向に0.009度、経度方向に0.008度である」
と言ってもなんのことだかちっともわかりません。そこで考えられたのが平面直角座標系です。
測量とか地図からわかる距離はこの平面直角座標系に基づいて成り立っています。
日本に19ヶ所の基準点を設けて、そこを(0,0)として範囲内を平面として扱う、というものになります。基準点については国土地理院のサイトを参照してください。

緯度経度→平面直角座標系変換

GPS端末などから得られる緯度経度から位置を計算するのに使えます。
実際には前章に記したように19ヶ所の基準点の考慮が必要です。座標原点の緯度経度(基準点)でも同一計算式で求めてその差分が原点からの距離(m)となります。

    gap :子午線弧長の差
    prc  :卯酉線曲率半径
    eta  :第二離心率に子午線曲率半径の余弦を乗じたもの

    lam :経度の差
    t      :子午線曲率半径の正接

    b     : 緯度
    m0  :座標系の原点における縮尺係数

と定義すると、Cで書くと、

  //—————– X座標
    x = m0 * (
                gap
                + pow(lam, 2) * prc * sin(b) * cos(b) / 2.0
                + pow(lam, 4) * prc * sin(b) * pow(cos(b), 3) / 24.0
                * (5.0 – pow(t, 2) + 9.0 * pow(eta, 2)
                + 4.0 * pow(eta, 4))
                + pow(lam, 6) * prc * sin(b) * pow(cos(b), 5) / 720.0
                * (61.0 – 58.0 * pow(t, 2)
                + pow(t, 4)
                + 270.0 * pow(eta, 2)
                – 330.0 * pow(t, 2) * pow(eta, 2))
    ) ;

    //—————– Y座標
    y = m0 * (
                lam * prc * cos(b) + pow(lam, 3) * pow(cos(b), 3) * prc / 6.0
                * (1.0 – pow(t, 2) + pow(eta, 2))
                + pow(lam, 5) * pow(cos(b), 5) * prc / 120.0
                * (5.0 – 18.0 * pow(t, 2) + pow(t, 4)
                + 14.0 * pow(eta, 2)
                – 58.0 * pow(t, 2) * pow(eta, 2))
    ) ;

ここで注意点は、xは緯度方向ということです。南北方向です。yが東西方向。

平面直角座標系変換→緯度経度

逆の変換です。

    DEG :マクロ。ANGLE→DEGREE変換
    phi  :子午線曲率半径
    ym  :y座標を縮率で割ったもの
    prc  :卯酉線曲率半径
    mrc :子午線曲率半径
    eta  :は第二離心率に子午線曲率半径の余弦を乗じたもの
    t      :子午線曲率半径の正接
    orgL:原点経度
    b     :求める緯度。DEGREE
    l      :求める経度。DEGREE

と定義すると、Cで書くと、

//—————– 緯度
    b = DEG(phi)
        – DEG(pow(ym, 2) * t / (2.0 * mrc * prc))
        + DEG(pow(ym, 4) * t * (5.0 + 3.0 * pow(t, 2)
        + pow(eta, 2)
        – 9.0 * pow(eta, 2) * pow(t, 2)
        – 4.0 * pow(eta, 4))
        / (24.0 * mrc * pow(prc, 3)))
        – DEG(pow(ym, 6) * t * (61.0 + 90.0 * pow(t, 2)
        + 45.0 * pow(t, 4)
        + 46.0 * pow(eta, 2)
        – 252.0 * pow(t, 2) * pow(eta, 2)
        – 90.0 * pow(t, 4) * pow(eta, 2))
        / (720.0 * mrc * pow(prc, 5))) ;

//—————– 経度
    l = orgL
      + DEG(ym / (prc * cos(phi)))
      – DEG(pow(ym, 3) * (1.0 + 2.0 * pow(t, 2) + pow(eta, 2))
      / (6.0 * pow(prc, 3) * cos(phi)))
      + DEG(pow(ym, 5) * (5.0 + 28.0 * pow(t, 2) + 24.0 * pow(t, 4)
      + 6.0 * pow(eta, 2)
      + 8.0 * pow(t, 2) * pow(eta, 2))
      / (120.0 * pow(prc, 5) * cos(phi))) ;

基準点がわかればそこからの位置で緯度経度が計算できるわけです。

 - C, 座標系変換, 計算