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

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,
}
``````

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.

• Check out first answer here. And, you mean "I am trying to calculate ray tracing," right? Mar 31, 2014 at 17:16

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.