How do I calculate the Azimuth (angle to north) between two WGS84 coordinates

I have got two WGS84 coordinates, latitude and longitude in degrees. These points are rather close together, e.g. only one metre apart.

Is there an easy way to calculate the azimuth of the line between these points, that is, the angle to north?

The naive approach would be to assume a Cartesian coordinate system (because these points are so close together) and just use

sin(a) = abs(L2-L1) / sqrt(sqr(L2-L1) + sqr(B2-B1))

a = azimuth L1, L2 = longitude B1, B2 = latitude

The error will be larger as the coordinates move away from the equator because there the distance between two longitudinal degrees becomes increasingly smaller than the one between two latitudinal degrees (which remains constant).

I found some quite complex formulas which I don't really want to implement because they seem to be overkill for points that are that close together and I don't need very high precision (two decimals are enough, one is probably fine either since there are other factors that reduce precision anyway, like the one the GPS returns).

Maybe I could just determine an approximate longitudinal correction factor depending on latitude and use somthing like this:

sin(a) = abs(L2*f-L1*f) / sqrt(sqr(L2*f-L1*f) + sqr(B2-B1))

where f is the correction factor

Any hints?

(I don't want to use any libraries for this, especially not ones that require runtime licenses. Any MPLed Delphi Source would be great.)

From stackoverflow

• I would recommend implementing a correction factor based on the longitude. I implemented a simular routine once to return all geocoded records within x miles of a specific spot and ran into simular issues. Unfortunately I don't have the code anymore, and can't seem to recall how I got to the correction number but you are on the right track.

• I found this link

  http://williams.best.vwh.net/avform.htm

  given in the answer to

  http://stackoverflow.com/questions/378281/lat-lon-distance-heading-lat-lon

  This looks promising, especially the flat earth approximation given near the end.

• The formulas that you refer to in the text are to calculate the great circle distance between 2 points. Here's how I calculate the angle between points:

  uses Math, ...;
  ...
  
  const
    cNO_ANGLE=-999;
  
  ...
  
  function getAngleBetweenPoints(X1,Y1,X2,Y2:double):double;
  var
    dx,dy:double;
  begin
    dx := X2 - X1;
    dy := Y2 - Y1;
  
    if (dx > 0) then  result := (Pi*0.5) - ArcTan(dy/dx)   else
    if (dx < 0) then  result := (Pi*1.5) - ArcTan(dy/dx)   else
    if (dy > 0) then  result := 0                          else
    if (dy < 0) then  result := Pi                         else
                      result := cNO_ANGLE; // the 2 points are equal
  
    result := RadToDeg(result);
  end;
  

  • Remember to handle the situation where 2 points are equal (check if the result equals cNO_ANGLE, or modify the function to throw an exception);

  • This function assumes that you're on a flat surface. With the small distances that you've mentioned this is all fine, but if you're going to be calculating the heading between cities all over the world you might want to look into something that takes the shape of the earth in count;

  • It's best to provide this function with coordinates that are already mapped to a flat surface. You could feed WGS84 Latitude directly into Y (and lon into X) to get a rough approximation though.

  Jader Dias : Great solution! I translated it to C# to replace my previous answer.

  Martin Beckett : Isn't that just atan2() ?

  MPelletier : @Martin Beckett: Yes, atan2 does this too with `Atan2(dy;dx)`, but of course a negative dx will return a negative value that will have to be added to 360. And 0's have to be dealt with separately.

• Here is the C# solution. Tested for 0, 45, 90, 135, 180, 225, 270 and 315 angles.

  Edit I replaced my previous ugly solution, by the C# translation of Wouter's solution:

  public double GetAzimuth(LatLng destination)
  {
      var longitudinalDifference = destination.Lng - this.Lng;
      var latitudinalDifference = destination.Lat - this.Lat;
      var azimuth = (Math.PI * .5d) - Math.Atan(latitudinalDifference / longitudinalDifference);
      if (longitudinalDifference > 0) return azimuth;
      else if (longitudinalDifference < 0) return azimuth + Math.PI;
      else if (latitudinalDifference < 0) return Math.PI;
      return 0d;
  }
  
  public double GetDegreesAzimuth(LatLng destination)
  {
      return RadiansToDegreesConversionFactor * GetAzimuth(destination);
  }
  

  Jonathan Leffler : Isn't it best to use Math.Atan2 rather than Math.Atan to avoid division, and in particular to avoid division by zero. The added advantage is it gives an answer with a range of -π .. +π, thus obviating the corrections you make at the end.

• This would work only for small differences. Otherwise you can't just "latitudinalDifference / longitudinalDifference".