地球上の点と線分の距離

第53回　座標にもっとも近い線分上の点を内積で求める｜gihyo.jp … 技術評論社

点$P(x_P, y_P)$と、線分$AB$（$A(x_A, y_A), B(x_B, y_B)$）の距離を求める。

• Pからおろした垂線の足がAB上にあれば垂線の長さ
• 外分点にあれば、AかB、近いほうの端点とPとの距離

なので、まずはその判定が必要になる。

ベクトルの内積を使う。$\vec{AP}$と$\vec{AB}$のなす角を$\theta(0^\circ<\theta<180^\circ)$として、

$$ \vec{AP}\cdot\vec{AB}=|\vec{AP}||\vec{AB}|\cos{\theta} $$

また、PからAB(の延長線上)におろした垂線の足をCとすると、以下の式も成立する。

$$(AC)=|\vec{AP}|\cos{\theta}$$

ここで(AC)は、$\vec{AB}$ 方向に沿った線分ACの距離だが、負になり得る値である。（そのため$|\vec{AC}|$とは表記していない）

負の時というのは、$|\vec{AP}|>=0$より$\cos{\theta}<0$、つまり$\theta>90^\circ$なので、CはAを挟んでBと逆方向にあるということである。逆に正なら同じ側にある。

2式を合わせ、さらに正規化し、パラメータtを得る。

$$ (AC) = \frac{\vec{AP}\cdot\vec{AB}}{|\vec{AB}|} \\ t = \frac{(AC)}{|\vec{AB}|} = \frac{\vec{AP}\cdot\vec{AB}}{|\vec{AB}|^2} $$

• $t<0$のとき、CはAに近い外分点。APが最短距離。
• $t>1$のとき、CはBに近い外分点。BPが最短距離。
  • $|\vec{AB}|<|\vec{AC}|$
• その他のとき、Cは線分ABの内分点。$C = A + (t * \vec{AB})$で求められ、CPが最短距離。

Python3 による実装例

ax, ay = 0, 0
bx, by = 10, 10
px, py = 9, 1

abx = bx - ax
aby = by - ay
apx = px - ax
apy = py - ay
t = (apx * abx + apy * aby) / (abx ** 2 + aby ** 2)
if t <= 0:
    cx = ax
    cy = ay
elif 1 <= t:
    cx = bx
    cy = by
else:
    cx = ax + t * abx
    cy = ay + t * aby

dist = ((px - cx) ** 2 + (py - cy) ** 2) ** 0.5

前置きはここまで。

緯度経度で示された2点A,Bを結ぶ線と、点Pとの地球表面上の距離を求めたい時。

線分ABの結び方は、複数考えられる。今回考慮するABは、厳密には等角航路だが、地球全体からすると非常に短いため(せいぜい1km)、どれを選んでも誤差は無視できるとする。

• 等角航路
  • メルカトル図法上での直線
• 大圏航路
  • 楕円球上の最短距離。正距図法での中心からの直線。「球心,A,Bの3点をとおる平面」と「楕円球表面」との共通部分
• 地球表面を無視して3次元上でABを結ぶ直線
• その他何らかの図法で平面地図に投影した時の、ABを結ぶ線分

国土交通省では地図を精確に表せるよう、ガウス・クリューゲル図法 - Wikipediaを用いた平面直角座標系を定義している。

これを利用して平面に落とし込み、平面上で上記の計算を行うと点と線分の距離が求まる。

• メリット
  • 計算が直感的に分かりやすい
  • ABを結ぶ線分が本来の意味を持つ
  • A,Bの変換後の座標を記憶しておけば、繰り返し処理する場合でも計算量は上記のベクトル計算と場合分けのみで少なく済む
• デメリット
  • 原点の取り方によって19の系に分かれていて、適切な座標系の選択が必要
  • あまり座標原点から離れると誤差が大きくなる（原則東西方向に130km以内）

平面直角座標系での系をわざと間違えてみると、誤差が発生することがわかる。誤差の割合は精度によらずだいたい一定。とはいえ、実用上は高が知れている。

• 6系(大阪付近)の緯度経度を12系(北海道)で平面座標に変換した時の、点と直線の距離の誤差
• 6系の平面直角座標に変換して求めた距離を真値とし、12系をそれと比較

取りうる緯度経度の範囲が広い時、どこであろうと系を気にすることなく一定の精度で距離を計測したい。緯度経度から系を割り出すことはできても、処理内容が北海道の線分と沖縄の点を比較する可能性があるとき、誤差は避けられない。

これは検索してもあまり資料が出てこない。実用上は平面直角座標で何も問題ないし、お役所のお墨付きもあるし、メリットの方が大きいからかもしれない。ただ、理論上は地球を楕円球としてみれば、3次元空間のまま一定の公式で導出できそうなものである。

