Many of the functions I'm working on consume arrays of points (GPS tracks) []float64{lng,lat} on which to run statistical analysis. Individual tracks can have 50,000+ points, describing a real journey from A to B.
Testing functions that process GPS tracks has been unexpectedly difficult. Test data of the form [1.0,2.0], to test logically is fine. But beyond that, I want to be able to test consistency in things like finding clusters, or coefficient breakpoints.
For testing I'm not concerned with where the locations are on earth, but it's handy to be able to view tracks on a map for visual confirmation. So the coordinates need to, sort of, coordinate.
I've created a function that generates semi-coherent GPS location data tracks. It's trivial to create a 5,000 data input track, designed to test something in particular. Being able to export it as geoJSON and view the shape on map is good for speedy intuition checks.
I can setup the data to be skewed one way or another, or insert a cluster of weirdness. Particularly useful at larger data sizes, where one analysis method may be better than another.
func CreateGPSData() {
line := orb.LineString{orb.Point{-3.188267, 55.953251}}
// general direction and distance per point
bearingRange := [2]int{Direction_SSE, Direction_SSW}
distanceRange := [2]int{10 * 1000, 15 * 1000}
// generate a test GPS track
for range 5000 {
newPoint := generateNewLocation(line[len(line)-1],
bearingRange,
distanceRange)
line = append(line, newPoint)
}
// add skewness in the data
bearingRange = [2]int{Direction_W, Direction_WNW}
distanceRange = [2]int{1000, 1500}
for range 100 {
newPoint := generateNewLocation(line[len(line)-1],
bearingRange,
distanceRange)
line = append(line, newPoint)
}
// do testing
}
What doesn't seem to be too common is generating a new point, some distance away in a particular direction. The opposite of point.DistanceTo(point2). That's really all that the following function does.
The breakdown starts with some readability. Using compass rose names really helps with reading the calling function.
// the compass rose, naming format for readability
const (
Direction_N = iota
Direction_NNE
Direction_NE
Direction_ENE
Direction_E
Direction_ESE
Direction_SE
Direction_SSE
Direction_S
Direction_SSW
Direction_SW
Direction_WSW
Direction_W
Direction_WNW
Direction_NW
Direction_NNW
)
const (
compassRoseDegrees = 22.5
)
Rather than simply a series of linear points, I want a GPS track that wiggles a little. So the parameters are a range of direction and distance, giving a small amount of randomness when calculating the next point.
// generateNewLocation returns a new point in the range of direction and
// distance. It is meant to build non-repetitive but predictable GPS tracks, to
// help generate test input cases.
//
// It's also meant to be readable code.
func generateNewLocation(start orb.Point, direction [2]int, distance [2]int) orb.Point {
Possibly the most common stupid mistakes I make are mixing up longitude and latitude indexes through typos. Named index values solve this.
// Mistakes with lon/lat indexing area easy to make, explicit index names
// helps
const (
Longitude = 0
Latitude = 1
)
I use the excellent orb package to hide much of the GPS specific calculations.
A point to note is that the Go maths package works in radians, so here's a helper function.
var (
latitudeOneDegreeOfDistance = 111000 // metres
newPoint orb.Point // []float64{Long, Lat}
// convert from degrees to radians
deg2rad = func(d float64) float64 { return d * math.Pi / 180 }
)
A while ago I went back to study maths, and have found lots of use for it. Here using trig to workout the distances on the ground.
// Use trigonometry of a right angled triangle to solve the distances on the ground.
// The hypotenuse is our desired distance to travel, and one angle
// is our desired bearing.
//
// now work out the vertical (longitude) and horizontal (latitude) sides in
// distance units.
hyp := (rand.Float64() * float64(distance[1]-distance[0])) + float64(distance[0])
// Get the compass bearing in degrees, with a little randomness between the
// general direction. Non-linear tracks are easier to troubleshoot visually.
bearingMin := float64(direction[0]) * 22.5
bearingMax := float64(direction[1]) * 22.5
angle := (rand.Float64() * (bearingMax - bearingMin)) + bearingMin
// Calulate the other side lengths using SOH CAH TOA. The Go math package
// works in radians
adj := math.Cos(deg2rad(angle)) * hyp // adjacent side of angle
opp := math.Sin(deg2rad(angle)) * hyp // opposite side of angle
Each degree change in latitude moves a fixed distance on the earth surface. So it is fairly simple to find the degree change we need to move our required distance on the ground.
// Each degree change in every latitude equates to ~111 km on the ground. So
// now find the degree change required for the length of adj
latitudeDelta := (1.0 / float64(latitudeOneDegreeOfDistance)) * adj
newPoint[Latitude] = start[Latitude] + latitudeDelta
Longitude 'distance moved on the earth surface' on the other hand changes depending on the corresponding latitude.
First need to find what one longitude degree of travel is, at the current latitude. This is the longitude equivalent to latitudeOneDegreeOfDistance = 111000.
// Distance on the ground for each degree of longitude changes depending on
// latitude because the earth is not perfectly spherical. So we need to
// calculate the distance of one degree longitude for our current latitude.
p1 := orb.Point{1.0, start[Latitude]}
p2 := orb.Point{2.0, start[Latitude]}
longitudeOneDegreeOfDistance := geo.Distance(p1, p2) // returns metres
After which it's the same calculation to find the longitude degree required to move our required distance on the ground.
// Now we can use this value to calculate the longitude degree change
// required to move opp distance (in a horizontal straight line) at this
// latitude.
longitudeDelta := (1.0 / longitudeOneDegreeOfDistance) * opp
Now we have the new point, located roughly at the direction and distance required.
// The new point is a vertical and horizontal shift to arrive at hyp
// distance from the start point on the required bearing.
newPoint[Longitude] = start[Longitude] + longitudeDelta
return newPoint
To print the geoJSON object, which can be viewed using the Mapbox viewer. The geojson package is part of orb.
fc := geojson.NewFeatureCollection()
f := geojson.NewFeature(line)
fc.Append(f)
rawJSON, _ := fc.MarshalJSON()
fmt.Println(string(rawJSON))
The geoJSON output, which can be dropped right into the Mapbox viewer.
{"features":[{"type":"Feature","geometry":{"type":"LineString","coordinates":[[-3.188267,55.953251],[-3.135746915373942,55.86790144687853],[-3.1203515713600174,55.77480489552086],[-3.1047549460145367,55.67255911652837],[-3.157869557083529,55.55060126862168],[-3.228498091310353,55.43475670091365],[-3.2496225715462885,55.436812333443044],[-3.2661004391340707,55.439707415244094],[-3.2813563020992547,55.44313063594729],[-3.3021255075874594,55.44562577470332],[-3.319279338837285,55.446998922512265]]},"properties":{"stroke":"#b464d6","stroke-width":3}}],"type":"FeatureCollection"}
The full code is in a Gist here.
Top comments (0)