Q: How do you calculate the distance in km between 2 lat/lng points?

[idc77]

package main

import (
	"fmt"

	"github.com/golang/geo/s1"
	"github.com/golang/geo/s2"
)

func main() {
	lat1 := 9.1829321
	lng1 := 48.7758459
	lat2 := 5.03131
	lng2 := 39.542448

	pg1 := s2.PointFromLatLng(s2.LatLng{
		Lat: s1.Angle(lat1),
		Lng: s1.Angle(lng1),
	})
	pg2 := s2.PointFromLatLng(s2.LatLng{
		Lat: s1.Angle(lat2),
		Lng: s1.Angle(lng2),
	})
	sdist := pg1.Distance(pg2)
	fmt.Printf("s2 distance: %f\n", sdist)

}

s2 distance: 1.499294

This can't be true. It's nearly 1080km.

Can you please tell how to do it?

[doodlesbykumbi]

@idc77 the method Point#Distance returns an Angle. To get the distance in km you'll need to multiply by the radius of the earth.

There's also some additional changes needed from your code to construct a Point from lat-long, see below.

package main

import (
	"fmt"

	"github.com/golang/geo/s2"
)

const earthRadiusKm = 6371.01 // See https://github.com/golang/geo/blob/master/s2/s2_test.go#L46-L51

func main() {
	lat1 := 9.1829321
	lng1 := 48.7758459
	lat2 := 5.03131
	lng2 := 39.542448

	pg1 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat1, lng1))
	pg2 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat2, lng2))

       sdist := pg2.Distance(pg1) * earthRadiusKm
	fmt.Printf("s2 distance %f\n", sdist)
}

[idc77]

Hi @doodlesbykumbi thank you very much for your reply.

Something is really not adding up.
I'm comparing the results of MongoDB's distance calculation and
"github.com/kellydunn/golang-geo"
and
this package.

package main

import (
	"fmt"

	"github.com/golang/geo/s2"
	geo2 "github.com/kellydunn/golang-geo"
)

const earthRadiusKm = 6371.01 // See https://github.com/golang/geo/blob/master/s2/s2_test.go#L46-L51

func main() {
	lat1 := 9.1829321
	lng1 := 48.7758459
	// lat2 := 5.03131
	// lng2 := 39.542448
	lat2 := 0.272091
	lng2 := 45.937435
	p1 := geo2.NewPoint(9.1829321, 48.7758459)
	p2 := geo2.NewPoint(5.03131, 39.542448)

	dist := p1.GreatCircleDistance(p2)
	fmt.Printf("great circle distance: %f\n", dist)

	pg1 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat1, lng1))
	pg2 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat2, lng2))
	sdist := pg1.Distance(pg2) * earthRadiusKm
	odist := pg2.Distance(pg1) * earthRadiusKm
	fmt.Printf("s2 distance: %f\n", sdist)
	fmt.Printf("s2 r-distance: %f\n", odist)
}

great circle distance: 1118.298549
s2 distance: 1039.471777
s2 r-distance: 1039.471777

however
a $near query in MongoDB resulsts in 741940.633007093 meters

			geoStage := bson.D{{
				"$geoNear", bson.D{
					{"near", profile.Location},
					{"distanceField", "distance"},
					{"minDistance", 1},
					{"maxDistance", radius},
					{"query", bson.D{
						{"subject_id", m.TargetID},
					}},
					{"spherical", true},
				},
			},
			}

This is just to get the distance of "my profile" to "target profile".
Maybe it's because of the "spherical" = true

https://www.mongodb.com/docs/manual/reference/operator/aggregation/geoNear/

spherical	boolean	Optional. Determines how MongoDB calculates the distance between two points:When true, MongoDB uses $nearSphere semantics and calculates distances using spherical geometry.When false, MongoDB uses $near semantics: spherical geometry for 2dsphere indexes and planar geometry for 2d indexes.Default: false.

spherical

boolean

Optional. Determines how MongoDB calculates the distance between two points:

When true, MongoDB uses [$nearSphere](https://www.mongodb.com/docs/manual/reference/operator/query/nearSphere/#mongodb-query-op.-nearSphere) semantics and calculates distances using spherical geometry.
When false, MongoDB uses [$near](https://www.mongodb.com/docs/manual/reference/operator/query/near/#mongodb-query-op.-near) semantics: spherical geometry for [2dsphere](https://www.mongodb.com/docs/manual/core/2dsphere/) indexes and planar geometry for [2d](https://www.mongodb.com/docs/manual/core/2d/) indexes.

Default: false.

[rsned]

It's definitely something on the mongo side. If I use Google Maps and "Measure Distance", I get close to 1100 km for
p1 = 9.1829321, 48.7758459
p2 = 5.03131, 39.542448

https://www.google.com/maps/dir/5.03131,+39.542448/9.1829321,48.7758459/@7.8185624,42.7915689,7z/data=!4m7!4m6!1m3!2m2!1d39.542448!2d5.03131!1m0!3e2

And ~1036 for
p1 = 9.1829321, 48.7758459
p2 = 0.272091, 45.937435

https://www.google.com/maps/dir/9.1829321,+48.7758459/0.272091,45.937435/@4.7679831,47.1160662,7z/data=!4m7!4m6!1m3!2m2!1d48.7758459!2d9.1829321!1m0!3e2

Are you sure the points you are using in the mongo query are the same as in the s2 call?

[idc77]

Sorry I posted the wrong source.
I will post a new source

package main

import (
	"fmt"

	"github.com/golang/geo/s2"
	geo2 "github.com/kellydunn/golang-geo"
)

const earthRadiusKm = 6371.01 // See https://github.com/golang/geo/blob/master/s2/s2_test.go#L46-L51

func main() {
	lat1 := 9.1829321
	lng1 := 48.7758459
	// lat2 := 5.03131
	// lng2 := 39.542448#
	lat2 := 15.39967
	lng2 := 45.749428
	p1 := geo2.NewPoint(lat1, lng1)
	p2 := geo2.NewPoint(lat2, lng2)

	dist := p1.GreatCircleDistance(p2)
	fmt.Printf("great circle distance: %f\n", dist)

	pg1 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat1, lng1))
	pg2 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat2, lng2))
	sdist := pg1.Distance(pg2) * earthRadiusKm
	odist := pg2.Distance(pg1) * earthRadiusKm
	fmt.Printf("s2 distance: %f\n", sdist)
	fmt.Printf("s2 r-distance: %f\n", odist)
}

great circle distance: 765.406106
s2 distance: 765.407307
s2 r-distance: 765.407307

The non-spherical result in mongodb is 577683.4757611528 aka about 577 km.

p1 is a real location
p2 is made up of randomly generated values

Using google maps I couldn't get the exact location of p1
coordinates should be reverse in mongodb

package models

// Location is a GeoJSON type.
type Location struct {
	Type        string    `json:"type" bson:"type"`
	Coordinates []float64 `json:"coordinates" bson:"coordinates"`
}

// NewPoint returns a GeoJSON Point with longitude and latitude.
func NewPoint(long, lat float64) Location {
	return Location{
		"Point",
		[]float64{long, lat},
	}
}

mongodb saves the points in reverse order long, lat.
Maybe that is the mystery, I need to go for a walk to clear my mind, will return later.

[idc77]

Okay, that is indeed the case.

package main

import (
	"fmt"

	"github.com/golang/geo/s2"
	geo2 "github.com/kellydunn/golang-geo"
)

const earthRadiusKm = 6371.01 // See https://github.com/golang/geo/blob/master/s2/s2_test.go#L46-L51

func main() {
	lat1 := 9.1829321
	lng1 := 48.7758459
	lat2 := 15.39967
	lng2 := 45.749428
	// p1 := geo2.NewPoint(lat1, lng1)
	p1 := geo2.NewPoint(lng1, lat1)
	// p2 := geo2.NewPoint(lat2, lng2)
	p2 := geo2.NewPoint(lng2, lat2)

	dist := p1.GreatCircleDistance(p2)
	fmt.Printf("great circle distance: %f\n", dist)

	// pg1 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat1, lng1))
	pg1 := s2.PointFromLatLng(s2.LatLngFromDegrees(lng1, lat1))
	// pg2 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat2, lng2))
	pg2 := s2.PointFromLatLng(s2.LatLngFromDegrees(lng2, lat2))
	sdist := pg1.Distance(pg2) * earthRadiusKm
	// odist := pg2.Distance(pg1) * earthRadiusKm
	fmt.Printf("s2 distance: %f\n", sdist)
	// fmt.Printf("s2 r-distance: %f\n", odist)
}