3次元の地心直交座標系に変換してから3次元ベクトルで同じように内外判定を行ってみた。

• 内外判定は3次元上のベクトル$\vec{AB}$と点Pで行う
• Pの垂線がAB上に無い場合は、近い方の端点との距離を算出する
• Pの垂線がAB上にある場合は、内分比によってAとBの緯度経度の数値をそのまま内分した点$C$を垂線の足とし、$CP$の距離を算出する
• 距離の算出にはHubenyの公式を用いる

3次元上のベクトル$\vec{AB}$と、楕円球表面をなぞった$\stackrel{\frown}{AB}$は異なるため、平面直角座標系と比較して誤差が出る。誤差はABやPが互いに離れているほど大きくなるが、近い場合はほぼ等しくなる。

• 誤差は、平面直角座標で求めた距離が真値であると仮定しての、地心直交座標で求めた距離の誤差
• 精度aとは、ABの緯度差・経度差が各2a度、CPの緯度差・経度差が各a度という条件でテストした結果を示す
• 座標系6系近傍で行った
• 原点に近い条件でテストしたので、比較対象となる平面直角座標の結果は大体正確なはず
  • だが、距離が離れると平面上の距離と地球表面上の距離にも差が出るため、まったく正確とは言えない
  • 結局、真値はわからない
• 結果を見ると、最大でも1km程度の距離を求めたいだけならほぼ誤差20cm未満で収まることがいえる
• だが、距離が離れると誤差は平面直角座標以上に出るわけで、積極的に使う理由があるかというとあまり無いか

平面直角座標と比較して、データ量は2次元→3次元で単純に1.5倍になるが、内外判定の計算量はあまり変わらない。ただ、最終的な距離計算にHubenyを用いるため、その点では遅くなる。より精確になるともいえる。

速度は、ランダムに100万回pythonで計算させると、以下の通り。確かに平面と比較すると遅いが、どちらも十分に速いため大きな問題にはならなさそう。

pythonではpyprojを使うと楽ちんですね

展開

import pyproj

class ConvertPlane:
    def __init__(self, from_epsg=4326, to_epsg=2448):
        self.from_epsg = pyproj.Proj('+init=EPSG:{}'.format(from_epsg))
        self.to_epsg = pyproj.Proj('+init=EPSG:{}'.format(to_epsg))

    def latlng2xy(self, lat, lng):
        return pyproj.transform(self.from_epsg, self.to_epsg, lng, lat)

    def p2l(self, alat, alng, blat, blng, plat, plng):
        a_ = complex(*self.latlng2xy(alat, alng))
        b_ = complex(*self.latlng2xy(blat, blng))
        p_ = complex(*self.latlng2xy(plat, plng))

        ab_, ap_ = b_ - a_, p_ - a_
        abp_prod = ab_.conjugate() * ap_  # 内積 + 外積j
        if abp_prod.real <= 0:
            return abs(ap_)
        ab = abs(ab_)
        ah = abp_prod.real / ab
        if ah >= ab:
            return abs(p_ - b_)
        return abs(abp_prod.imag) / ab


class Convert3d:
    gi = pyproj.Geod(ellps='WGS84').inv

    def __init__(self, from_epsg=4326, to_epsg=5332):
        self.from_epsg = pyproj.Proj('+init=EPSG:{}'.format(from_epsg))
        self.to_epsg = pyproj.Proj('+init=EPSG:{}'.format(to_epsg))

    def latlng2xyz(self, lat, lng):
        return pyproj.transform(self.from_epsg, self.to_epsg, lng, lat)

    def p2l(self, alat, alng, blat, blng, plat, plng):
        a_ = np.array(latlng2xyz(alat, alng))
        b_ = np.array(latlng2xyz(blat, blng))
        p_ = np.array(latlng2xyz(plat, plng))
        ab_, ap_ = b_ - a_, p_ - a_
        dot_prod = np.dot(ab_, ap_)
        if dot_prod <= 0:
            return self.gi(alng, alat, plng, plat)[2]
        ab = np.linalg.norm(ab_)
        ah = dot_prod / ab
        if ah >= ab:
            return self.gi(blng, blat, plng, plat)[2]
        
        # 内分点なら、垂線の足の座標はABの緯度経度を数値的に内分した値とする
        # ここを地球楕円体に沿うように改良するともう少し精度が高まるか？
        rh = ah / ab
        hlat = alat + rh * (blat - alat)
        hlng = alng + rh * (blng - alng)
        return self.gi(plng, plat, hlng, hlat)[2]