師匠の散歩

VBAと測地学

地球楕円体の半径

展望計算を行ううえで地球半径の決定が不可欠でした。ネットで出てくるいろいろな楕円体の半径に関する情報を防備録としてまとめておきます。

国土地理院のページで楕円体の書公式は今は削除されており閲覧ができません。ネットで関係資料を探すと次のようなものが見つかりました。

付録6 計算式集 (平成25年3月29日一部改正)
http://psgsv2.gsi.go.jp/koukyou/jyunsoku/pdf/H25_junsoku_furoku6.pdf

Vincentyの式の導出
http://www1.kaiho.mlit.go.jp/GIJUTSUKOKUSAI/KENKYU/report/rhr50/rhr50-sa02.pdf

Wikipedia 扁球、オイラーの定理 (微分幾何学)
https://ja.wikipedia.org/wiki/オイラーの定理 (微分幾何学)

地球緯度と地心緯度との違い
http://www.infra.kochi-tech.ac.jp/takagi/Geomatics/2GeoCoord.pdf

楕円体の原子と諸公式

測量法(昭和24年法律第188号)第34条の規定に基づく作業規程の準則(平成20年国土交通省告示第413号)

楕円体の原子

  測地系:GRS80

  長半径  Semi-mejor axis
    a = 6378137 m
  
  扁平率  Flattening
    f = 1 / 298.257222101

楕円体の諸公式

  短半径  Semi-minor axis
    b = a * (1 - f)

  離心率  eccentricity
    e = sqrt(2 * f - f2)

  緯度:φにおいて 
    W = sqrt(1 - e2 * sin2(φ))
  
  とすると、曲率半径M,Nは以下のようになる。

  子午線曲率半径 Meridian radius of curvature
    M = a * (1 - e2) / W3

  卯酉線曲率半径  Prime vertical radius of curvature
    N = a / W

  平均半径
    R = sqrt( M * N )
  

曲率半径

山岳展望では緯度φにおける曲率半径Rが必要でした。さらに言えば、地点(φ1,λ1),(φ2,λ2)を通る曲率半径があるとGoodでした。

中心角φのときの楕円の半径R

  Radius_ellipsoid = a * b / sqrt( (a * sin(φ))2 + (b * cos(φ))2 )

回転楕円体と体積が同じ球の曲率半径

  Radius_elipsoid_ball = ( a2 * b ) ^ (1/3)

オイラーの式の曲率半径

緯度:φにおける、方位角:αのときの曲率半径は次の式で表される。

  Euler's function オイラーの式を使った半径
  Radius_euler = M * N / (N * cos2α + M * sin2α)  Wikipedia オイラーの定理

  地点1(φ1,λ1)、地点2(φ2,λ2)、方位角αの曲率半径の算術平均を求める。
  Radius_euler_mean = [Radius_euler_mean(φ1,α) + Radius_euler_mean(φ2,α) ] / 2

師匠の散歩の山岳展望で採用する半径

師匠が山岳展望で採用する半径は、2地点の曲率半径を平均し、展望係数をかけたものを採用しています。


  地点1(φ1,λ1)、地点2(φ2,λ2)、方位角α、展望係数κ の曲率半径
  Radius_tenbou(φ1) = M(φ1) * N(φ1) / (N(φ1) * cos2α + M(φ1) * sin2α) 
  Radius_tenbou(φ2) = M(φ2) * N(φ2) / (N(φ2) * cos2α + M(φ2) * sin2α) 

  Radius_tenbou = [ Radius_tenbou(φ1) + Radius_tenbou(φ2) ] / 2 * κ

地球上の曲率半径

測地系をGRS80としたとき、地球上の曲率半径をグラフに表すと下記のようになります。横軸は緯度[度]、縦軸は半径[m]です。

さらなる深みへ

展望の計算を突き詰めていくと、どんどん測地学の領域に足を突っ込むことになってきました。

今、私たちが使っている「標高」は回転楕円体上の一点となる3次元データの値ではなく、「重力の等ポテンシャル面」もしくは「水準面」の高さをジオイド高と呼び、回転楕円体上の高度からジオイド高を差し引いた数値になっています(国土地理院:ジオイドとは)。

山岳展望において、山の重力が大きいからといって実際にへこんでくれるわけはなく、3次元的な位置関係が重要なのですから、標高ではなくジオイド高を足した楕円体高で議論しないといけないことがわかります。国土地理院からマップデータを入手してデータ処理するには労力が大きすぎるし、あくまでも簡易的な計算で山名同定を行うのが師匠の目的ですから、今後も標高を使い続けるこにとします。

あとがき

2019年8月に 指摘していただきましたが、修正するのに半年伸びてしまいました。すみません。Wの項の2乗を修正しました。2019/12/09


Topに戻る ・・・ 戻る
Copyright(C) Grandmaster since 2010 最終更新:2019/12/09