Open mbecker opened 5 years ago
Hi @mbecker
First of all thank you for filling this issue and spending your time researching the problem.
I agree with all your points.
Do you maybe have the time to implement them yourself. If not, can I use the DistanceVincenty()
code you provided? I think all except 2. are easy to implement. For 2, I'd probably create a new DistanceWithOptions(point, point, DistanceOpts)
function, and the DistanceOpts
would contain the distance algorithm type, the number of iterations (for Vincenti's,) etc. I'd leave the existing functions for backward compatibility, but internally they would call this new function (with options so that the result is still the same).
For example:
type DistanceAlgorithm func(Point, Point) (float64, error)
type DistanceOpts struct {
Algorithm DistanceAlgorithm
IgnoreElevations bool
Iterations int
}
What do you think?
Ps. Regarding:
I still don't kow "how strava is calculating the distance"
Neither do I. I tried to understand it (not only for Strava) times, but no luck. :( At the end of the day every software uses its own heuristics to make a good guesstimate.
Hi @tkrajina, thanks or the feedback! Would be glad to help. Unfortunately I'm quite new to go and do not understand how the "DistanceOpts" and "DistanceAlgorithm" should be used by the user. Any hint would be nice!
I've implemented it in a more easy way. The requirement should be at first not to break your API and that's why I've just added a new function: func (g *GPX) LengthVincenty() float64 (I'm igoring the error and would return 0.0)
gpx.go
// LengthVincenty returns the Vincenty length of all tracks
func (g *GPX) LengthVincenty() (float64, error) {
var lengthVincenty float64
for _, trk := range g.Tracks {
length, err := trk.LengthVincenty()
if err != nil {
return 0, err
}
lengthVincenty += length
}
return lengthVincenty, nil
}
// LengthVincenty returns the Vincenty length of a GPX track.
func (trk *GPXTrack) LengthVincenty() (float64, error) {
var lengthVincenty float64
for _, seg := range trk.Segments {
length, err := seg.LengthVincenty()
if err != nil {
return 0, err
}
lengthVincenty += length
}
return lengthVincenty, nil
}
// LengthVincenty returns the Vincenty length of a GPX segment.
func (seg *GPXTrackSegment) LengthVincenty() (float64, error) {
points := make([]Point, len(seg.Points))
for pointNo, point := range seg.Points {
points[pointNo] = point.Point
}
lengthVincenty, err := LengthVincenty(points)
// If the Vincenty formula can not calculate the distance between two points then return "0.0" (inseatd of an error)
if err != nil {
return 0.0, err
}
return lengthVincenty, nil
}
geo.go
const (
oneDegree = 1000.0 * 10000.8 / 90.0
earthRadius float64 = 6378137 // WGS-84 ellipsoid; See https://en.wikipedia.org/wiki/World_Geodetic_System
flattening float64 = 1 / 298.257223563
semiMinorAxisB float64 = 6356752.314245
//
epsilon = 1e-12
maxIterations = 200
)
// Vincenty formula
func toRadians(deg float64) float64 {
return deg * (math.Pi / 180)
}
func lengthVincenty(locs []Point) (float64, error) {
var previousLoc Point
var res float64
for k, v := range locs {
if k > 0 {
previousLoc = locs[k-1]
d, err := v.DistanceVincenty(&previousLoc)
if err != nil {
return 0, err
}
res += d
}
}
return res, nil
}
// LengthVincenty returns the geographical distance in km between the points p1 (lat1, long1) and p2 (lat2, long2) using Vincenty's inverse formula.
// The surface of the Earth is approximated by the WGS-84 ellipsoid.
// This method may fail to converge for nearly antipodal points.
// https://github.com/asmarques/geodist/blob/master/vincenty.go
func DistanceVincenty(lat1, long1, lat2, long2 float64) (float64, error) {
if lat1 == lat2 && long1 == long2 {
return 0, nil
}
U1 := math.Atan((1 - flattening) * math.Tan(toRadians(lat1)))
U2 := math.Atan((1 - flattening) * math.Tan(toRadians(lat2)))
L := toRadians(long2 - long1)
sinU1 := math.Sin(U1)
cosU1 := math.Cos(U1)
sinU2 := math.Sin(U2)
cosU2 := math.Cos(U2)
lambda := L
result := math.NaN()
for i := 0; i < maxIterations; i++ {
curLambda := lambda
sinSigma := math.Sqrt(math.Pow(cosU2*math.Sin(lambda), 2) +
math.Pow(cosU1*sinU2-sinU1*cosU2*math.Cos(lambda), 2))
cosSigma := sinU1*sinU2 + cosU1*cosU2*math.Cos(lambda)
sigma := math.Atan2(sinSigma, cosSigma)
sinAlpha := (cosU1 * cosU2 * math.Sin(lambda)) / math.Sin(sigma)
cosSqrAlpha := 1 - math.Pow(sinAlpha, 2)
cos2sigmam := 0.0
if cosSqrAlpha != 0 {
cos2sigmam = math.Cos(sigma) - ((2 * sinU1 * sinU2) / cosSqrAlpha)
}
C := (flattening / 16) * cosSqrAlpha * (4 + flattening*(4-3*cosSqrAlpha))
lambda = L + (1-C)*flattening*sinAlpha*(sigma+C*sinSigma*(cos2sigmam+C*cosSigma*(-1+2*math.Pow(cos2sigmam, 2))))
if math.Abs(lambda-curLambda) < epsilon {
uSqr := cosSqrAlpha * ((math.Pow(earthRadius, 2) - math.Pow(semiMinorAxisB, 2)) / math.Pow(semiMinorAxisB, 2))
k1 := (math.Sqrt(1+uSqr) - 1) / (math.Sqrt(1+uSqr) + 1)
A := (1 + (math.Pow(k1, 2) / 4)) / (1 - k1)
B := k1 * (1 - (3*math.Pow(k1, 2))/8)
deltaSigma := B * sinSigma * (cos2sigmam + (B/4)*(cosSigma*(-1+2*math.Pow(cos2sigmam, 2))-
(B/6)*cos2sigmam*(-3+4*math.Pow(sinSigma, 2))*(-3+4*math.Pow(cos2sigmam, 2))))
s := semiMinorAxisB * A * (sigma - deltaSigma)
result = s / 1000
break
}
}
if math.IsNaN(result) {
return result, fmt.Errorf("failed to converge for Point(Latitude: %f, Lomgitude: %f) and Point(Latitude: %f, Lomgitude: %f)", lat1, long1, lat2, long2)
}
return result, nil
}
//Length3D calculates the lenght of given points list including elevation distance
func LengthVincenty(locs []Point) (float64, error) {
return lengthVincenty(locs)
}
I've adapted your style to use the funcs but would be happy to help that a user may provide a custom calculation method.
P.s. The used method for Vincenty is from: https://github.com/asmarques/geodist/blob/master/vincenty.go (Both repos are aligned with Apache 2.0 License; so shouldn't be a problem to use)
Is there anything I can do to assist whatever getting merged in?
Hi @EliDavis3D , I'm sorry, the posts are quite old and I do not remember all stuff. But I see that I've more or less clearly described what to do ;-) Looking at my repos and my commits for the forked repo at https://github.com/tkrajina/gpxgo/commit/42d8413575fa83909dfb5b8013bac7ab06ba8c42 I would say that the implementation is already done and must be only pushed back to the original repo.
Because I have in my implementation some references to my forked repo I would say do the following:
(1) Fork this repo (2) Look at my implementation / commit; copy the code in your forked repo (3) Test the implementation with the referenced files (4) Push the commit back to the original repo
Does this make sense? :-)
Hi @tkrajina , thanks for the library! Good work!
Referenced gpx file: https://gist.github.com/mbecker/a44881bfe29b0982fac6c69cae498125
Looking at the length2d and length3d of the referenced gpx file I've noticed that Strava is showing me different values.
The Strava values are as follows:
The gpxgo values are as follows:
So looking into your code I've noticed the following:
The method "HaversineDistance" is only used if the distance between the two points are too distant: https://github.com/tkrajina/gpxgo/blob/1b1f71eb1b590891d4ce4e1325f8bacbb20d8d4c/gpx/geo.go#L188 -> Valid aproach since you are saying in gpxpy that the calculation is too complex (https://github.com/tkrajina/gpxpy/issues/9) -> Maybe the user should decide which operation to use? Just an idea!
In gpxy you've changed the value for EARTH_RADIUS as follows: https://github.com/jedie/gpxpy/commit/b371de31826cfecc049a57178d168b04a7f6d0d8
-> In gpxgo it's still the "old" value: https://github.com/tkrajina/gpxgo/blob/1b1f71eb1b590891d4ce4e1325f8bacbb20d8d4c/gpx/geo.go#L14
-> Since the value of "earthRadius" is only in the func "HaversineDistance" it's not used for calculation -> I've manually changed the function to use "HaversineDistance" and gets the follwing results:
(Not sure about this) In the func "HaversineDistance" the elevation is not used. Is that by design of the formula?
Digging into the different formulas for calculating the distance between two points I've read that the formula "Vincenty" should be more accurate (https://en.wikipedia.org/wiki/Vincenty%27s_formulae). Using that formula I do get the following value:
(5) I've inserted all points into a postgresql postgis database and created a LINESTRING() with the points. Postgresql Potsgis is calculating the length of the LINESTRING as follows:
So to summarize my observation of the different calculations are as follows:
Referencing to the gpxpy issue 123 (https://github.com/tkrajina/gpxpy/issues/123) and using the generated generated_10km,gpx and half_equator.gpx the results are as follows:
Based on the test results I would propose to do the following:
Starting from my original question, I still don't kow "how strava is calculating the distance" :-D
What do you think of the small updates?
Thanks again for your work and best regards.