回転楕円体面上（地球表面上）の距離などの計算

こんにちは。
回転楕円体面上（地球表面上、WGS84）の距離などの各種計算（近似）を使いやすくまとめたものを見つけたので^{1 2}、少し手直ししてみました（cheap-ruler.js）³。

下記が使用例です（これの二点間距離計算値は、12.4604705 km で、近似誤差は 2.5 mm ⁴）。

usage.js

const lat_ruler = 60.0;
const ruler = cheapRuler(lat_ruler, 'kilometers');
const distance = ruler.distance([135.5, 59.95], [135.6, 60.05]); // #=> 12.4604705 km

ソース（抜粋）は、

cheap-ruler.js

// WGS84 ellipsoidal earth
const consts = cheapRuler.constants = {
    RE: 6378137.0, // IS-GPS
    FE: (1/298.257223563), //IS-GPS
    DEGREE: Math.PI/180
};
consts.E2 = consts.FE * (2 - consts.FE);

function cheapRuler(lat /*: number */, units /*: ?string */) {
    if (lat === undefined) throw new Error('No latitude given.');
    if (units && !factors[units]) throw new Error('Unknown unit ' + units + '. Use one of: ' + Object.keys(factors).join(', '));
    return new CheapRuler(Math.cos(lat * consts.DEGREE), units);
}

// see https://en.wikipedia.org/wiki/Geographical_distance
function CheapRuler(coslat, units) {
    const m = consts.DEGREE * consts.RE / factors[units ? units : 'kilometers'];
    const w2 = 1 / (1 - consts.E2 * (1 - coslat * coslat));
    const w = Math.sqrt(w2);
    this.kx = m * w * coslat;
    this.ky = m * w * w2 * (1 - consts.E2);
}

CheapRuler.prototype = {
    /**
     * Given two points of the form [longitude, latitude], returns the distance.
     *
     * @param {Array<number>} a point [longitude, latitude]
     * @param {Array<number>} b point [longitude, latitude]
     * @returns {number} distance
     * @example
     * const distance = ruler.distance([30.5, 50.5], [30.51, 50.49]);
     * //=distance
     */
    distance: function (a, b) {
        const dx = (a[0] - b[0]) * this.kx;
        const dy = (a[1] - b[1]) * this.ky;
        return Math.hypot(dx, dy);
    }
}

---

1. これは曲率を使った一次近似計算式です（「地理経緯度の変換式」(Wikipedia)）。「回転楕円体面上（地球表面上）の正多角形」と同様の計算です。 ↩

2. polygon領域内部に点が含まれるかどうかの判定については、今回の近似（曲率係数を掛ける）を行っても行わなくても結果には差がありません。したがって「二次元多角形領域内部に点が含まれるかどうかを判定」のような方法をそのまま経緯度値に適用しても近似的には良いでしょう。 ↩

3. 手直し前のオリジナルのcheap-ruler.jsは、$e^2$で展開し係数は数値化した計算式を使っています。 ↩

4. このような距離条件では 「回転楕円体面上（地球表面上）の二点間距離高精度近似」を使う必要はあまりなく、今回の近似方法で十分そうです。 ↩