はじめに

地理座標（経度,緯度）で与えられた2点間の距離や、ポリゴンの面積は、最近はライブラリ^[1]が充実しているので簡単に求められるようになった。

しかし、ライブラリを使うにしても、準拠楕円体は何を採用しているか、どのようなアルゴリズムを採用しているかを知らないと、いざというときに恥をかいてしまうかも。

距離や面積をより正確に計算する方法を探す中で調べた諸々をまとめておきます。

作ったパッケージ:geo4326

準拠楕円体

計算時に使用する半径などの値は、地球をどの楕円体に近似^[2]するかで決まる。

	長半径a [m]	扁平率f	平均半径 mean radius [m]	備考
WGS84	6378137	1 / 298.257223563	6371008.771	GPSの測地系に使用されている
GRS80	6378137	1 / 298.257222101	6371008.771	世界測地系（日本測地系2011）に採用されている
ベッセル	6377397.155	1 / 299.152813	6370291.091	古い日本測地系で採用されていた

極半径 b = a(1 - f)
平均半径 r = (2a - b)/3

地球を球体として扱う

公式が使えて計算が楽。
半径には平均半径を使う。

距離

球体上の最短距離である大円距離(大圏距離)を求める。
英語で調べるときはgreat circleやgeodesic lineあたりでググる。
公式はwikipediaにある通り。
ライブラリはだいたいこの公式を採用している。

面積

指定した座標が囲む球体の表面の面積を求める。

球上の面積を求めるリュイリエの公式(l’Huiller’s formula)は複雑でやってられないから簡易式を使おう、というNASAの論文のアルゴリズムを使っているケースが多い。

Openlayers
Leaflet(プラグイン) Leaflet.MeasureAreaControl
Turf.js

地球を楕円体として扱う

より正確。
球と比べて計算がまあめんどくさい。

距離

Hubenyの式

参考: http://tancro.e-central.tv/grandmaster/excel/hubenystandard.html
比較的軽量な計算で楕円体上の最短距離を計算できる。

Vincenty法

参考: https://www.movable-type.co.uk/scripts/latlong-vincenty.html
Hubenyの式よりも精度が高い（誤差0.5mm以下まで攻めることができるらしい）。
ただし2点のとり方によっては計算が収束しないことがある。

国土地理院の方法

参考: https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/bl2st/bl2st.htm

距離と方位角の計算ページや地理院地図電子国土Webで使われているアルゴリズム。
Vincenty法と比べてどういったメリットがあるのかは、元の論文を読んだわけではないのでわかっていない。しかし国土地理院がわざわざ採用しているのだからきっと良いに違いない。

JavaScriptで愚直に実装してみた^[3]。変数名やフォーマットは各々いい感じにしてほしい。

EPSILON = 1e-8
SEMEMAJOR_AXIS_GRS80 = 6378137
FLATTENING_GRS80 = 298.257222101

function _eq(a, b) {
  return Math.abs(a - b) < EPSILON;
}

function _toRadians(degree) {
  return degree * (Math.PI / 180);
}

