diff --git a/geo/distance.go b/geo/distance.go index 44d81b5..a04aa9b 100644 --- a/geo/distance.go +++ b/geo/distance.go @@ -69,3 +69,50 @@ func Midpoint(p, p2 orb.Point) orb.Point { return r } + +// PointAtBearingAndDistance returns the point at the given bearing and distance in meters from the point +func PointAtBearingAndDistance(p orb.Point, bearing, distance float64) orb.Point { + aLat := deg2rad(p[1]) + aLon := deg2rad(p[0]) + + bearingRadians := deg2rad(bearing) + + distanceRatio := distance / orb.EarthRadius + bLat := math.Asin(math.Sin(aLat)*math.Cos(distanceRatio) + math.Cos(aLat)*math.Sin(distanceRatio)*math.Cos(bearingRadians)) + bLon := aLon + math.Atan2(math.Sin(bearingRadians)*math.Sin(distanceRatio)*math.Cos(aLat), math.Cos(distanceRatio)-math.Sin(aLat)*math.Sin(bLat)) + + r := orb.Point{ + rad2deg(bLon), + rad2deg(bLat), + } + + return r +} + +func PointAtDistanceAlongLine(ls orb.LineString, distance float64) (orb.Point, float64) { + if len(ls) == 0 { + panic("empty LineString") + } + + if distance < 0 || len(ls) == 1 { + return ls[0], 0.0 + } + + travelled := 0.0 + i := 1 + var from, to orb.Point + for ; i < len(ls); i++ { + from = ls[i-1] + to = ls[i] + actualSegmentDistance := DistanceHaversine(from, to) + expectedSegmentDistance := distance - travelled + if expectedSegmentDistance < actualSegmentDistance { + bearing := Bearing(from, to) + return PointAtBearingAndDistance(from, bearing, expectedSegmentDistance), bearing + } + travelled += actualSegmentDistance + } + + bearing := Bearing(from, to) + return to, bearing +} diff --git a/geo/distance_test.go b/geo/distance_test.go index af82a4b..f13933e 100644 --- a/geo/distance_test.go +++ b/geo/distance_test.go @@ -89,3 +89,112 @@ func TestMidpoint(t *testing.T) { t.Errorf("expected %v, got %v", answer, m) } } + +func TestPointAtBearingAndDistance(t *testing.T) { + expected := orb.Point{-0.841153, 52.68179432} + bearing := 127.373 + distance := 85194.89 + actual := PointAtBearingAndDistance(orb.Point{-1.8444, 53.1506}, bearing, distance) + + if d := DistanceHaversine(actual, expected); d > 1 { + t.Errorf("expected %v, got %v (%vm away)", expected, actual, d) + } +} +func TestMidpointAgainstPointAtBearingAndDistance(t *testing.T) { + a := orb.Point{-1.8444, 53.1506} + b := orb.Point{0.1406, 52.2047} + bearing := Bearing(a, b) + distance := DistanceHaversine(a, b) + acceptableTolerance := 1e-06 // unit is meters + + p1 := PointAtBearingAndDistance(a, bearing, distance/2) + p2 := Midpoint(a, b) + + if d := DistanceHaversine(p1, p2); d > acceptableTolerance { + t.Errorf("expected %v to be within %vm of %v", p1, acceptableTolerance, p2) + } +} + +func TestPointAtDistanceAlongLineWithEmptyLineString(t *testing.T) { + defer func() { + if r := recover(); r == nil { + t.Errorf("PointAtDistanceAlongLine did not panic") + } + }() + + line := orb.LineString{} + PointAtDistanceAlongLine(line, 90000) +} + +func TestPointAtDistanceAlongLineWithSinglePoint(t *testing.T) { + expectedPoint := orb.Point{-1.8444, 53.1506} + line := orb.LineString{ + expectedPoint, + } + actualPoint, actualBearing := PointAtDistanceAlongLine(line, 90000) + + if actualPoint != expectedPoint { + t.Errorf("expected %v but got %v", expectedPoint, actualPoint) + } + if actualBearing != 0.0 { + t.Errorf("expected %v but got %v", actualBearing, 0.0) + } +} + +func TestPointAtDistanceAlongLineWithMinimalPoints(t *testing.T) { + expected := orb.Point{-0.841153, 52.68179432} + acceptableDistanceTolerance := 1.0 // unit is meters + line := orb.LineString{ + orb.Point{-1.8444, 53.1506}, + orb.Point{0.1406, 52.2047}, + } + acceptableBearingTolerance := 0.01 // unit is degrees + expectedBearing := Bearing(line[0], line[1]) + actual, actualBearing := PointAtDistanceAlongLine(line, 85194.89) + + if d := DistanceHaversine(expected, actual); d > acceptableDistanceTolerance { + t.Errorf("expected %v to be within %vm of %v (%vm away)", actual, acceptableDistanceTolerance, expected, d) + } + if b := math.Abs(actualBearing - expectedBearing); b > acceptableBearingTolerance { + t.Errorf("expected bearing %v to be within %v degrees of %v", actualBearing, acceptableBearingTolerance, expectedBearing) + } +} + +func TestPointAtDistanceAlongLineWithMultiplePoints(t *testing.T) { + expected := orb.Point{-0.78526, 52.65506} + acceptableTolerance := 1.0 // unit is meters + line := orb.LineString{ + orb.Point{-1.8444, 53.1506}, + orb.Point{-0.8411, 52.6817}, + orb.Point{0.1406, 52.2047}, + } + acceptableBearingTolerance := 0.01 // unit is degrees + expectedBearing := Bearing(line[1], line[2]) + actualPoint, actualBearing := PointAtDistanceAlongLine(line, 90000) + + if d := DistanceHaversine(expected, actualPoint); d > acceptableTolerance { + t.Errorf("expected %v to be within %vm of %v (%vm away)", expected, acceptableTolerance, actualPoint, d) + } + if b := math.Abs(actualBearing - expectedBearing); b > acceptableBearingTolerance { + t.Errorf("expected bearing %v to be within %v degrees of %v", actualBearing, acceptableBearingTolerance, expectedBearing) + } +} + +func TestPointAtDistanceAlongLinePastEndOfLine(t *testing.T) { + expected := orb.Point{0.1406, 52.2047} + line := orb.LineString{ + orb.Point{-1.8444, 53.1506}, + orb.Point{-0.8411, 52.6817}, + expected, + } + acceptableBearingTolerance := 0.01 // unit is degrees + expectedBearing := Bearing(line[1], line[2]) + actualPoint, actualBearing := PointAtDistanceAlongLine(line, 200000) + + if actualPoint != expected { + t.Errorf("expected %v but got %v", expected, actualPoint) + } + if b := math.Abs(actualBearing - expectedBearing); b > acceptableBearingTolerance { + t.Errorf("expected bearing %v to be within %v degrees of %v", actualBearing, acceptableBearingTolerance, expectedBearing) + } +}