Distance and bearings between GPS coordinates in SQL

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.

3 thoughts on “Distance and bearings between GPS coordinates in SQL

  1. says:

    Really great post!

    I’d been pondering some of the concepts demonsrated here for a ship routing idea I had a while back; I wanted to programmatically determine the shortest path between two seaports, given their Lat/Lon. The hard part is figuring out waypoints to get around obstacles. I might be able to use common shipping lanes like sortcuts or Zelda-esque warp points. If I can get a route (or set of feasible routes), then I can weight them by the ocean’s current vector and have a really slick routing system.

    Was this post a part of something you did for work? Planning the shortest path for cattle drives or something?

  2. Robinhood says:

    Can u post the vb6 codes for above to calculate Bearing between two coordinates. your reply highly appreciated.Thank you.

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.