3 ポイント 投稿者 GN⁺ 2023-08-28 | 1件のコメント | WhatsAppで共有
  • 緯度・経度座標にメートル単位の x,y 移動量を加える場合、移動距離が 数 km 以下で極点付近でなければ、単純な近似式だけでも高速に計算できる
  • 基本計算は、y方向の 111,111mを緯度1度、x方向の 111,111 * cos(latitude)mを経度1度とみなす方法で、100m北へ移動するなら 100 / 111111 度を足せばよい
  • 同じ考え方を地球半径 R=6378137 の球として dLat=dn/R, dLon=de/(R*cos(lat)) で計算することもでき、緯度51度で dn=100, de=100 なら latO=51.00089832, lonO=0.001427437 になる
  • 要求精度が 10m以内でオフセットが1kmまでなら、Aviation Formulary のようなより複雑な式を使うこともできるが、単純な平面近似でも1kmのオフセットで誤差は50m未満と見込まれる
  • 緯度によって1度の長さが変わる効果まで扱う必要があるなら、meters per degree の式を使うか、地域投影座標系に変換して移動量を加えたあと、再び緯度・経度に戻すほうが安全

短距離の移動なら 111,111m/度 の近似で十分

  • 小さな移動量は次の近似で緯度・経度の変化量を計算できる
    • y方向 111,111m ≈ 緯度1度
    • x方向 111,111 * cos(latitude)m ≈ 経度1度
  • 新しい座標はおおよそ次のように求める
    • lat_new = lat + dy / 111111
    • lon_new = lon + dx / (111111 * cos(latitude))
  • cos(latitude) には実行環境に合った単位を入れる必要がある
    • ラジアンが必要な環境では latitude * pi / 180 への変換が必要
  • この近似は、移動量が大きすぎず、極点のすぐ近くではなく、要求精度が非常に高くない場合に適している

111,111m という数字の根拠と誤差範囲

  • 111,111 という値はメートルの歴史的定義と結びついている
    • フランスがもともとメートルを、赤道から北極までパリ子午線に沿って測った距離の 10^7 分の1として定義したため
    • 10^7 / 90 = 111,111.1m が緯度1度に相当する
  • コメントでの検証では、x,y それぞれ1400m、全体変位2kmを UTM 計算と比較したとき、結果は 8.6m以下 の誤差に収まった
    • その条件で最悪の緯度は81度だった
    • 誤差は89.6度を超える地点まで10m未満に保たれた
  • 単純式は、経度が極地方へ行くほど狭くなる効果を cos(latitude) で反映している
    • 経度1度の実際の距離が短くなるため、同じ x方向のメートル移動は高緯度ほどより大きな経度変化に変換される

地球半径を使う同じ計算

  • 同じ計算は地球半径ベースの式でも表現できる
//Position, decimal degrees
lat = 51.0
lon = 0.0

//Earth’s radius, sphere
R=6378137

//offsets in meters
dn = 100
de = 100

//Coordinate offsets in radians
dLat = dn/R
dLon = de/(R*Cos(Pi*lat/180))

//OffsetPosition, decimal degrees
latO = lat + dLat * 180/Pi
lonO = lon + dLon * 180/Pi
  • この例は次の結果を返す
latO = 51,00089832
lonO = 0,001427437
  • この方法は 111,111m/度 の近似とほぼ同じ解法で、111,319.5m に近い半径ベースの値を使う点が違う
  • x方向の移動は 真の東西方向、y方向の移動は 南北方向 に近い必要がある
    • 地域投影座標系の easting/northing が回転しているなら、先に東西・南北成分へ変換する必要がある

より高い精度が必要なときの選択肢

  • Aviation Formulary の “lat/long given radial and distance” の式は、距離と方位角から新しい緯度・経度を計算するときに使える
    • 三角関数の使用を減らしたい組み込み環境ではやや複雑かもしれない
    • 距離パラメータは distance / earth radius 形式の ラジアン 値として扱う
  • 地域に適した平面座標系へ投影してからオフセットを加える方法も可能