export function distance(p1, p2, userOptions = {}) {
  // https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/bl2st/bl2st.htm
  // https://www.tandfonline.com/doi/abs/10.1179/sre.1996.33.261.461
  const options = Object.assign(
    {
      semimajorAxis: SEMEMAJOR_AXIS_GRS80,
      flattening: FLATTENING_GRS80,
      truncation: 1e-14,
      maxCount: 100,
    },
    userOptions
  );
  const lon1 = _toRadians(p1[0]);
  const lat1 = _toRadians(p1[1]);
  const lon2 = _toRadians(p2[0]);
  const lat2 = _toRadians(p2[1]);
  const a = options.semimajorAxis;
  const f = 1 / options.flattening;
  const lonDiff = lon2 - lon1;
  let lonDiffWarp = lonDiff;
  if (lonDiff > Math.PI) lonDiffWarp -= 2 * Math.PI;
  else if (lonDiff < -Math.PI) lonDiffWarp += 2 * Math.PI;

  const L = Math.abs(lonDiffWarp);
  const L_d = Math.PI - L;
  const Sigma = lat1 + lat2;
  const u1 =
    lonDiffWarp >= 0
      ? Math.atan((1 - f) * Math.tan(lat1))
      : Math.atan((1 - f) * Math.tan(lat2));
  const u2 =
    lonDiffWarp >= 0
      ? Math.atan((1 - f) * Math.tan(lat2))
      : Math.atan((1 - f) * Math.tan(lat1));
  const Sigma_d = u1 + u2;
  const Delta_d = u2 - u1;
  const xi = Math.cos(Sigma_d / 2);
  const xi_d = Math.sin(Sigma_d / 2);
  const eta = Math.sin(Delta_d / 2);
  const eta_d = Math.cos(Delta_d / 2);
  const x = Math.sin(u1) * Math.sin(u2);
  const y = Math.cos(u1) * Math.cos(u2);
  const c = y * Math.cos(L) + x;
  const epsilon = (f * (2 - f)) / (1 - f) ** 2;
  let theta, zone;
  if (c >= 0) {
    theta = L * (1 + f * y);
    zone = 1;
  } else if (c >= -Math.cos(_toRadians(3) * Math.cos(u1))) {
    theta = L_d;
    zone = 2;
  } else {
    zone = 3;
    const R =
      f *
      Math.PI *
      Math.cos(u1) ** 2 *
      (1 -
        (f * (1 + f) * Math.sin(u1) ** 2) / 4 +
        (3 / 16) * f ** 2 * Math.sin(u1) ** 4);
    const d1 = L_d * Math.cos(u1) - R;
    const d2 = Math.abs(Sigma_d) + R;
    const q = L_d / (f * Math.PI);
    const f1 = (f * (1 + f / 2)) / 4;
    const gamma0 = q + f1 * q - f1 * q ** 3;

    if (_eq(Sigma, 0)) {
      if (d1 > 0) {
        //b1
        theta = L_d;
      } else if(_eq(d1, 0)) {
        // b2
        const Gamma = Math.sin(u1) ** 2;
        const n0 =
          (epsilon * Gamma) / (Math.sqrt(1 + epsilon * Gamma) + 1) ** 2;
        const A = (1 + n0) * (1 + (5 / 4) * n0 ** 2);
        return (1 - f) * a * A * Math.PI;
      }else {
        let cnt = 0;
        let gamma0 = 0;
        let Gamma = 0;
        while (cnt < options.maxCount) {
          Gamma = 1 - gamma0 ** 2;
          const D = f * (1 + f) / 4 - 3 / 16 * f ** 2 * Gamma;
          const gamma1 = q * (1 - D * Gamma);
          if (Math.abs(gamma0 - gamma1) < options.truncation) break;
          gamma0 = gamma1;
          cnt += 1;
        }
        if (cnt >= options.maxCount) throw new Error();
        const n0 =
          (epsilon * Gamma) / (Math.sqrt(1 + epsilon * Gamma) + 1) ** 2;
        const A = (1 + n0) * (1 + (5 / 4) * n0 ** 2);
        return (1 - f) * a * A * Math.PI;      
      }
    } else {
      // a
      const A0 = Math.atan(d1 / d2);
      const B0 = Math.asin(R / Math.sqrt(d1 ** 2 + d2 ** 2));
      const psi = A0 + B0;
      const j = gamma0 / Math.cos(u1);
      const k =
        ((1 + f1) * Math.abs(Sigma_d) * (1 - f * y)) / (f * Math.PI * y);
      const j1 = j / (1 + k * Math.sin(psi));
      const psi_d = Math.asin(j1);
      const psi_dd = Math.asin((j1 * Math.cos(u1)) / Math.cos(u2));
      theta =
        2 *
        Math.atan(
          (Math.tan((psi_d + psi_dd) / 2) * Math.sin(Math.abs(Sigma_d) / 2)) /
            Math.cos(Delta_d / 2)
        );
    }
  }

  let cnt = 0;
  let Gamma = 0,
    sigma = 0,
    zeta = 0,
    J = 0,
    K = 0;
  while (cnt < options.maxCount) {
    const g =
      zone === 1
        ? Math.sqrt(
            eta ** 2 * Math.cos(theta / 2) ** 2 +
              xi ** 2 * Math.sin(theta / 2) ** 2
          )
        : Math.sqrt(
            eta ** 2 * Math.sin(theta / 2) ** 2 +
              xi ** 2 * Math.cos(theta / 2) ** 2
          );
    const h =
      zone === 1
        ? Math.sqrt(
            eta_d ** 2 * Math.cos(theta / 2) ** 2 +
              xi_d ** 2 * Math.sin(theta / 2) ** 2
          )
        : Math.sqrt(
            eta_d ** 2 * Math.sin(theta / 2) ** 2 +
              xi_d ** 2 * Math.cos(theta / 2) ** 2
          );
    sigma = 2 * Math.atan(g / h);
    J = 2 * g * h;
    K = h ** 2 - g ** 2;
    const gamma = (y * Math.sin(theta)) / J;
    Gamma = 1 - gamma ** 2;
    zeta = Gamma * K - 2 * x;
    const zeta_d = zeta + x;
    const D = (f * (1 + f)) / 4 - (3 / 16) * f ** 2 * Gamma;
    const E =
      (1 - D * Gamma) *
      f *
      gamma *
      (sigma + D * J * (zeta + D * K * (2 * zeta ** 2 - Gamma ** 2)));
    const F = zone === 1 ? theta - L - E : theta - L_d + E;
    const G =
      f * gamma ** 2 * (1 - 2 * D * Gamma) +
      ((f * zeta_d * sigma) / J) * (1 - D * Gamma + (f * gamma ** 2) / 2) +
      (f ** 2 * zeta * zeta_d) / 4;
    theta -= F / (1 - G);
    if (Math.abs(F) < options.truncation) break;
    cnt += 1;
  }
  if (cnt >= options.maxCount) throw new Error();

  const n0 = (epsilon * Gamma) / (Math.sqrt(1 + epsilon * Gamma) + 1) ** 2;
  const A = (1 + n0) * (1 + (5 / 4) * n0 ** 2);
  const B =
    (epsilon * (1 - (3 * n0 ** 2) / 8)) /
    (Math.sqrt(1 + epsilon * Gamma) + 1) ** 2;
  return (
    (1 - f) *
    a *
    A *
    (sigma -
      B *
        J *
        (zeta -
          (B *
            (K * (Gamma ** 2 - 2 * zeta ** 2) -
              (B * zeta * (1 - 4 * K ** 2) * (3 * Gamma ** 2 - 4 * zeta ** 2)) /
                6)) /
            4))
  );
}

