0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@kkdd

Bowring for short lines（精度良好な地理的距離の短距離近似法）をさらに精度改善

0
Last updated at Posted at 2025-01-04

こんにちは。
精度良好な地理的距離の短距離近似法である Bowring's method for short lines (Wikipedia)1 を、さらに精度改善し、精度評価2を行いました。

これは、"Bowring for short lines" (Karney; geographiclib) をそっくり踏襲したものです（下記）。

精度比較評価

オリジナルの

  • Bowring (1981) （bowring_short）

に比べ、Karneyによる下記の改良式は精度が向上します：

  • 「地理緯度＋中間緯度」を用いる計算（bowring_short_）。
  • さらに「パラメトリック緯度＋中間緯度」を用いた計算の精度はより向上します（bowring_short_beta_；100kmでは10ppt3）。

mini1736353739.jpg

計算速度評価

formula mean calculation time [ms]
flatsurface 23.68
hubeny 38.82
spherical 31.44
straightline 24.79
gaussmidlatitude 41.57
bowring_short 37.55
bowring_short_ 45.41
bowring_short_beta_ 61.89
andoyer 30.54
andoyerlambert 47.42
andoyerlambertthomas 52.63
vincenty 191.15
bowring-inverse 205.84
bowring-inverse_ 204.67
karney 519.06

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
}
bowring_short_.js
// 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
}
bowring_short_beta_.js
// 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
}
straightline.js
function straightline(lat1, lon1, lat2, lon2){  // in degrees
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;

  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const latmid = (lat1 + lat2) / 2 * degree;
  const sinlat = sin(latmid);
  const coslat = cos(latmid);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio	

  return 2 * n * a * hypot(
    sin(londiffhalf) * coslat,
    cos(londiffhalf) * latdiffhalf * m_by_n
  );
}
gaussmidlatitude.js
function gaussmidlatitude(lat1, lon1, lat2, lon2){  // in degrees
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  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 latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const latmid = (lat1 + lat2) / 2 * degree;
  const sinlat = sin(latmid);
  const coslat = cos(latmid);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio	

  // great circle formula
  return n * a * asin2(hypot(
    sin(londiffhalf) * coslat,
    cos(londiffhalf) * sin(latdiffhalf * m_by_n)
  ));
}

  1. Bowring, B. R., "The direct and inverse problems for short geodesic lines on the ellipsoid". Surveying and Mapping (1981), p.135. ↩

  2. 参考： 二点間測地線距離（地球上）：各計算方法の誤差評価#精度評価用入力データ生成 ↩

  3. 誤差は $| \Delta s_\text{approx} - \Delta s_\text{true}| \lessapprox \frac{\Delta s^4}{6000 a^3}$。 ↩

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?