feat: add implementation for haversine distance for geo-coordinates

geo/haversine_distance_from.go

@@ -0,0 +1,163 @@
+package geo
+
+import (
+	"fmt"
+	"math"
+
+	"github.com/paulmach/orb"
+)
+
+// 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]
+	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
+		}
+	}
+
+	return DistanceHaversine(point,orb.Point{x,y})
+}
+
+// HaversineDistanceFrom returns the distance on earth from the boundary of the geometry in
+// 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 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) {
+	if g == nil {
+		return math.Inf(1), -1
+	}
+
+	switch g := g.(type) {
+	case orb.Point:
+		return DistanceHaversine(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 := DistanceHaversine(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
+		}
+	}
+
+	return DistanceHaversine(point,orb.Point{x,y})
+}

planar/distance.go

@@ -19,3 +19,4 @@ func DistanceSquared(p1, p2 orb.Point) float64 {
 	d1 := (p1[1] - p2[1])
 	return d0*d0 + d1*d1
 }
+

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]