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:

enter image description here

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? — mkennedy, Mar 31 '14 at 17:16

Adam Matan · Accepted Answer · 2014-04-02T06:49:14.960

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#)

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

1 Answers1