Never been to DZone Snippets before?

Snippets is a public source code repository. Easily build up your personal collection of code snippets, categorize them with tags / keywords, and share them with the world

About this user

Steve Clay http://mrclay.org/

« Newer Snippets
Older Snippets »
Showing 1-1 of 1 total  RSS 

GPS distance and initial bearing between points (MySQL)

Assume you have a table of locations with Latitude and Longitude for each one. In my case the table is "station" and the primary key is "LocID".

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
« Newer Snippets
Older Snippets »
Showing 1-1 of 1 total  RSS