boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-01 08:32:10 +00:00
Code Issues Projects Releases Wiki Activity

[formula] Fix errors in inverse formulas (manifesting near poles).

Browse Source 
vincenty - fix error in formula (missing sqr)
differential_quantities - fix error in formula (wrong equation and lack of normalization)
andoyer - wrong azimuth at south pole
This commit is contained in:
Adam Wulkiewicz
2017-03-02 16:54:36 +01:00
parent e55968fa7b
commit bdd2d2c60c
3 changed files with 48 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
include/boost/geometry/formulas/andoyer_inverse.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2015-2016 Oracle and/or its affiliates.
// Copyright (c) 2015-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -137,7 +137,14 @@ public:
CT A = c0;
CT U = c0;
if ( ! math::equals(cos_lat2, c0) )
if (math::equals(cos_lat2, c0))
{
if (sin_lat2 < c0)
{
A = pi;
}
}
else
{
CT const tan_lat2 = sin_lat2/cos_lat2;
CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon;
@@ -148,7 +155,14 @@ public:
CT B = c0;
CT V = c0;
if ( ! math::equals(cos_lat1, c0) )
if (math::equals(cos_lat1, c0))
{
if (sin_lat1 < c0)
{
B = pi;
}
}
else
{
CT const tan_lat1 = sin_lat1/cos_lat1;
CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon;

19
include/boost/geometry/formulas/differential_quantities.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2016 Oracle and/or its affiliates.
// Copyright (c) 2016-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -64,8 +64,8 @@ public:
CT const c1 = 1;
CT const one_minus_f = c1 - f;
CT const sin_bet1 = one_minus_f * sin_lat1;
CT const sin_bet2 = one_minus_f * sin_lat2;
CT sin_bet1 = one_minus_f * sin_lat1;
CT sin_bet2 = one_minus_f * sin_lat2;
// equator
if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0))
@@ -89,14 +89,17 @@ public:
CT const e2 = f * (c2 - f);
CT const ep2 = e2 / math::sqr(one_minus_f);
CT const cos_bet1 = cos_lat1;
CT const cos_bet2 = cos_lat2;
CT const sin_alp1 = sin(azimuth);
CT const cos_alp1 = cos(azimuth);
//CT const sin_alp2 = sin(reverse_azimuth);
CT const cos_alp2 = cos(reverse_azimuth);
CT cos_bet1 = cos_lat1;
CT cos_bet2 = cos_lat2;
normalize(sin_bet1, cos_bet1);
normalize(sin_bet2, cos_bet2);
CT sin_sig1 = sin_bet1;
CT cos_sig1 = cos_alp1 * cos_bet1;
CT sin_sig2 = sin_bet2;
@@ -112,8 +115,8 @@ public:
J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) :
J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ;
CT const dn1 = math::sqrt(c1 + e2 * math::sqr(sin_lat1));
CT const dn2 = math::sqrt(c1 + e2 * math::sqr(sin_lat2));
CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1));
CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2));
if (BOOST_GEOMETRY_CONDITION(EnableReducedLength))
{

35
include/boost/geometry/formulas/vincenty_inverse.hpp
View File

@@ -2,8 +2,8 @@
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// This file was modified by Oracle on 2014, 2016.
// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates.
// This file was modified by Oracle on 2014, 2016, 2017.
// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -40,7 +40,7 @@ namespace boost { namespace geometry { namespace formula
\brief The solution of the inverse problem of geodesics on latlong coordinates, after Vincenty, 1975
\author See
- http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
- http://www.icsm.gov.au/gda/gdav2.3.pdf
- http://www.icsm.gov.au/gda/gda-v_2.4.pdf
\author Adapted from various implementations to get it close to the original document
- http://www.movable-type.co.uk/scripts/LatLongVincenty.html
- http://exogen.case.edu/projects/geopy/source/geopy.distance.html
@@ -98,10 +98,10 @@ public:
CT const radius_a = CT(get_radius<0>(spheroid));
CT const radius_b = CT(get_radius<2>(spheroid));
CT const flattening = formula::flattening<CT>(spheroid);
CT const f = formula::flattening<CT>(spheroid);
// U: reduced latitude, defined by tan U = (1-f) tan phi
CT const one_min_f = c1 - flattening;
CT const one_min_f = c1 - f;
CT const tan_U1 = one_min_f * tan(lat1); // above (1)
CT const tan_U2 = one_min_f * tan(lat2); // above (1)
@@ -112,8 +112,9 @@ public:
CT const cos_U1 = c1 / temp_den_U1;
CT const cos_U2 = c1 / temp_den_U2;
// sin = tan / sqrt(1 + tan^2)
CT const sin_U1 = tan_U1 / temp_den_U1;
CT const sin_U2 = tan_U2 / temp_den_U2;
// sin = tan * cos
CT const sin_U1 = tan_U1 * cos_U1;
CT const sin_U2 = tan_U2 * cos_U2;
// calculate sin U and cos U directly
//CT const U1 = atan(tan_U1);
@@ -129,7 +130,8 @@ public:
CT sin_sigma;
CT sin_alpha;
CT cos2_alpha;
CT cos2_sigma_m;
CT cos_2sigma_m;
CT cos2_2sigma_m;
CT sigma;
int counter = 0; // robustness
@@ -143,12 +145,13 @@ public:
CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
cos2_alpha = c1 - math::sqr(sin_alpha);
cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
cos_2sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
cos2_2sigma_m = math::sqr(cos_2sigma_m);
CT C = flattening/c16 * cos2_alpha * (c4 + flattening * (c4 - c3 * cos2_alpha)); // (10)
CT C = f/c16 * cos2_alpha * (c4 + f * (c4 - c3 * cos2_alpha)); // (10)
sigma = atan2(sin_sigma, cos_sigma); // (16)
lambda = L + (c1 - C) * flattening * sin_alpha *
(sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11)
lambda = L + (c1 - C) * f * sin_alpha *
(sigma + C * sin_sigma * (cos_2sigma_m + C * cos_sigma * (-c1 + c2 * cos2_2sigma_m))); // (11)
++counter; // robustness
@@ -181,8 +184,10 @@ public:
CT A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
CT B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
CT delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
- (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6)
CT const cos_sigma = cos(sigma);
CT const sin2_sigma = math::sqr(sin_sigma);
CT delta_sigma = B * sin_sigma * (cos_2sigma_m + (B/c4) * (cos_sigma* (-c1 + c2 * cos2_2sigma_m)
- (B/c6) * cos_2sigma_m * (-c3 + c4 * sin2_sigma) * (-c3 + c4 * cos2_2sigma_m))); // (6)
result.distance = radius_b * A * (sigma - delta_sigma); // (19)
}
@@ -205,7 +210,7 @@ public:
typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
quantities::apply(lon1, lat1, lon2, lat2,
result.azimuth, result.reverse_azimuth,
radius_b, flattening,
radius_b, f,
result.reduced_length, result.geodesic_scale);
}