How do I find the lat/long that is x km north of a given lat/long

Question

I have some C# code that generates google maps. This codes looks at all the Points I need to plot on the map and then works out the Bounds of a rectangle to include those points. It then passes this bounds to the Google Maps API to set the zoom level appropriately to show all of the points on the map.

This code is working fine however I have a new requirement.

One of the points may have a precision associated with it. If this is the case then I draw a circle around the point with the radius set to the precision value. Again this works fine however my bounds checking is now not doing what I want it to do. I want to have the bounding box include the complete circle.

This requires an algorithm to take a point x and calculate the point y that would be z metres north of x and also z metres south of x.

Does anyone have this algorithm, preferably in C#. I did find a generic algorithm here but I appear to have not implemented this correctly as the answers I am getting are 1000s of km adrift.

This is the Generic example

Lat/lon given radial and distance

A point {lat,lon} is a distance d out on the tc radial from point 1 if:

     lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
     IF (cos(lat)=0)
        lon=lon1      // endpoint a pole
     ELSE
        lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
     ENDIF

And this is my C# translation.

  // Extend a Point North/South by the specified distance
    public static Point ExtendPoint(Point _pt, int _distance, int _bearing )
    {
        Decimal lat = 0.0;
        Decimal lng = 0.0;

        lat = Math.Asin(Math.Sin(_pt.Lat) * Math.Cos(_distance) + Math.Cos(_pt.Lat) * 
            Math.Sin(_distance) * Math.Cos(_bearing));

         if (Math.Cos(lat) == 0)
         {
            lng = _pt.Lng;      // endpoint a pole
         }
         else 
         {
             lng = (
                 (_pt.Lng - Math.Asin(Math.Sin(_bearing) * Math.Sin(_distance) / Math.Cos(lat)) 
                 + Math.PI) % (2 * Math.PI)) - Math.PI;
         }

         ret = new Point(lat,lng);
         return ret;
    }

I am calling this function with a bearing of 0 to calculate the new northerly position and a value of 180 to calculate the new southerly position.

Can anyone either see what I have done wrong or perhaps provide a known working algorithm.

Answer 1

If you have a given latitude and longitude you can calculate the correct latitude and longitude of an x-km change in latitude like so:

new-lat = ((old-km-north + x-km-change)/40,075) * 360)
           ^ is the ratio of the                  ^ times the ratio of the circle
           of the earth the change                by 360 to get the total ratio 
           covers.                                covered in degrees.

The same can apply to longitude. If you have the total distance plus the change you can calculate the total degrees in a similar fashion.

new-long = ((old-km-east + x-km-change)/40,075) * 360)
           ^ is the ratio of the                  ^ times the ratio of the circle
           of the earth the change                by 360 to get the total ratio 
           covers.                                covered in degrees.

Again, these calculations should work, but I'm running off pure intuition here, but the logic does seem to hold true.

Edit: As pointed out by Skizz 40,075 needs to be adjusted to the circumference of the earth at any given latitude using 2.pi.r.cos(lat) or 40074.cos(lat)

Answer 2

It is more accurate if you first reproject it to UTM and then check the distance.

Hope this helps

Answer 3

I have a very similar piece of code. It got me very close results when compared to another implementation.

I think the problem with yours is that you are using "distance" as linear distance in meters instead of angular distance in radians.

/// <summary>
/// Calculates the end-point from a given source at a given range (meters) and bearing (degrees).
/// This methods uses simple geometry equations to calculate the end-point.
/// </summary>
/// <param name="source">Point of origin</param>
/// <param name="range">Range in meters</param>
/// <param name="bearing">Bearing in degrees</param>
/// <returns>End-point from the source given the desired range and bearing.</returns>
public static LatLonAlt CalculateDerivedPosition(LatLonAlt source, double range, double bearing)
{
    double latA = source.Latitude * UnitConstants.DegreesToRadians;
    double lonA = source.Longitude * UnitConstants.DegreesToRadians;
    double angularDistance = range / GeospatialConstants.EarthRadius;
    double trueCourse = bearing * UnitConstants.DegreesToRadians;

    double lat = Math.Asin(
        Math.Sin(latA) * Math.Cos(angularDistance) + 
        Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse));

    double dlon = Math.Atan2(
        Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA), 
        Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat));

    double lon = ((lonA + dlon + Math.PI) % UnitConstants.TwoPi) - Math.PI;

    return new LatLonAlt(
        lat * UnitConstants.RadiansToDegrees, 
        lon * UnitConstants.RadiansToDegrees, 
        source.Altitude);
}

Where

public const double EarthRadius = 6378137.0;   //  WGS-84 ellipsoid parameters

and LatLonAlt is in degrees/meters (conversion takes place internally). Adjust as needed.

I assume you can figure out what the value for UnitConstants.DegreesToRadians is :)

Answer 4

I'm not sure if I'm missing something here, but I think the question could be rephrased as, "I have a lat/lon point, and I want to find the point x meters north and x meters south of that point."

If that's the question then you don't need to find a new longitude (which makes things simpler), you just need a new latitude. A degree of latitude is roughly 60 nautical miles long anywhere on Earth, and a nautical mile is 1,852 meters. So, for new latitudes x meters north and south:

north_lat = lat + x / (1852 * 60)
north_lat = min(north_lat, 90)

south_lat = lat - x / (1852 * 60)
south_lat = max(south_lat, -90)

This is not completely accurate because the Earth is not a perfect sphere with exactly 60 nautical miles between each degree of latitude. However, the other answers assume that lines of latitude are equidistant, so I'm assuming you don't care about that. If you're interested in how much error that might introduce, there is a nice table on Wikipedia that shows "Surface distance per 1° change in latitude" for different latitudes at this link:

http://en.wikipedia.org/wiki/Latitude#Degree_length