Eaxcact distance

[mbecker]

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:

Strava - Before correction of GPS: 68.69 km
Strava - After correction of GPS : 68.93 km (manual correction of Strava)

The gpxgo values are as follows:

gpxgo - length2d: 68.828110 km
gpxgo - length3d: 68.840900 km

So looking into your code I've noticed the following:

1. The method "HaversineDistance" is only used if the distance between the two points are too distant:

   gpxgo/gpx/geo.go

   Line 188 in 1b1f71e

   if haversine || absLat > 0.2 || absLon > 0.2 {

   
   -> Valid aproach since you are saying in gpxpy that the calculation is too complex (Correct distance between points with haversine formula gpxpy#9)
   -> Maybe the user should decide which operation to use? Just an idea!

2. In gpxy you've changed the value for EARTH_RADIUS as follows: jedie/gpxpy@b371de3

EARTH_RADIUS = 6378.137 * 1000)

-> In gpxgo it's still the "old" value:

gpxgo/gpx/geo.go

Line 14 in 1b1f71e

const earthRadius = 6371 * 1000

const earthRadius = 6371 * 1000)

-> 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:

gpxgo - lengt2d HaversineDistance: 68.874521 km
gpxgo - lengt2d HaversineDistance (updated new value for earthRadius): 68.951676 km

3. (Not sure about this)
   In the func "HaversineDistance" the elevation is not used. Is that by design of the formula?

4. 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:

func DistanceVincenty(p1 Point, p2 Point, withElevation bool) (float64, error) {
	if p1.Lat == p2.Lat && p1.Long == p2.Long {
		return 0, nil
	}

	U1 := math.Atan((1 - f) * math.Tan(toRadians(p1.Lat)))
	U2 := math.Atan((1 - f) * math.Tan(toRadians(p2.Lat)))
	L := toRadians(p2.Long - p1.Long)
	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 := (f / 16) * cosSqrAlpha * (4 + f*(4-3*cosSqrAlpha))
		lambda = L + (1-C)*f*sinAlpha*(sigma+C*sinSigma*(cos2sigmam+C*cosSigma*(-1+2*math.Pow(cos2sigmam, 2))))

		if math.Abs(lambda-curLambda) < epsilon {
			uSqr := cosSqrAlpha * ((math.Pow(a, 2) - math.Pow(b, 2)) / math.Pow(b, 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 := b * A * (sigma - deltaSigma)
			result = s

			break
		}
	}

	if math.IsNaN(result) {
		return result, fmt.Errorf("failed to converge for %v and %v", p1, p2)
	}

	if !withElevation {
		return result, nil
	}

	eleDiff := 0.0
	if p1.ElevationNotNull && p2.ElevationNotNull {
		eleDiff = p1.Elevation - p2.Elevation
	}

	return math.Sqrt(math.Pow(result, 2) + math.Pow(eleDiff, 2)), nil
}
Distance Vincenty: 69.040838 km (without elevation)
Distance Vincenty: 69.053588 km (with elevation)

(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:

select ST_LengthSpheroid(line,'SPHEROID["WGS 84",6378137,298.257223563]')/1000 from strava_test.strava_tracks where gpxid = 1040;
-- Result: 69.0408387739 km

So to summarize my observation of the different calculations are as follows:

Strava - Before correction of GPS: 68.69 km
Strava - After correction of GPS : 68.93 km (manual correction of Strava)

gpxgo - length2d: 68.828110 km
gpxgo - length3d: 68.840900 km

gpxgo - lengt2d HaversineDistance: 68.874521 km (without elevation)
gpxgo - lengt2d HaversineDistance (updated new value for earthRadius): 68.951676 km (without elevation)

Distance Vincenty: 69.040838 km (without elevation)
Distance Vincenty: 69.053588 km (with elevation)

postgresql postgis - 69.0408387739 km

Referencing to the gpxpy issue 123 (tkrajina/gpxpy#123) and using the generated generated_10km,gpx and half_equator.gpx the results are as follows:

-- File: generated_10km.gpx --
gpxgo - length2d: 9.988810
gpxgo - lenght3d: 9.988810
gpxgo - lengt2d HaversineDistance: 9.988810
gpxgo - lengt2d HaversineDistance (updated new value for earthRadius):: 10.000000 !!! As in the python test !!!
Distance Vincenty: 10.000000 (without elevation) !!! As in the python test !!!
Distance Vincenty: 10.000000 (with elevation) !!! As in the python test !!!

-- File: half_equator.gpx --
gpxgo - length2d: 20015.086796
gpxgo - lenght3d: 20015.086796
gpxgo - lengt2d HaversineDistance: 20015.086796
gpxgo - lengt2d HaversineDistance (updated new value for earthRadius):: 20037.508343 !!! As in the python test !!!
Distance Vincenty: 20037.508343 (without elevation) !!! As in the python test !!!
Distance Vincenty: 20037.508343 (with elevation) !!! As in the python test !!!

Based on the test results I would propose to do the following:

1. Update the value for "earthRadius" (

   gpxgo/gpx/geo.go

   Line 14 in 1b1f71e

   const earthRadius = 6371 * 1000

   ) to be in line with gpxpy
2. Let the user decide which formula to use (either the fast one by yourself or Haversine)
3. Since the result of postgresql postgis and the formula of Vincenty are the same and the Vincenty and python are the same, implement the formula and let the user decide which formula to use

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.

[tkrajina]

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.

[mbecker]

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)

[EliDavis3D]

Is there anything I can do to assist whatever getting merged in?

[mbecker]

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 42d8413 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? :-)