feat: add implementation for haversine distance for geo-coordinates

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
+}

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]

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})
+}