I swapped longitude and latitude and the resulsts are drumroll

great circle distance: 577.040408
s2 distance: 577.041313

MongoDB value: 577683.4757611528

So MongoDB requires long, lat format, but does the calculation wrong internall by mistaking lat for long and the other way around.

That's if I didn't do anything wrong.
I can't believe no one hit this issue before.

Thank you all for you help :) I'll report it to MongoDB.
Well... I would if they didn't make it such a PITA to do.
I guess it will remain an issue until someone with a MongoDB Jira account that doesn't use 2FA stumbles upon this.

[idc77]

It wasn't mongo. It was me. I'm the mongo ;P Self-knowledge is the 1st step to enlightenment
But there are still slightly different results.

great circle distance: 577.040408
s2 distance: 577.041313
mongodb distance: 577683.475761

package main

import (
	"context"
	"fmt"
	"log"
	"time"

	"github.com/golang/geo/s2"
	geo "github.com/kellydunn/golang-geo"
	"go.mongodb.org/mongo-driver/bson"
	"go.mongodb.org/mongo-driver/bson/primitive"
	"go.mongodb.org/mongo-driver/mongo"
	mopts "go.mongodb.org/mongo-driver/mongo/options"
	"go.mongodb.org/mongo-driver/mongo/readpref"
	"go.mongodb.org/mongo-driver/x/bsonx"
)

const earthRadiusKm = 6371.01 // See https://github.com/golang/geo/blob/master/s2/s2_test.go#L46-L51

var (
	Database          = "testgeo"
	DefaultTimeout    = time.Second * 10
	ProfileCollection = "profile"
)

type (
	Profile struct {
		ID       primitive.ObjectID `bson:"_id,omitempty" json:"id"`
		Location Location           `bson:"location" json:"location"`
	}
	ProfileResult struct {
		ID       primitive.ObjectID `bson:"_id,omitempty" json:"id"`
		Location Location           `bson:"location" json:"location"`
		Distance float64            `bson:"distance" json:"distance"`
	}
	Location struct {
		Type        string    `json:"type" bson:"type"`
		Coordinates []float64 `json:"coordinates" bson:"coordinates"`
	}
)

func main() {
	// lat1 := 9.1829321
	// lng1 := 48.7758459
	// lat2 := 45.749428
	//lng2 := 15.39967
	// lat1 := 9.653194
	// lng1 := 48.709797
	// lat2 := 9.182313
	// lng2 := 48.783187
	lng1 := 9.1829321
	lat1 := 48.7758459

	lng2 := 15.39967
	lat2 := 45.749428

	p1 := geo.NewPoint(lat1, lng1)
	// p1 := geo2.NewPoint(lng1, lat1)
	p2 := geo.NewPoint(lat2, lng2)
	// p2 := geo2.NewPoint(lng2, lat2)

	dist := p1.GreatCircleDistance(p2)
	fmt.Printf("great circle distance: %f\n", dist)

	pg1 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat1, lng1))
	// pg1 := s2.PointFromLatLng(s2.LatLngFromDegrees(lng1, lat1))
	pg2 := s2.PointFromLatLng(s2.LatLngFromDegrees(lat2, lng2))
	// pg2 := s2.PointFromLatLng(s2.LatLngFromDegrees(lng2, lat2))
	sdist := pg1.Distance(pg2) * earthRadiusKm
	// odist := pg2.Distance(pg1) * earthRadiusKm
	fmt.Printf("s2 distance: %f\n", sdist)
	// fmt.Printf("s2 r-distance: %f\n", odist)
	// fmt.Printf("mongo result: %f\n", 5165738.082454552)

	mp1 := new(Profile)
	mp1.Location = NewPoint(lng1, lat1)

	mp2 := new(Profile)
	mp2.Location = NewPoint(lng2, lat2)

	client, e := mongoInit()
	if e != nil {
		log.Fatal(e.Error())
	}
	defer func() {
		if e = client.Disconnect(context.Background()); e != nil {
			panic(e)
		}
	}()

	res1, e := CreateProfile(client, mp1)
	if e != nil {
		log.Fatal(e.Error())
	}

	res2, e := CreateProfile(client, mp2)
	if e != nil {
		log.Fatal(e.Error())
	}

	mp1.ID = res1.InsertedID.(primitive.ObjectID)
	mp2.ID = res2.InsertedID.(primitive.ObjectID)

	var radius int32
	radius = int32(6371) * 1000
	geoStage := bson.D{{
		"$geoNear", bson.D{
			{"near", mp1.Location},
			{"distanceField", "distance"},
			{"minDistance", 1},
			{"maxDistance", radius},
			{"spherical", true},
		},
	},
	}
	result, e := AggregateProfile(client, mongo.Pipeline{geoStage})
	if e != nil {
		return
	}
	fmt.Printf("mongodb distance: %f\n", result[0].Distance)

}

