Opened 2 years ago

Closed 2 years ago

## #11917 closed Bugs (fixed)

# geometry default geographic (andoyer) antipodal distance is zero

Reported by: | Owned by: | Barend Gehrels | |
---|---|---|---|

Milestone: | Boost 1.61.0 | Component: | geometry |

Version: | Boost 1.60.0 | Severity: | Problem |

Keywords: | wgs84 distance andoyer | Cc: |

### Description

When calling the default (andoyer) `boost::geometry::distance()`

function,
with `boost::geometry::cs::geographic`

points that are close to being antipodal
(e.g. the North and South poles or points that are diametrically opposed around the equator),
it returns reasonable results. However, when called with antipodal points it returns zero, e.g.:

Andoyer distance near poles: 20003917.164970 near opposite equatorial points: 20037508.151445 Andoyer distance at poles: 0.000000 opposite equatorial points: 0.000000

Here the code I used to obtain the results above:

std::cout.setf(std::ios::fixed); typedef boost::geometry::cs::geographic<boost::geometry::radian> Wgs84Coords; typedef boost::geometry::model::point<double, 2, Wgs84Coords> GeographicPoint; // Min value for delta. 0.000000014 causes Andoyer to fail. const double DELTA(0.000000015); // Note boost points are Long & Lat NOT Lat & Long GeographicPoint near_north_pole (0.0, M_PI_2 - DELTA); GeographicPoint near_south_pole (0.0, -M_PI_2 + DELTA); double andoyer_minor(boost::geometry::distance(north_pole, south_pole)); std::cout << "Andoyer distance near poles: " << andoyer_minor << "\n"; GeographicPoint near_equator_east ( M_PI_2 - DELTA, 0.0); GeographicPoint near_equator_west (-M_PI_2 + DELTA, 0.0); double andoyer_major(boost::geometry::distance(near_equator_east, near_equator_west)); std::cout << "near opposite equatorial points: " << andoyer_major << "\n\n"; GeographicPoint north_pole (0.0, M_PI_2); GeographicPoint south_pole (0.0, -M_PI_2); andoyer_minor = boost::geometry::distance(north_pole, south_pole); std::cout << "Andoyer distance at poles: " << andoyer_minor << "\n"; GeographicPoint equator_east ( M_PI_2, 0.0); GeographicPoint equator_west (-M_PI_2, 0.0); andoyer_major = boost::geometry::distance(equator_east, equator_west); std::cout << "opposite equatorial points: " << andoyer_major << std::endl;

This is clearly a bug. The distance between the antipodal points should be slightly larger than the distance between the other points, and definitely NOT zero!

### Attachments (0)

### Change History (5)

### comment:1 Changed 2 years ago by

### comment:2 Changed 2 years ago by

I agree with Charles Karney's comment above.

I originally raised this issue after finding an issue with the `Vincenty`

function, see: http://stackoverflow.com/questions/34767736/why-is-boostgeometry-geographic-vincenty-distance-inaccurate-around-the-equato

On further investigation (and with a bit of help from Charles) I came to realise that both the `Vincenty`

and `Andoyer`

algorithms calculate incorrect distances for antipodal points on the equator. Whilst [GeographicLib?](http://geographiclib.sourceforge.net/) calculates the correct distance.

It should be included in boost geometry and replace `Andoyer`

as the default algorithm.

BTW, please note that Charles's link to `GeographicLib`

is incorrect, please use the one in this comment.

### comment:3 Changed 2 years ago by

It might also be worth pointing out that "Algorithms for Geodesics" also provides a method for computing the area of a geodesic polygon (a polygon whose sides are geodesics) on an ellipsoid.

### comment:4 Changed 2 years ago by

Thanks, both, and we welcome the papers and involvement of Charles Karney! We were already aware of your library and your papers, also with respect to area, this is very helpful.

### comment:5 Changed 2 years ago by

Milestone: | To Be Determined → Boost 1.61.0 |
---|---|

Resolution: | → fixed |

Status: | new → closed |

Version: | Boost 1.61.0 → Boost 1.60.0 |

**Note:**See TracTickets for help on using tickets.

I recommend that boost::geometry include the algorithm for geodesic distance given in my paper, Algorithms for geodesics. This returns an accurate result for the distance for arbitrary points on terrestrially relevant ellipsoids. GeographicLib (version 1.25 and later) also includes the solution in terms of elliptic integrals which is accurate for ellipsoids of revolution with 1/100 < b/a < 100.