A basic example of creating a geotrajectory visualization in KML, with a great circle drawn between two points of interest. Requires geodesic.py, a script that provides a function "tweensegs."

geodesic.py

import math

def tweensegs(longitude1, latitude1, longitude2, latitude2, num_of_segments):
    ptlon1 = longitude1
    ptlat1 = latitude1
    ptlon2 = longitude2
    ptlat2 = latitude2

    numberofsegments = num_of_segments
    onelessthansegments = numberofsegments - 1
    fractionalincrement = (1.0/onelessthansegments)

    ptlon1_radians = math.radians(ptlon1)
    ptlat1_radians = math.radians(ptlat1)
    ptlon2_radians = math.radians(ptlon2)
    ptlat2_radians = math.radians(ptlat2)

    distance_radians=2*math.asin(math.sqrt(math.pow((math.sin((ptlat1_radians-ptlat2_radians)/2)), 2) + math.cos(ptlat1_radians)*math.cos(ptlat2_radians)*math.pow((math.sin((ptlon1_radians-ptlon2_radians)/2)), 2)))
    # 6371.009 represents the mean radius of the earth
    #   shortest path distance
    distance_km = 6371.009 * distance_radians

    mylats = []
    mylons = []

    # Write the starting coordinates
    mylats.append([])
    mylons.append([])
    mylats[0] = ptlat1
    mylons[0] = ptlon1 

    f = fractionalincrement
    icounter = 1
    while (icounter <  onelessthansegments):
            icountmin1 = icounter - 1
            mylats.append([])
            mylons.append([])
            # f is expressed as a fraction along the route from point 1 to point 2
            A=math.sin((1-f)*distance_radians)/math.sin(distance_radians)
            B=math.sin(f*distance_radians)/math.sin(distance_radians)
            x = A*math.cos(ptlat1_radians)*math.cos(ptlon1_radians) + B*math.cos(ptlat2_radians)*math.cos(ptlon2_radians)
            y = A*math.cos(ptlat1_radians)*math.sin(ptlon1_radians) +  B*math.cos(ptlat2_radians)*math.sin(ptlon2_radians)
            z = A*math.sin(ptlat1_radians) + B*math.sin(ptlat2_radians)
            newlat=math.atan2(z, math.sqrt(math.pow(x, 2)+math.pow(y, 2)))
            newlon=math.atan2(y, x)
            newlat_degrees = math.degrees(newlat)
            newlon_degrees = math.degrees(newlon)
            mylats[icounter] = newlat_degrees
            mylons[icounter] = newlon_degrees
            icounter += 1
            f = f + fractionalincrement

    # Write the ending coordinates
    mylats.append([])
    mylons.append([])
    mylats[onelessthansegments] = ptlat2
    mylons[onelessthansegments] = ptlon2

    # Now,  the array mylats[] and mylons[] have the coordinate pairs for intermediate points along the geodesic
    #   My mylat[0], mylat[0] and mylat[num_of_segments-1], mylat[num_of_segments-1] are the geodesic end points
    buff = ''

    # Write a kml of the results
    zipcounter = 0
    kmlheader = "<?xml version=\"1.0\" encoding=\"UTF-8\"?><kml xmlns=\"http://www.opengis.net/kml/2.2\"><Document><name>LineString.kml</name><open>1</open><Placemark><name>unextruded</name><LineString><extrude>1</extrude><tessellate>1</tessellate><coordinates>"
    buff += kmlheader
    while (zipcounter < numberofsegments):
            outputstuff = repr(mylons[zipcounter]) + "," + repr(mylats[zipcounter]) + ",0 "
            buff += outputstuff
            zipcounter += 1
    kmlfooter = "</coordinates></LineString></Placemark></Document></kml>"
    buff += kmlfooter

    return buff

geotrajectory.py

import ogr, pykml
from pykml.parser import fromstring as from_kml_string
from pykml.factory import KML_ElementMaker as KML
from lxml import etree
from geodesic import tweensegs
# For geodesic.py and the function "tweensegs" see this article:
# http://gis.stackexchange.com/questions/47/what-tools-in-python-are-available-for-doing-great-circle-distance-line-creati/316#316

# These coordinates must be strings to keep precision
PARIS = ['2.35', '48.86', 'Paris, France']
ROME = ['12.50', '41.90', 'Rome, Italy']
NYC = ['-73.94', '40.66', 'New York City, U.S.A.']
JERUSALEM = ['35.22', '31.78', 'Jerusalem']
PT_A = PARIS
PT_B = JERUSALEM

pt_a = ogr.CreateGeometryFromWkt('POINT(%s %s)' % tuple(PT_A[0:2]))
pt_b = ogr.CreateGeometryFromWkt('POINT(%s %s)' % tuple(PT_B[0:2]))

# Buffer the points, an ellipse with radius 0.1 degrees, with 20 points per arc
buff_a = pt_a.Buffer(10.0, 20)
buff_b = pt_b.Buffer(15.0, 20)

features = [
    buff_a,
    buff_b
]

placemarks = [
    KML.Placemark(
        KML.name(PT_A[2]),
        KML.Point(
            KML.coordinates('%s, %s, 0' % tuple(PT_A[0:2]))
        ),
        KML.Style(
            KML.IconStyle(
                KML.Icon(
                    KML.href('http://maps.google.com/mapfiles/kml/paddle/A.png')
                ),
                KML.scale(4)
            )
        )
    ),
    KML.Placemark(
        KML.name(PT_B[2]),
        KML.Point(
            KML.coordinates('%s, %s, 0' % tuple(PT_B[0:2]))
        ),
        KML.Style(
            KML.IconStyle(
                KML.Icon(
                    KML.href('http://maps.google.com/mapfiles/kml/paddle/B.png')
                ),
                KML.scale(4)
            )
        )
    ),
    KML.Placemark(
        KML.Style(
            KML.LineStyle(
                KML.color('ffffffff'),
                KML.width(3)
            )
        ),
        KML.LineString(
            KML.coordinates(str(from_kml_string(tweensegs(float(PT_A[0]), float(PT_A[1]), float(PT_B[0]), float(PT_B[1]),
                20)).xpath('kml:Document/kml:Placemark/kml:LineString/kml:coordinates',
                namespaces={'kml': "http://www.opengis.net/kml/2.2"})[0]).strip())
        )
    )
]

for feature in features:
    # Parse each feature as a KML string and append the object to a sequence
    placemarks.append(KML.Placemark(
        KML.Style(
            KML.LineStyle(
                KML.color('ff0000ff'),
                KML.width(3)
            ),
            KML.PolyStyle(
                KML.color('770000ff')
            )
        ),
        from_kml_string(feature.ExportToKML())
    ))
        
document = KML.Document(*placemarks)

with open('GeoTrajectory.kml', 'wb') as kml_doc:
    kml_doc.write(etree.tostring(document, pretty_print=True))