func mongoInit() (*mongo.Client, error) {
	ctx, cancel := context.WithTimeout(context.Background(), DefaultTimeout)
	defer cancel()
	client, e := mongo.Connect(ctx, mopts.Client().ApplyURI("mongodb://localhost:27017"))
	if e != nil {
		return nil, e
	}

	ctx, cancel = context.WithTimeout(context.Background(), DefaultTimeout)
	defer cancel()
	e = client.Ping(ctx, readpref.Primary())
	if e != nil {
		return nil, e
	}

	// drop database
	_ = client.Database(Database).Drop(context.Background())

	// create geoJSON 2dsphere index
	ctx, cancel = context.WithTimeout(context.Background(), DefaultTimeout)
	defer cancel()

	db := client.Database(Database)
	indexOpts := mopts.CreateIndexes().SetMaxTime(DefaultTimeout)

	// Index to location 2dsphere type.
	pointIndexModel := mongo.IndexModel{
		Options: mopts.Index(),
		Keys: bsonx.MDoc{
			"location": bsonx.String("2dsphere"),
		},
	}

	profileIndexes := db.Collection("profile").Indexes()

	// _, e = profileIndexes.CreateOne(ctx, pointIndexModel, indexOpts)
	_, e = profileIndexes.CreateOne(ctx, pointIndexModel, indexOpts)
	if e != nil {
		return nil, e
	}

	return client, nil
}

// NewPoint returns a GeoJSON Point with longitude and latitude.
func NewPoint(long, lat float64) Location {
	return Location{
		"Point",
		[]float64{long, lat},
	}
}

func CreateProfile(client *mongo.Client, m *Profile) (*mongo.InsertOneResult, error) {
	ctx, cancel := context.WithTimeout(context.Background(), DefaultTimeout)
	defer cancel()
	c := client.Database(Database).Collection(ProfileCollection)
	return c.InsertOne(ctx, m)
}

func AggregateProfile(client *mongo.Client, pipeline interface{}) ([]*ProfileResult, error) {
	ctx, cancel := context.WithTimeout(context.Background(), DefaultTimeout)
	defer cancel()
	c := client.Database(Database).Collection(ProfileCollection)
	cursor, e := c.Aggregate(ctx, pipeline)
	if e != nil {
		return nil, e
	}
	ctx, cancel = context.WithTimeout(context.Background(), DefaultTimeout)
	defer cancel()
	var m []*ProfileResult
	if e := cursor.All(ctx, &m); e != nil {
		return nil, e
	}
	return m, nil
}

But it's close enough and the users won't need very high accuracy in their results.
Would be interesting to know what postgis returns there.
Anyhow another day wasted spent on self-enlightenment ;P

[idc77]

If anyone is still reading, the postgis 3.2.1 on postgresql 14 result ahaha

SELECT ST_Distance(
               ST_GeographyFromText('Point(9.1829321 48.7758459)'),
               ST_GeographyFromText('Point(15.39967 45.749428)'))
           AS geography_distance;

578123.95280377

4 approaches, 4 different results
Very entertaining
In reality it's probably somewhere around 490km. I took 2 real life points and ran them through the above program and the result was 52km. Using Google Maps and walking,the "not a straight line" point to point distance was 43km. Maybe Google Maps take the distance from the sea level into account or from the center of the earth. Maybe it uses vectors to calculate.
But those calculations are all 2d on a sphere, right?

Somewhat interesting

1 row retrieved starting from 1 in 39 ms (execution: 3 ms, fetching: 36 ms)

That's about 200ms faster than MongoDB (but then again I assume MongoDB was going through 100k datasets and comparing the distance to min and max distance).

Have a nice day. This was my last post in this issue.