GPS distance and initial bearing between points (MySQL)
First we create a view to help with the 3D geometry (6378 = Earth's radius in km):
1 2 CREATE VIEW gpsGlb AS 3 SELECT 4 LocID 5 ,6378 * COS(RADIANS(Latitude)) * COS(RADIANS(Longitude)) AS x 6 ,6378 * COS(RADIANS(Latitude)) * SIN(RADIANS(Longitude)) AS y 7 ,6378 * SIN(RADIANS(Latitude)) AS z 8 FROM station;
Here I query for distances to all my locations that are NOT LocID = 405 (rounded miles in my case):
1 2 SELECT 3 LocID 4 ,ROUND((2 * 6378 * ASIN(d / 2 / 6378)) * 0.621371192) AS dist_mi 5 FROM 6 (SELECT 7 SQRT(dx * dx + dy * dy + dz * dz) AS d 8 ,LocID 9 FROM 10 (SELECT 11 p1.x - p2.x AS dx 12 ,p1.y - p2.y AS dy 13 ,p1.z - p2.z AS dz 14 ,p2.LocID 15 FROM gpsGlb p1 16 JOIN gpsGlb p2 ON (p1.LocID = 405 AND p2.LocID != 405) 17 ) t1 18 ) t2 19 ORDER BY dist_mi
Here I get the initial bearing to the locations. The "boxed" calculation will come in handy later.
1 2 SELECT 3 LocID 4 ,(360 + DEGREES(ATAN2(y, x))) % 360 AS initBearing_deg 5 ,ROUND(((360 + DEGREES(ATAN2(y, x))) % 360) / 22.5) * 22.5 6 AS initBearingBoxed_deg 7 FROM 8 (SELECT 9 SIN(RADIANS(s2.Longitude - s1.Longitude)) * COS(RADIANS(s2.Latitude)) 10 AS y 11 ,COS(RADIANS(s1.Latitude)) * SIN(RADIANS(s2.Latitude)) 12 - SIN(RADIANS(s1.Latitude)) * COS(RADIANS(s2.Latitude)) 13 * COS(RADIANS(s2.Longitude - s1.Longitude)) 14 AS x 15 ,s2.LocID 16 FROM station s1 17 JOIN station s2 ON (s1.LocID = 405 AND s2.LocID != 405) 18 ) q1
Here's the combined query plus boxed degrees converted to 'NNE', etc. I've also added a limit for the distance in the qq1 subquery.
1 2 SELECT 3 qq2.LocID 4 ,dist_mi 5 ,CASE initBearingBoxed_deg 6 WHEN 22.5 THEN 'NNE' WHEN 45 THEN 'NE' 7 WHEN 67.5 THEN 'ENE' WHEN 90 THEN 'E' 8 WHEN 112.5 THEN 'ESE' WHEN 135 THEN 'SE' 9 WHEN 157.5 THEN 'SSE' WHEN 180 THEN 'S' 10 WHEN 202.5 THEN 'SSW' WHEN 225 THEN 'SW' 11 WHEN 247.5 THEN 'WSW' WHEN 270 THEN 'W' 12 WHEN 292.5 THEN 'WNW' WHEN 315 THEN 'NW' 13 WHEN 337.5 THEN 'NNW' ELSE 'N' 14 END AS bearing 15 FROM ( 16 SELECT 17 LocID 18 ,ROUND((2 * 6378 * ASIN(d / 2 / 6378)) * 0.621371192) AS dist_mi 19 FROM 20 (SELECT 21 SQRT(dx * dx + dy * dy + dz * dz) AS d 22 ,LocID 23 FROM 24 (SELECT 25 p1.x - p2.x AS dx 26 ,p1.y - p2.y AS dy 27 ,p1.z - p2.z AS dz 28 ,p2.LocID 29 FROM gpsGlb p1 30 JOIN gpsGlb p2 ON (p1.LocID = 405 AND p2.LocID != 405) 31 ) t1 32 ) t2 33 ) qq1 34 JOIN ( 35 SELECT 36 LocID 37 ,(360 + DEGREES(ATAN2(y, x))) % 360 AS initBearing_deg 38 ,(360 + ROUND((DEGREES(ATAN2(y, x))) / 22.5) * 22.5) % 360 39 AS initBearingBoxed_deg 40 FROM 41 (SELECT 42 SIN(RADIANS(s2.Longitude - s1.Longitude)) * COS(RADIANS(s2.Latitude)) 43 AS y 44 ,COS(RADIANS(s1.Latitude)) * SIN(RADIANS(s2.Latitude)) 45 - SIN(RADIANS(s1.Latitude)) * COS(RADIANS(s2.Latitude)) 46 * COS(RADIANS(s2.Longitude - s1.Longitude)) 47 AS x 48 ,s2.LocID 49 FROM station s1 50 JOIN station s2 ON (s1.LocID = 405 AND s2.LocID != 405) 51 ) q1 52 ) qq2 ON (qq1.LocID = qq2.LocID 53 AND qq1.dist_mi <= 60) 54 ORDER BY dist_mi