flat_coordinate = latlon_to_utm(original_coordinate)
new_flat_coordinate = flat_coordinate + (x,y)
result_coordinate = utm_to_latlon(new_flat_coordinate)
  • この方法は UTM に限らず、その地域に適した平面座標系であれば利用できる
  • ただし、UTM ゾーン境界で移動後に別の UTM ゾーンへまたがる場合は、そのまま適用しにくい

言語別の実装例と緯度ごとの高精度式

  • Python の例は 111,111m/度 の近似をそのまま関数化している
from math import cos, radians

def meters_to_lat_lon_displacement(m, origin_latitude):
    lat = m / 111111
    lon = m / (111111 * cos(radians(origin_latitude)))
    return lat, lon
  • R の例も同じ計算を行う
deg2rad = function(deg) {(deg * pi) / (180)}

meters_to_lat_lon_displacement = function(m, origin_latitude){
  lat = m / 111111
  lon = m / (111111 * cos((deg2rad(origin_latitude))))
  return(list(lat=lat,lon=lon))
}
  • より正確な緯度別の meters per degree の式は次のように使える
meters_per_degree_lat = (111132.92 - 559.82 * np.cos(2 * lat0_rad) +
                             1.175 * np.cos(4 * lat0_rad) - 0.0023 * np.cos(6 * lat0_rad))

meters_per_degree_lon = (111412.84 * np.cos(lat0_rad) -
                            93.5 * np.cos(3 * lat0_rad) + 0.118 * np.cos(5 * lat0_rad))
  • この高精度式は、緯度によって 緯度1度と経度1度の長さ が継続的に変化する点を反映している
  • Swift の例では、緯度に応じた地球半径を計算し、距離と方位角から新しい CLLocationCoordinate2D を求める形を使っている

