Python: Ray-trace a WGS84 point with heading and distance

Question

One of my services defines geofences for a mobile app. The geofences consists of a lat/lon pair and radius in meters, in JSON, e.g.:

{
    "geofenceId":"geofence_b39da38cff7be4b4693cbc31b2394bbb",
    "latitude":32.083611,
    "longitude":34.806111,
    "radius":15000.0
}

I would like to display these geofences in KML. Alas, KML files don't feature circles, so I have to approximate using a polygon with n vertices. Each vertex location is the head of a vector originating at the centre of the circle, whose length is the circle radius and has a given heading:

What's the Pythonic way to calculate P2 from P1, θ and radius?

Update: I've found some obvious solutions (like this), which assume that the earth is a perfect sphere. I calculated a circle with a radius of 15 Kilometers and 72 points. For each point, I've calculated the exact vincenty distance from the center. The minimum and maximum are:

14941.7356806
15014.5884245

The error is significant - a few dozen meters. I would like to find a solution that uses a better approximation.

Second Update:

I am trying to calculate ray tracing, which is (P1_lat, P1_lon, radius, heading) -> (P2_lat, P2_lon). This is different from the distance between two given points, which can be easily done with geopy's vincenty formula.

Answer 1

The C++ library geographiclib has been implemented in Python. Ray tracing, or more accurately - the Geodesic Direct problem, is as easy as:

geodesic.Geodesic.WGS84.Direct(lat1, lon1, azimuth_degrees, distance_meters)

Let's test the accuracy:

from geographiclib import geodesic
from geopy.distance import vincenty

lat1, lon1 = 32.074322, 34.792081         # Azrieli Centre, Tel Aviv

distances=[]
for degree in range(360):
    result=geodesic.Geodesic.WGS84.Direct(lat1, lon1, degree, 15000)
    lat2, lon2 = result["lat2"], result["lon2"]
    distances.append(vincenty( (lat1, lon1), (lat2, lon2)).meters)

print min(distances), max(distances)

Gives:

14999.9999876 14999.9999999

The accuracy is between 5 and 7 digits, which is great for my needs.

Links:

• geographiclib in PyPi
• A list of geographiclib implementations (C, Fortran, Java, Python, JavaScript, Matlab, IDL and C#)