From 3ef56d5fc669258314ba0445871cb24e1afd340f Mon Sep 17 00:00:00 2001 From: sumsam Date: Thu, 26 Aug 2021 18:48:13 +0200 Subject: [PATCH 1/4] feat: add implementation for haversine distance for geo-coordinates --- planar/distance.go | 24 +++++ planar/distance_from.go | 2 +- planar/haversine_distance_from.go | 169 ++++++++++++++++++++++++++++++ 3 files changed, 194 insertions(+), 1 deletion(-) create mode 100644 planar/haversine_distance_from.go diff --git a/planar/distance.go b/planar/distance.go index a0304a7..adaf6ce 100644 --- a/planar/distance.go +++ b/planar/distance.go @@ -6,6 +6,11 @@ import ( "github.com/paulmach/orb" ) +const ( + // EarthRadius is the Earth's radius in kilometers + EarthRadius = 6371 +) + // Distance returns the distance between two points in 2d euclidean geometry. func Distance(p1, p2 orb.Point) float64 { d0 := (p1[0] - p2[0]) @@ -19,3 +24,22 @@ func DistanceSquared(p1, p2 orb.Point) float64 { d1 := (p1[1] - p2[1]) return d0*d0 + d1*d1 } + +// HaversineDistance calculates the distance between two points on earth in kilometers +// Implementation from: http://www.movable-type.co.uk/scripts/latlong.html +func HaversineDistance(p1, p2 orb.Point) float64 { + dLat := (p2[1] - p1[1]) * (math.Pi / 180.0) + dLon := (p2[0] - p1[0]) * (math.Pi / 180.0) + + lat1 := p1[1] * (math.Pi / 180.0) + lat2 := p2[1] * (math.Pi / 180.0) + + a1 := math.Sin(dLat/2) * math.Sin(dLat/2) + a2 := math.Sin(dLon/2) * math.Sin(dLon/2) * math.Cos(lat1) * math.Cos(lat2) + + a := a1 + a2 + + c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a)) + + return EarthRadius * c +} diff --git a/planar/distance_from.go b/planar/distance_from.go index da4fffc..1cd49c7 100644 --- a/planar/distance_from.go +++ b/planar/distance_from.go @@ -12,7 +12,7 @@ func DistanceFromSegment(a, b, point orb.Point) float64 { return math.Sqrt(DistanceFromSegmentSquared(a, b, point)) } -// DistanceFromSegmentSquared returns point's squared distance from the segement [a, b]. +// DistanceFromSegmentSquared returns point's squared distance from the segment [a, b]. func DistanceFromSegmentSquared(a, b, point orb.Point) float64 { x := a[0] y := a[1] diff --git a/planar/haversine_distance_from.go b/planar/haversine_distance_from.go new file mode 100644 index 0000000..bdc240e --- /dev/null +++ b/planar/haversine_distance_from.go @@ -0,0 +1,169 @@ +package planar + +import ( + "fmt" + "math" + + "github.com/paulmach/orb" +) + +// HaversineDistanceFromSegment returns point's haversine distance on earth from the segment [a, b] in kilometers. +func HaversineDistanceFromSegment(a, b, point orb.Point) float64 { + x := a[0] + y := a[1] + dx := b[0] - x + dy := b[1] - y + + if dx != 0 || dy != 0 { + t := ((point[0]-x)*dx + (point[1]-y)*dy) / (dx*dx + dy*dy) + + if t > 1 { + x = b[0] + y = b[1] + } else if t > 0 { + x += dx * t + y += dy * t + } + } + + dx = point[0] - x + dy = point[1] - y + + return HaversineDistance(point,orb.Point{x,y}) +} + +// HaversineDistanceFrom returns the distance on earth from the boundary of the geometry in +// kilometers +func HaversineDistanceFrom(g orb.Geometry, p orb.Point) float64 { + d, _ := HaversineDistanceFromWithIndex(g, p) + return d +} + + +// HaversineDistanceFromWithIndex returns the minimum haversine distance on earth in kilometers +// from the boundary of the geometry plus the index of the sub-geometry +// that was the match. +func HaversineDistanceFromWithIndex(g orb.Geometry, p orb.Point) (float64, int) { + if g == nil { + return math.Inf(1), -1 + } + + switch g := g.(type) { + case orb.Point: + return HaversineDistance(g, p), 0 + case orb.MultiPoint: + return multiPointHaversineDistanceFrom(g, p) + case orb.LineString: + return lineStringHaversineDistanceFrom(g, p) + case orb.MultiLineString: + dist := math.Inf(1) + index := -1 + for i, ls := range g { + if d, _ := lineStringHaversineDistanceFrom(ls, p); d < dist { + dist = d + index = i + } + } + + return dist, index + case orb.Ring: + return lineStringHaversineDistanceFrom(orb.LineString(g), p) + case orb.Polygon: + return polygonHaversineDistanceFrom(g, p) + case orb.MultiPolygon: + dist := math.Inf(1) + index := -1 + for i, poly := range g { + if d, _ := polygonHaversineDistanceFrom(poly, p); d < dist { + dist = d + index = i + } + } + + return dist, index + case orb.Collection: + dist := math.Inf(1) + index := -1 + for i, ge := range g { + if d, _ := HaversineDistanceFromWithIndex(ge, p); d < dist { + dist = d + index = i + } + } + + return dist, index + case orb.Bound: + return HaversineDistanceFromWithIndex(g.ToRing(), p) + } + + panic(fmt.Sprintf("geometry type not supported: %T", g)) +} + +func multiPointHaversineDistanceFrom(mp orb.MultiPoint, p orb.Point) (float64, int) { + dist := math.Inf(1) + index := -1 + + for i := range mp { + if d := HaversineDistance(mp[i], p); d < dist { + dist = d + index = i + } + } + + return dist, index +} + +func lineStringHaversineDistanceFrom(ls orb.LineString, p orb.Point) (float64, int) { + dist := math.Inf(1) + index := -1 + + for i := 0; i < len(ls)-1; i++ { + if d := segmentHaversineDistanceFrom(ls[i], ls[i+1], p); d < dist { + dist = d + index = i + } + } + + return dist, index +} + +func polygonHaversineDistanceFrom(p orb.Polygon, point orb.Point) (float64, int) { + if len(p) == 0 { + return math.Inf(1), -1 + } + + dist, index := lineStringHaversineDistanceFrom(orb.LineString(p[0]), point) + for i := 1; i < len(p); i++ { + d, i := lineStringHaversineDistanceFrom(orb.LineString(p[i]), point) + if d < dist { + dist = d + index = i + } + } + + return dist, index +} + +func segmentHaversineDistanceFrom(p1, p2, point orb.Point) float64 { + x := p1[0] + y := p1[1] + dx := p2[0] - x + dy := p2[1] - y + + if dx != 0 || dy != 0 { + t := ((point[0]-x)*dx + (point[1]-y)*dy) / (dx*dx + dy*dy) + + if t > 1 { + x = p2[0] + y = p2[1] + } else if t > 0 { + x += dx * t + y += dy * t + } + } + + dx = point[0] - x + dy = point[1] - y + + return HaversineDistance(point,orb.Point{x,y}) +} From 540148b5cd08a8d72d8fd3eed8128582f83b03f2 Mon Sep 17 00:00:00 2001 From: sumsam Date: Tue, 28 Dec 2021 17:51:37 +0100 Subject: [PATCH 2/4] move geo calculations in geo package --- {planar => geo}/haversine_distance_from.go | 16 +++++----------- planar/distance.go | 18 ------------------ 2 files changed, 5 insertions(+), 29 deletions(-) rename {planar => geo}/haversine_distance_from.go (92%) diff --git a/planar/haversine_distance_from.go b/geo/haversine_distance_from.go similarity index 92% rename from planar/haversine_distance_from.go rename to geo/haversine_distance_from.go index bdc240e..457f89d 100644 --- a/planar/haversine_distance_from.go +++ b/geo/haversine_distance_from.go @@ -1,4 +1,4 @@ -package planar +package geo import ( "fmt" @@ -26,10 +26,7 @@ func HaversineDistanceFromSegment(a, b, point orb.Point) float64 { } } - dx = point[0] - x - dy = point[1] - y - - return HaversineDistance(point,orb.Point{x,y}) + return DistanceHaversine(point,orb.Point{x,y}) } // HaversineDistanceFrom returns the distance on earth from the boundary of the geometry in @@ -50,7 +47,7 @@ func HaversineDistanceFromWithIndex(g orb.Geometry, p orb.Point) (float64, int) switch g := g.(type) { case orb.Point: - return HaversineDistance(g, p), 0 + return DistanceHaversine(g, p), 0 case orb.MultiPoint: return multiPointHaversineDistanceFrom(g, p) case orb.LineString: @@ -104,7 +101,7 @@ func multiPointHaversineDistanceFrom(mp orb.MultiPoint, p orb.Point) (float64, i index := -1 for i := range mp { - if d := HaversineDistance(mp[i], p); d < dist { + if d := DistanceHaversine(mp[i], p); d < dist { dist = d index = i } @@ -162,8 +159,5 @@ func segmentHaversineDistanceFrom(p1, p2, point orb.Point) float64 { } } - dx = point[0] - x - dy = point[1] - y - - return HaversineDistance(point,orb.Point{x,y}) + return DistanceHaversine(point,orb.Point{x,y}) } diff --git a/planar/distance.go b/planar/distance.go index adaf6ce..a1ef3a2 100644 --- a/planar/distance.go +++ b/planar/distance.go @@ -25,21 +25,3 @@ func DistanceSquared(p1, p2 orb.Point) float64 { return d0*d0 + d1*d1 } -// HaversineDistance calculates the distance between two points on earth in kilometers -// Implementation from: http://www.movable-type.co.uk/scripts/latlong.html -func HaversineDistance(p1, p2 orb.Point) float64 { - dLat := (p2[1] - p1[1]) * (math.Pi / 180.0) - dLon := (p2[0] - p1[0]) * (math.Pi / 180.0) - - lat1 := p1[1] * (math.Pi / 180.0) - lat2 := p2[1] * (math.Pi / 180.0) - - a1 := math.Sin(dLat/2) * math.Sin(dLat/2) - a2 := math.Sin(dLon/2) * math.Sin(dLon/2) * math.Cos(lat1) * math.Cos(lat2) - - a := a1 + a2 - - c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a)) - - return EarthRadius * c -} From 4043e59daa51afb26a30fcdbf37c5dd68d94d359 Mon Sep 17 00:00:00 2001 From: sumsam Date: Tue, 28 Dec 2021 17:54:06 +0100 Subject: [PATCH 3/4] update the comments: km->m --- geo/haversine_distance_from.go | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/geo/haversine_distance_from.go b/geo/haversine_distance_from.go index 457f89d..038f77b 100644 --- a/geo/haversine_distance_from.go +++ b/geo/haversine_distance_from.go @@ -7,7 +7,7 @@ import ( "github.com/paulmach/orb" ) -// HaversineDistanceFromSegment returns point's haversine distance on earth from the segment [a, b] in kilometers. +// HaversineDistanceFromSegment returns point's haversine distance on earth from the segment [a, b] in meters. func HaversineDistanceFromSegment(a, b, point orb.Point) float64 { x := a[0] y := a[1] @@ -30,14 +30,14 @@ func HaversineDistanceFromSegment(a, b, point orb.Point) float64 { } // HaversineDistanceFrom returns the distance on earth from the boundary of the geometry in -// kilometers +// meters func HaversineDistanceFrom(g orb.Geometry, p orb.Point) float64 { d, _ := HaversineDistanceFromWithIndex(g, p) return d } -// HaversineDistanceFromWithIndex returns the minimum haversine distance on earth in kilometers +// HaversineDistanceFromWithIndex returns the minimum haversine distance on earth in meters // from the boundary of the geometry plus the index of the sub-geometry // that was the match. func HaversineDistanceFromWithIndex(g orb.Geometry, p orb.Point) (float64, int) { From b2e5441fae4fc75488c6bb8f66a0f63b4b8bb5a0 Mon Sep 17 00:00:00 2001 From: sumsam Date: Tue, 28 Dec 2021 18:01:57 +0100 Subject: [PATCH 4/4] remove unused constant --- planar/distance.go | 5 ----- 1 file changed, 5 deletions(-) diff --git a/planar/distance.go b/planar/distance.go index a1ef3a2..10871f9 100644 --- a/planar/distance.go +++ b/planar/distance.go @@ -6,11 +6,6 @@ import ( "github.com/paulmach/orb" ) -const ( - // EarthRadius is the Earth's radius in kilometers - EarthRadius = 6371 -) - // Distance returns the distance between two points in 2d euclidean geometry. func Distance(p1, p2 orb.Point) float64 { d0 := (p1[0] - p2[0])