search

golang / geo

S2 geometry library in Go
Apache License 2.0
1.69k stars 182 forks source link

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

Closed idc77 closed 2 years ago

idc77 commented 2 years ago 
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 commented 2 years ago

@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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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.