Another post from my sql abominations series. I’m running this on MySQL 5 particularly.
Assume you have a table of locations with Latitude and Longitude for each one. In my case the table is “station”, primary key being “LocID”. With help from this article, first we create a view to get 3D coordinates (6378 = Earth’s radius in km):
CREATE VIEW gpsGlb AS
SELECT
LocID
,6378 * COS(RADIANS(Latitude)) * COS(RADIANS(Longitude)) AS x
,6378 * COS(RADIANS(Latitude)) * SIN(RADIANS(Longitude)) AS y
,6378 * SIN(RADIANS(Latitude)) AS z
FROM station;
Now we can query for great-circle distance (I want rounded miles) to all locations from, say, LocID 405:
SELECT
LocID
,ROUND((2 * 6378 * ASIN(d / 2 / 6378)) * 0.621371192) AS dist_mi
FROM
(SELECT
SQRT(dx * dx + dy * dy + dz * dz) AS d
,LocID
FROM
(SELECT
p1.x - p2.x AS dx
,p1.y - p2.y AS dy
,p1.z - p2.z AS dz
,p2.LocID
FROM gpsGlb p1
JOIN gpsGlb p2 ON (p1.LocID = 405 AND p2.LocID != 405)
) t1
) t2
ORDER BY dist_mi
With help from this article, we can query for the initial bearing to each location. The “boxed” calculation will come in handy later.
SELECT
LocID
,(360 + DEGREES(ATAN2(y, x))) % 360 AS initBearing_deg
,ROUND(((360 + DEGREES(ATAN2(y, x))) % 360) / 22.5) * 22.5
AS initBearingBoxed_deg
FROM
(SELECT
SIN(RADIANS(s2.Longitude - s1.Longitude)) * COS(RADIANS(s2.Latitude))
AS y
,COS(RADIANS(s1.Latitude)) * SIN(RADIANS(s2.Latitude))
- SIN(RADIANS(s1.Latitude)) * COS(RADIANS(s2.Latitude))
* COS(RADIANS(s2.Longitude - s1.Longitude))
AS x
,s2.LocID
FROM station s1
JOIN station s2 ON (s1.LocID = 405 AND s2.LocID != 405)
) q1
What you’ve all been waiting for! The combined query plus boxed compass directions (like ‘NNE’), etc. I’ve also added a limit for the distance in the qq1 subquery since I only want close locations.
SELECT
qq2.LocID
,dist_mi
,CASE initBearingBoxed_deg
WHEN 22.5 THEN 'NNE' WHEN 45 THEN 'NE'
WHEN 67.5 THEN 'ENE' WHEN 90 THEN 'E'
WHEN 112.5 THEN 'ESE' WHEN 135 THEN 'SE'
WHEN 157.5 THEN 'SSE' WHEN 180 THEN 'S'
WHEN 202.5 THEN 'SSW' WHEN 225 THEN 'SW'
WHEN 247.5 THEN 'WSW' WHEN 270 THEN 'W'
WHEN 292.5 THEN 'WNW' WHEN 315 THEN 'NW'
WHEN 337.5 THEN 'NNW' ELSE 'N'
END AS bearing
FROM (
SELECT
LocID
,ROUND((2 * 6378 * ASIN(d / 2 / 6378)) * 0.621371192) AS dist_mi
FROM
(SELECT
SQRT(dx * dx + dy * dy + dz * dz) AS d
,LocID
FROM
(SELECT
p1.x - p2.x AS dx
,p1.y - p2.y AS dy
,p1.z - p2.z AS dz
,p2.LocID
FROM gpsGlb p1
JOIN gpsGlb p2 ON (p1.LocID = 405 AND p2.LocID != 405)
) t1
) t2
) qq1
JOIN (
SELECT
LocID
,(360 + DEGREES(ATAN2(y, x))) % 360 AS initBearing_deg
,(360 + ROUND((DEGREES(ATAN2(y, x))) / 22.5) * 22.5) % 360
AS initBearingBoxed_deg
FROM
(SELECT
SIN(RADIANS(s2.Longitude - s1.Longitude)) * COS(RADIANS(s2.Latitude))
AS y
,COS(RADIANS(s1.Latitude)) * SIN(RADIANS(s2.Latitude))
- SIN(RADIANS(s1.Latitude)) * COS(RADIANS(s2.Latitude))
* COS(RADIANS(s2.Longitude - s1.Longitude))
AS x
,s2.LocID
FROM station s1
JOIN station s2 ON (s1.LocID = 405 AND s2.LocID != 405)
) q1
) qq2 ON (qq1.LocID = qq2.LocID
AND qq1.dist_mi <= 60)
ORDER BY dist_mi
Result set is something like:
LocID |
dist_mi |
bearing |
250 |
25 |
E |
260 |
30 |
NNE |
240 |
42 |
ENE |
Hope this is useful to someone.