面積

国土地理院の方法

参考: https://maps.gsi.go.jp/help/pdf/calc_area.pdf
楕円体を球に正積投影し、球上の面積の計算も簡易式でなくリュイリエの公式をきちんと計算する。
地理院地図電子国土Webで採用されている。

こちらもJavaScriptでの実装例を貼っておく。

const _latitudeBelt = (
  latRad,
  e
) => {
  return (
    (e * Math.sin(latRad)) / (1 - (e * Math.sin(latRad)) ** 2) +
    Math.atanh(e * Math.sin(latRad))
  );
};
const _authalicLatitude = (latRad: number, e) => {
  return Math.asin(_latitudeBelt(latRad, e) / _latitudeBelt(Math.PI / 2, e));
};

const _haversine = (
  lonRad1,
  latRad1,
  lonRad2,
  latRad2
): number => {
  return (
    2 *
    Math.asin(
      Math.sqrt(
        Math.abs(
          Math.sin((latRad1 - latRad2) / 2) ** 2 +
            Math.cos(latRad1) *
              Math.cos(latRad2) *
              Math.sin((lonRad1 - lonRad2) / 2) ** 2
        )
      )
    )
  );
};

export function area(linearRing, userOptions = {}) {
  // https://maps.gsi.go.jp/help/pdf/calc_area.pdf
  const options = Object.assign(
    {
      semimajorAxis: SEMEMAJOR_AXIS_GRS80,
      flattening: FLATTENING_GRS80,
    },
    userOptions
  );
  const f = 1 / options.flattening;
  const e2 = f * (2 - f);
  const e = Math.sqrt(e2);
  const r2 =
    (_latitudeBelt(Math.PI / 2, e) * options.semimajorAxis ** 2 * (1 / e - e)) /
    2;

  let area = 0;
  const lon0 = _toRadians(linearRing[0][0]);
  const lat0 = _authalicLatitude(_toRadians(linearRing[0][1]), e);
  let lon1 = _toRadians(linearRing[1][0]);
  let lat1 = _authalicLatitude(_toRadians(linearRing[1][1]), e);

  const l = linearRing.length - 2;
  for (let i = 1; i < l; i++) {
    const lon2 = _toRadians(linearRing[i + 1][0]);
    const lat2 = _authalicLatitude(_toRadians(linearRing[i + 1][1]), e);

    const prime =
      Math.cos(lat0) *
        Math.cos(lon0) *
        Math.cos(lat1) *
        Math.sin(lon1) *
        Math.sin(lat2) +
        Math.cos(lat1) *
          Math.cos(lon1) *
          Math.cos(lat2) *
          Math.sin(lon2) *
          Math.sin(lat0) +
        Math.cos(lat2) *
          Math.cos(lon2) *
          Math.cos(lat0) *
          Math.sin(lon0) *
          Math.sin(lat1) -
        Math.cos(lat2) *
          Math.cos(lon2) *
          Math.cos(lat1) *
          Math.sin(lon1) *
          Math.sin(lat0) -
        Math.cos(lat1) *
          Math.cos(lon1) *
          Math.cos(lat0) *
          Math.sin(lon0) *
          Math.sin(lat2) -
        Math.cos(lat0) *
          Math.cos(lon0) *
          Math.cos(lat2) *
          Math.sin(lon2) *
          Math.sin(lat1) >
      0
        ? 1
        : -1;
    const d1 = _haversine(lon0, lat0, lon1, lat1);
    const d2 = _haversine(lon1, lat1, lon2, lat2);
    const d3 = _haversine(lon2, lat2, lon0, lat0);
    const s = (d1 + d2 + d3) / 2;
    area +=
      prime *
      Math.atan(
        Math.sqrt(
          Math.abs(
            Math.tan(s / 2) *
              Math.tan((s - d1) / 2) *
              Math.tan((s - d2) / 2) *
              Math.tan((s - d3) / 2)
          )
        )
      );
    lon1 = lon2;
    lat1 = lat2;
  }
  return Math.abs(4 * r2 * area);
}

