日本測地系から世界測地系に変換する(けっこうガチで計算するVer)

ruby_map_service.001

先月、世界測地系から日本測地系に変換するエントリを書いた。
世界測地系から日本測地系に変換する
けど、その式ではめっちゃ誤差が出てしまう事象が起こった。
(日本測地系→世界測地系をやってる時でした。)
元データの緯度経度がどういうわけか、普通(Googleとかでもらえるやつ)よりも2桁くらい少ないものだったが、当初は
「まぁけど普通に使えるっしょ」
と、たかをくくって例の近似式で変換して使おうとしたところ
700m前後のズレが発生
これはマズイ・・ と、思い原因をさぐるも 普通に考えて元データがダメなんだろ? と疑っていた。
ところが、実はガチで計算するとドンピシャでズレがおさまるという事を知る。

直行座標変換

ガチで計算すると言っても、最初そんな計算式が存在する事を微塵も知らなかったのだが、会社の人に
「このPHPのソースではうまくいってるらしいよ」
と、ソースコードだけ渡されたので はぁ・・ と思い中を見てみると、けっこうややこしい式がツラツラと。

そのコードを書いた人は既に会社にいないのか知らないが、とにかくドキュメントはソースです 状態だったので、何をやっているのか解剖すると共に、それをRubyに移植した。

計算のやり方は、おそらくこの「師匠の散歩」の中で紹介されている直行座標変換を利用する式を使っていると思われる。
(けど計算式全部のってるわけじゃないので、どうやって調べたんだろう〜〜)

ざっくりした流れはこう。

1.元々の日本測地系で使用しているベッセル楕円球の各定数を使いながら緯度経度を一度X,y座標に変換
2.x,y をそれぞれ世界測地系の方へ移動(移動距離は一定)
3.今度は世界測地系の定数(WGS84)を使いながら再び緯度経度に変換

試してみたところ完璧にズレが収まったのでビビりながらも、PHPのソースを書いたどっかの人には感謝を禁じ得ない。

Ruby移植版

ソースは最終的にこんな具合になった。

module JapanGeodeticSystem
  class << self
    RD = Math::PI / 180 # ラジアン(度)

    # 日本測地系の定数(ベッセル楕円体)
    R_JP  = 6377397.155             # 赤道半径
    F_JP  = 1 / 299.1528128         # 扁平率
    E2_JP = 2 * F_JP - F_JP * F_JP  # 第一離心率

    # 世界測地系の定数(WGS84)
    R_WS  = 6378137.0               # 赤道半径
    F_WS  = 1 / 298.257223563       # 扁平率
    E2_WS = 2 * F_WS - F_WS * F_WS  # 第一離心率

    # 並行移動量(m)
    DX = -148.0; DY = 507.0; DZ = 681.0

    # 楕円体高
    HEIGHT = 0

    # 日本測地系から世界測地系に変換する
    def convert_into_wgs(lat, lng)
      x, y, z = llh2xyz(lat, lng, HEIGHT, R_JP, E2_JP, RD)
      x += DX; y += DY; z += DZ
      lat, lng, h = xyz2llh(x, y, z, R_WS, E2_WS, RD)
      [lat, lng, h]
    end

    # 座標系の変換(緯度経度 -> xyz)
    def llh2xyz(lat, lng, h, a, e2, rd)
      lat *= rd
      lng *= rd
      sb = Math.sin(lat)
      cb = Math.cos(lat)
      rn = a / Math.sqrt(1 - e2 * sb * sb)
      x = (rn + h) * cb * Math.cos(lng)
      y = (rn + h) * cb * Math.sin(lng)
      z = (rn * (1 - e2) + h) * sb
      [x, y, z]
    end

    # 座標系の変換(xyz -> 緯度経度)
    def xyz2llh(x, y, z, a, e2, rd)
      bda = Math.sqrt(1 - e2)
      p = Math.sqrt(x * x + y * y)
      t = Math.atan2(z, p * bda)
      st = Math.sin(t)
      ct = Math.cos(t)
      b = Math.atan2(z + e2 * a / bda * st * st * st,
                      p - e2 * a * ct * ct * ct)
      l = Math.atan2(y, x)
      sb = Math.sin(b)
      rn = a / Math.sqrt(1 - e2 * sb * sb)
      h = p / Math.cos(b) - rn
      l1 = b / rd
      l2 = l / rd
      [l1, l2, h]
    end
  end
end

# 使い方
JapanGeodeticSystem.convert_into_wgs({日本測地系の緯度}、{おなじく経度}) => lat, lng

ちなみにHEIGHT(楕円の高さ)が0だったので、(なんとなく理系の性か知りませんが)ホントにいいのかな?と思ってたんですが、たぶん国土地理院の提供してる変換プログラムのQAに書いてるのがアンサーっぽくて、100mずれても誤差6mmらしい。

というわけでこのModuleを使って事なきを得ました。
めでたしめでたし。

 

2016-05-05 | Posted in RubyNo Comments » 


関連記事

Comment





Comment



*