こんにちは。
精度良好な地理的距離の短距離近似法である Bowring's method for short lines (Wikipedia)1 を、さらに精度改善しました。
用いた計算方法は、"Bowring for short lines"(Karney; geographiclib) をそっくり踏襲し、中間緯度(mid-latitude)を用いて精度改善しました。
精度評価2を行うと、
-
「パラメトリック緯度+中間緯度」を用いた計算は、短距離近似法としては確かに最も精度良好のようで(下記のプロット内の
bowring_short_beta_
;100kmでは10ppt)、 -
「地理緯度+中間緯度」を用いる計算はこれに次いで精度良好のようです(
bowring_short_
)。
精度評価
source code
bowring_short.js
// Bowring, B. R. (1981). "The direct and inverse problems for short geodesic lines on the ellipsoid". Surveying and Mapping. 41 (2): 135–141.
function bowring_short(lat1, lon1, lat2, lon2){ // in degrees
"use strict";
const a = 6378137.0; // GRS80
const f = 1 / 298.257223563; // WGS84
const e2 = f * (2 - f); // first eccentricity squared
const ep2 = e2 / (1 - f) ** 2; // second eccentricity squared; epsilon
const degree = Math.PI / 180.0;
const sin = Math.sin;
const cos = Math.cos;
const sqrt = Math.sqrt;
const hypot = Math.hypot;
const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;
const londiffhalf = (lon2 - lon1) / 2 * degree;
lat1 = lat1 * degree;
lat2 = lat2 * degree;
const latdiffhalf = (lat2 - lat1) / 2;
const lat0 = lat1; // Bowring
const cos_lat0 = cos(lat0);
const cos_lat0_squared = cos_lat0 ** 2;
const C = sqrt(1 + ep2);
const B2 = 1 + ep2 * cos_lat0_squared;
const B = sqrt(B2);
const A = sqrt(1 + ep2 * cos_lat0_squared ** 2);
// radius, longitude, and latitude on the conformal sphere
const Rp = C / B2;
const lonpdiffhalf = A * londiffhalf;
const latpdiffhalf = latdiffhalf / B * (1 + (3 / 2 * ep2 / B2) * latdiffhalf * sin(2 * lat0 + 4 / 3 * latdiffhalf)); // D
const cos_latp0 = (B * cos_lat0 * cos(latpdiffhalf) - sin(lat0) * sin(latpdiffhalf)) / A;
// great circle formula
const E = cos(lonpdiffhalf) * sin(latpdiffhalf);
const F = sin(lonpdiffhalf) * cos_latp0;
return a * Rp * asin2(hypot(E, F)); // geodesic distance
}
// by Karney
// https://geographiclib.sourceforge.io/C++/doc/geodesic.html#bowring
function bowring_short_(lat1, lon1, lat2, lon2){ // in degrees
"use strict";
const a = 6378137.0; // GRS80
const f = 1 / 298.257223563; // WGS84
const e2 = f * (2 - f); // first eccentricity squared
const ep2 = e2 / (1 - f) ** 2; // second eccentricity squared; epsilon
const degree = Math.PI / 180.0;
const sin = Math.sin;
const cos = Math.cos;
const tan = Math.tan;
const atan = Math.atan;
const sqrt = Math.sqrt;
const hypot = Math.hypot;
const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;
const londiffhalf = (lon2 - lon1) / 2 * degree;
lat1 = lat1 * degree, lat2 = lat2 * degree;
const latm = (lat1 + lat2) / 2; // mid-latitude by Karney
const cosm2 = cos(latm) ** 2;
const C = sqrt(1 + ep2);
const B2 = 1 + ep2 * cosm2;
const B = sqrt(B2);
const A = sqrt(1 + ep2 * cosm2 ** 2);
// radius, longitude, and latitude on the conformal sphere
const Rp = C / B2;
const lonpdiffhalf = A * londiffhalf;
const dlat = lat1 - latm; // = -(lat2 - lat)
const phipm = atan(tan(latm)/B) + dlat*dlat * 3*ep2/(4*B*B2) * sin(2*latm) * cos(2/3*dlat);
const phipdiffhalf = dlat/B * (1 + (3*ep2*dlat)/(4*B2) * cos(2*latm) * sin(2/3*dlat));
// great circle formula
const E = cos(lonpdiffhalf) * sin(phipdiffhalf);
const F = sin(lonpdiffhalf) * cos(phipm);
return a * Rp * asin2(hypot(E, F)); // geodesic distance
}
// by Karney
// https://geographiclib.sourceforge.io/C++/doc/geodesic.html#bowring
function bowring_short_beta_(lat1, lon1, lat2, lon2){ // in degrees
"use strict";
const a = 6378137.0; // GRS80
const f = 1 / 298.257223563; // WGS84
const f1 = 1 - f;
const e2 = f * (2 - f); // first eccentricity squared
const ep2 = e2 / (f1*f1); // second eccentricity squared; epsilon
const degree = Math.PI / 180.0;
const sin = Math.sin;
const cos = Math.cos;
const tan = Math.tan;
const atan = Math.atan;
const sqrt = Math.sqrt;
const hypot = Math.hypot;
const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;
const londiffhalf = (lon1 - lon2) / 2 * degree;
const bet1 = atan(tan(lat1 * degree) * f1); // reduced latitude
const bet2 = atan(tan(lat2 * degree) * f1);
const betm = (bet1 + bet2)/2; // mid-latitude by Karney
const sbetm = sin(betm);
const cbetm = cos(betm);
const phim = atan(tan(betm) / f1);
const cosm2 = cos(phim) ** 2;
const C = sqrt(1 + ep2);
const B2 = 1 + ep2 * cosm2;
const B = sqrt(B2);
const A = sqrt(1 + ep2 * cosm2 ** 2);
// radius, longitude, and latitude on the conformal sphere
const Rp = C / B2;
const lonpdiffhalf = A * londiffhalf;
const dnm2 = 1 + ep2*sbetm*sbetm;
const dnm = sqrt(dnm2);
const dbet = (bet1-betm)/dnm; // = -(bet2-betm)/dnm
const phipm = atan(tan(phim)/B) + dbet * dbet * ep2/2 * cbetm*sbetm/dnm;
const phipdiffhalf = dbet * (1 + dbet * dbet * ep2/6 * (cbetm*cbetm/dnm2 - sbetm*sbetm) );
// great circle formula
const E = cos(lonpdiffhalf) * sin(phipdiffhalf);
const F = sin(lonpdiffhalf) * cos(phipm);
return a * Rp * asin2(hypot(E, F)); // geodesic distance
}