余談すごいぞ地理院地図電子国土Web

上記の通り、地理院地図電子国土Webは距離も面積も楕円体で計算を行ってくれる（電子国土WebはLeafletを利用しているが、計算するようのプラグインを自前で開発している）。

またポリゴンを地図上に描画するのにも大圏航路(楕円体上の大円の一部として2点間の線を描画)と等角航路(横軸経度、縦軸緯度の平面上の直線として地図上に描画)を選択することができる。

大圏航路の表示にはLeafletの既存のプラグインを利用しているが、わざわざ表示の選択ができるようにしてある。
測地のプロとしての気概が感じられ、調べている中で大変好感だった。

なお、OpenlayersやMapboxでは基本的に等角航路での表示となる。
軽く調べた範囲では大圏航路で表示する機能やプラグインは見つけられなかった。
これらのライブラリで大圏航路を表示したいときはarc.jsを使って点の補間をしてやるとよいだろう。

まとめ

主なライブラリはどれも地球を球に近似して距離や面積を計算している。
楕円体近似で計算することもできるので、コードを参考にやってみよう。
地理院地図は匠のこだわりが詰め込まれている。

脚注

ここでは主にJavaScriptのライブラリを対象とする。 ↩︎
そもそも近似でしか計算できない。地球は本当は凸凹のじゃがいもみたいな形であり、真の距離や面積が知りたければ実測する他無い。 ↩︎
方位角の計算は省いている。 ↩︎

Discussion