diff --git a/util/neighbour/nearest.go b/util/neighbour/nearest.go new file mode 100644 index 000000000..dbc3c82ee --- /dev/null +++ b/util/neighbour/nearest.go @@ -0,0 +1,71 @@ +// Package neighbour contains a basic non-optimized nearest-neighbour +// algorithm implementation for terrestrial GPS coordinates. +package neighbour + +import ( + "math" + + "golang.org/x/exp/slices" +) + +// Location is a latitude/longitude pair representing a point on Earth. +type Location struct { + Latitude float64 + Longitude float64 +} + +// Distance calculates the great-circle distance between two points using the +// Haversine formula, in kilometers. +// +// This is also known as the "as-the-crow-flies" distance. +func (l Location) Distance(other Location) float64 { + // For the following variable definitions: + // φ is latitude ("phi") + // λ is longitude ("lambda") + // R is earth’s radius (mean radius = 6,371km) + // + // We can calculate the distance using the haversine formula as such: + // a = sin²((φB - φA)/2) + cos φA * cos φB * sin²((λB - λA)/2) + // c = 2 * atan2( √a, √(1−a) ) + // d = R * c + + // Convert our latitude/longitude to radians, since the various math + // functions take radians but latitude/longitude are in degrees. + lat1, lon1 := degreesToRadians(l.Latitude), degreesToRadians(l.Longitude) + lat2, lon2 := degreesToRadians(other.Latitude), degreesToRadians(other.Longitude) + + deltaPhi := lat2 - lat1 + deltaLambda := lon2 - lon1 + + // Haversine + a := math.Pow(math.Sin(deltaPhi/2), 2) + + math.Cos(lat1)*math.Cos(lat2)*math.Pow(math.Sin(deltaLambda/2), 2) + c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a)) + + // Return the distance in km. + const earthRadiusKm = 6371 + return c * earthRadiusKm +} + +func degreesToRadians(d float64) float64 { + return d * math.Pi / 180 +} + +// Neighbours returns the nearest n neighbours to the provided point from the +// set of candidates. +func Neighbours(point Location, n int, candidates []Location) []Location { + // Calculate all distances up-front to avoid recalculating during the + // sort below. + distances := map[Location]float64{} + for _, candidate := range candidates { + distances[candidate] = candidate.Distance(point) + } + + // Sort the candidates slice by their distance to the provided point. + candidates = slices.Clone(candidates) + slices.SortFunc(candidates, func(a, b Location) bool { + return distances[a] < distances[b] + }) + + return candidates[:n] +} diff --git a/util/neighbour/nearest_test.go b/util/neighbour/nearest_test.go new file mode 100644 index 000000000..e550f8a9b --- /dev/null +++ b/util/neighbour/nearest_test.go @@ -0,0 +1,51 @@ +package neighbour + +import ( + "math" + "reflect" + "testing" +) + +func TestHaversine(t *testing.T) { + one := Location{51.510357, -0.116773} // King's College, London + two := Location{38.889931, -77.009003} // The White House + + dist := one.Distance(two) + want := 5897.658 + + if math.Abs(want-dist) > 0.001 { + t.Fatalf("distance mismatch; got %v, want %v", dist, want) + } +} + +func TestNeighbours(t *testing.T) { + // Provincial capitals + capitals := []Location{ + {48.4283182, -123.3649533}, // Victoria, BC, Canada + {60.721571, -135.054932}, // Whitehorse, YT, Canada + {53.5462055, -113.491241}, // Edmonton, AB, Canada + {62.4540807, -114.377385}, // Yellowknife, NT, Canada + {50.44876, -104.61731}, // Regina, SK, Canada + {49.8955367, -97.1384584}, // Winnipeg, MB, Canada + {63.74944, -68.521857}, // Iqaluit, NU, Canada + {43.6534817, -79.3839347}, // Toronto, ON, Canada + {45.5031824, -73.5698065}, // Montreal, QC, Canada + {45.94780155, -66.6534707}, // Fredericton, NB, Canada + {44.648618, -63.5859487}, // Halifax, NS, Canada + {46.234953, -63.132935}, // Charlottetown, PE, Canada + {47.5614705, -52.7126162}, // St. John’s, NL, Canada + } + + // Thunder Bay, Ontario, Canada + point := Location{48.382221, -89.246109} + nearest := Neighbours(point, 4, capitals) + want := []Location{ + {49.8955367, -97.1384584}, // Winnipeg, MB, Canada + {43.6534817, -79.3839347}, // Toronto, ON, Canada + {50.44876, -104.61731}, // Regina, SK, Canada + {45.5031824, -73.5698065}, // Montreal, QC, Canada + } + if !reflect.DeepEqual(nearest, want) { + t.Errorf("nearest points mismatch\ngot: %v\nwant: %v", nearest, want) + } +}