1件のコメント

 
GN⁺ 2023-08-28
Hacker Newsのコメント
  • メートルは1791年に、パリを通る子午線の四分円、つまり90度の弧の長さの1000万分の1として再定義された
    したがって 1° ≡ 1/90 × 10^7 m = 111,111.111... m となり、地球の周長もおよそ4000万m、つまり40,000 kmになる
    初期のメートルの定義は秒振り子、つまり周期が2秒の振り子の長さで、T ≈ 2π√(L/g) に T = 2, L = 1 を代入すると 1 = π√(1/g)、1 = π²/g となる
    だからgがπ²に近いことも単なる偶然ではなく、1 cm³の水が1 gであることも長い間グラムの定義だったためである

    • 秒はメートルの再定義より古い単位で、1日を「見栄えよく」分けたことに由来するだけに、今でも多少の偶然は残っているように見える
      メートルが秒振り子で定義されていた時代は、秒の定義とgの値に完全に結び付いており、式で書けば 1 m = 1 s² × g / π² である
      g ≈ π² は自然に出てくるが、地球の周長が40,000 kmに十分近かったため、メートルを10の累乗で再定義しても変化が大きくなかったのは偶然のように見える
      https://en.wikipedia.org/wiki/Second#Fraction_of_solar_day
    • 32ビット整数1つで、緯度や経度をおよそ1cmの精度で表現できる :D
    • メートルは実質的にパリ・フィート3つ分にも近く、約0.97mになる
      英国フィート3つ分は約0.91mにすぎない
      当時の人々が真空中で最も原理的、あるいは宇宙的に美しい長さの単位を導き出したというより、すでに使っていた単位を「そこにある棒の長さ」ではない方法で定義しようとした、というのに近い
    • GPS座標系が単にキロメートルだったらいいのにと思う
      360度の代わりに40,000 kmを使い、実際の計算では実距離を使いつつ、近似値は十分に近い
      そうすれば、少なくともメートル法の利用者には距離へ変換する必要がなくなる
      度単位の問題は、有用な距離に変換しにくいことにあり、こうしたコツは役に立つとはいえ、そもそも変換がないほうがよい
    • フランスは角度にも十進化を適用したので、実際には 1 gon = 100 km で、1 km は 1 centigon にすぎない
  • 海里1つ、約6076 ftは、地球の赤道上で正確に弧1分に相当する
    航海する立場では、すべてのマイルが海里だったらいいのにと思う
    海里には実際の意味があるが、5280 ftにはいったい何の意味があるのかと思う

    • マイルが5280 ftである理由は、80チェーンの長さだからである
      チェーンの長さは、エーカーを基準に課税していた英国の土地税法の副産物である
      ローマ・マイルは1000歩、つまり5000 ftだったので、そちらのほうが少し理にかなっていた
      https://en.wikipedia.org/wiki/Gunter%27s_chain
    • 他のほとんどの「マイル」はローマ・マイルから派生しており、フィート・ヤード・インチ・バーリーコーンのような英国単位とはある程度独立して発展したため、換算係数がおかしくなった
      元のローマ・マイルはローマ・フィート5000個分だった
      実際には 1 nmi ≡ 1.852 km と正確に定義されている
      メートルの元の定義でも 1/60 × 1/90 × 10^7 = 1851.85185185... m となる
      SIと、その前身であるMKS、CGSの中心的な特徴は、最初から単位間の換算可能性にあり、そのため 1 m ≡ 1 s ≡ 1 kg ≡ 1 N ≡ 1 Pa ≡ 1 J ≡ 1 A ≡ 1 C ≡ 1 V ≡ 1 Ω ≡ 1 F ≡ 1 W ≡ 1 Wb ≡ 1 T ≡ 1 H ≡ 1 Hz のような関係がある
      ここでの ≡ は厳密な同値ではなく、換算係数を緩く指す意味で使っている
      SIで例外に近いのはケルビン、モル、カンデラとその派生単位だが、前の2つはボルツマン定数とアボガドロ定数できれいに扱える
      カンデラがSIに入っているのは個人的に不満だ
    • 興味深いことに、マイルはもともとそれほど変ではない5000 ftだった
      ところが1500年代に英国が、当時の農業測量の計算をずっと簡単にするため、マイルを8 furlongに変えた
      https://en.m.wikipedia.org/wiki/Furlong
    • 中世の人々がマイルを5040 ftと定義してくれていればよかったのに。そうすれば 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 28, 30, 35, 36, 40, 42, 45, 48, 56, 60, 63, 70, 72, 80, 84, 90, 105, 112, 120, 126, 140, 144, 168, 180, 210, 240, 252, 280, 315, 336, 360, 420, 504, 560, 630, 720, 840, 1008, 1260, 1680, 2520 で割り切れる
    • 「海里には実際の意味がある」というときの実際が何を意味するのか気になる
      10を整数だけで3で割れないから慣用単位のほうが優れている、という類いの主張のように感じる
      円を360個の弧に分けると、焦点から特定の距離だけ離れた場所で、その弧1つに何らかの意味があると見ているようだ
      しかし約2000年前にギリシャがバビロニアの360の使用を取り入れ、バビロニアもその前の2000年間に天文学で使っていた、おおよその1年の日数の測定を洗練させてその数に到達したことを考えると、海里の意味は「実際」というより、派生的で偶然的な意味に近い
      さらに地球が扁球である点まで考慮すると、海里の長さは位置によって変わる
  • 米国に10年以上住んでいるが、いまだにヤード・ポンド法には慣れないし、これからもそうだと思う
    まったく筋が通らない
    メートル法は純金のようなもの: 1cm = 10mm, 1m = 100cm, 1km = 1000m, 1kg = 1000g, 1ton = 1000kg
    ヤード・ポンド法は「ちょっと待って」と言って、1in = ???, 1ft = 12in, 1yd = 3ft, 1mile = 5280ft, 1lb = 16oz みたいな感じ
    いったい誰がこんな狂気を作ったのかわからない

    • 米国に住んで気づいたのは、ほとんどの米国人は本当に必要でなければ単位同士を換算しないという点
      だから思ったほど問題が頻繁には表面化しない
      たまたまメートル法で表示されている場合でも、単位を換算しないのが目につく
      たとえば 1L の代わりに 1000mL、3.5kg の代わりに 3500g と書くようなもの
      欧州人なら「こちらは 600m、あちらは 1.2km」と言えるが、米国人は「こちらは 800ヤード、あちらは 1マイル」とはほとんど言わない
      欧州人なら「水を 4L 持っていく必要があるので、バッグが 4kg 重くなった」と言える
      米国人は「自分のボトルは 24液量オンスだから、だいたい 24重量オンスくらい」とは言えても、ガロンなら単に重さもガロンくらいだと言う可能性が高い
      結局、単位換算の問題は想像していたより少なかった。米国人が毎文ごとに単位換算しながら生活しているわけではないからだ
    • 米国で育った人の中にも、この体系を完全には理解していない人はかなり多いと思う
      コップ1杯の水が何オンスか、1マイルが何フィートかを知っている人が人口の50%を超えていたら驚く
      それでも米国でも科学界ではメートル法を標準として使っているので幸いだ
    • 中間の単位をいくつか飛ばしているが、それを知ると少しは説明がつく
      測量道具に由来するチェーンは 22ヤード
      1チェーンは 4 rod でもあるので、1 rod は 5½ヤードになる。実に妙だ
      10チェーンが furlong、8 furlong が 1マイル
      ちなみにエーカーは 1 furlong × 1 chain
      狂気のようでも、その中にはそれなりの体系がある
    • インチはいったい何に換算する必要があるのか?
      なぜインチをマイルに換算することがあるのか?
      生きていてインチ、あるいはフィートとインチをマイルに変えなければならないことはない
    • 単位換算も不便だが、本当に最悪なのは分数による長さ表記だと思う
      木工や手工芸では由来上、理にかなっているかもしれないが、他の用途ではどうだろう?
      ヤード・ポンド法の定規で 2 3/16" を 5.6cm と同じくらい素早く読んでみればいい
      ねじのサイズもこの影響を受けている
  • 光の速度で1ナノ秒の間に進む距離も、おおよそ 1フィート

    • 実際の値に驚くほど近い
      印象的 :)
      $ units c ft/ns の結果は * 0.98357106
    • 音の速度で1ミリ秒の間に進む距離も、おおよそ 1フィート
    • 時間を変えて、1 chrono = 光が 1e9 m を進むのにかかる時間、と定義すると格好よさそう
      1 kilochrono は 55分で、太陽日に単位を依存できない宇宙旅行のような状況ではかなり使えそう
    • 1年の秒数も π*10^7 に十分近い
    • 1 googol femtobarns も、おおよそ 1 square teraparsecs
  • 地球が扁球なら、緯度1度の実際の弧の長さは変わるのでは?
    「信頼できる」というのが単に「使える程度には近い」という意味なのか気になる
    地理と関係のない仕事を長くしすぎて、昔は知っていたことを忘れてしまったようだ

    • 信頼できるというのは、使える程度には近いという意味で、人口が多く住む地域では推定もかなり正確
      この事実を知ることになった私のユースケースは、ここにもある程度書いておいた: https://twitter.com/mholt6/status/1695685022710477043
      たとえ自分の場合に数 km 単位の誤差があったとしても、極地付近ではない可能性が高く、もし極地なら「はいはい、君は極点にいるんだね」くらいで済ませればいい
    • 地球の赤道直径は極直径より43 km大きい
      地球の軌道も同様
      学校では楕円だと習うが、実際の形の感覚はほとんど得られず、ほとんどの図は完全に間違った印象を与える
    • 記事タイトルだけを見ると、緯度1度の長さの変化が 1デシメートルにも満たないと思うかもしれないが、もちろんそうではない
      それでも多くの実用目的には十分近い
  • この記事には、111,111 * cos(latitude) m が経度1度だという良い経験則もある
    補正が気に入った
    実際には単純な定数を使ってもよく、25° は約 100,000m、44° は約 80,000m、57° は約 60,000m