Matplotlib Basemap tutorial 08: Shooting Great Circles

Following a question from Ricardo Gama (see his comment), I’ve prepared this new tutorial. He wondered if Basemap has a function similar to the track1 function in matlab (you know, that crappy costly thing…)…

Here is what I obtained :

The goal:

Plotting great circles with Basemap, but knowing only the longitude, latitude, the azimuth and a distance. Only the origin point is known.

The process:

Instead of modifying (for now) the basemap package, I’ve translated a routine I found on the internet. Ed Williams wrote very good scripts for shooting great circles (see his website : http://williams.best.vwh.net/gccalc.htm). Note, I haven’t contacted Ed, so, please note this code is the property of Ed, I’ve just translated the script.

The process:

I defined a function call “shoot” that takes longitude, latitude, azimuth and the distance until which  the great circle needs to be drawn. The function returns the longitude and latitude of the second point, allowing you to call the regular m.drawgreatcircle()

The only thing to do is to call that function, get the second point’s lon/lat and plot the draw the great circle :

glon1 = 4.360515
glat1 = 50.79747
azimuth = 270.
maxdist = 20000

glon2, glat2, backazimuth = shoot(glon1, glat1, azimuth, maxdist)
m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,color='red', lw=4.)

Note, with this code, if you input a too long maxdist, the smallest great circle is plotted…

To plot the great circle until you cross the -180/+180° limit, one can set up a function like:

def great(m, startlon, startlat, azimuth,*args, **kwargs):
    glon1 = startlon
    glat1 = startlat
    glon2 = glon1
    glat2 = glat1

    step = 50

    glon2, glat2, baz = shoot(glon1, glat1, azimuth, step)
    if azimuth-180 >= 0:
        while glon2 <= startlon:
            m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,**kwargs)
            azimuth = baz + 180.
            glat1, glon1 = (glat2, glon2)

            glon2, glat2, baz = shoot(glon1, glat1, azimuth, step)
    elif azimuth-180 < 0:
        while glon2 >= startlon:
            m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,**kwargs)
            azimuth = baz + 180.
            glat1, glon1 = (glat2, glon2)

            glon2, glat2, baz = shoot(glon1, glat1, azimuth, step)

for which you only have to provide the starting lat/lon, the starting azimuth. This function expects the basemap (m) as first arg. It will compute the locations along the great circle every “step” distance.

The code is after this break:

#
# BaseMap example by geophysique.be
# tutorial 08
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

import numpy as np

### PARAMETERS FOR MATPLOTLIB :
import matplotlib as mpl
mpl.rcParams['font.size'] = 10.
mpl.rcParams['font.family'] = 'Comic Sans MS'
mpl.rcParams['axes.labelsize'] = 8.
mpl.rcParams['xtick.labelsize'] = 6.
mpl.rcParams['ytick.labelsize'] = 6.

def shoot(lon, lat, azimuth, maxdist=None):
    """Shooter Function
    Original javascript on http://williams.best.vwh.net/gccalc.htm
    Translated to python by Thomas Lecocq
    """
    glat1 = lat * np.pi / 180.
    glon1 = lon * np.pi / 180.
    s = maxdist / 1.852
    faz = azimuth * np.pi / 180.

    EPS= 0.00000000005
    if ((np.abs(np.cos(glat1))<EPS) and not (np.abs(np.sin(faz))<EPS)):
        alert("Only N-S courses are meaningful, starting at a pole!")

    a=6378.13/1.852
    f=1/298.257223563
    r = 1 - f
    tu = r * np.tan(glat1)
    sf = np.sin(faz)
    cf = np.cos(faz)
    if (cf==0):
        b=0.
    else:
        b=2. * np.arctan2 (tu, cf)

    cu = 1. / np.sqrt(1 + tu * tu)
    su = tu * cu
    sa = cu * sf
    c2a = 1 - sa * sa
    x = 1. + np.sqrt(1. + c2a * (1. / (r * r) - 1.))
    x = (x - 2.) / x
    c = 1. - x
    c = (x * x / 4. + 1.) / c
    d = (0.375 * x * x - 1.) * x
    tu = s / (r * a * c)
    y = tu
    c = y + 1
    while (np.abs (y - c) > EPS):

        sy = np.sin(y)
        cy = np.cos(y)
        cz = np.cos(b + y)
        e = 2. * cz * cz - 1.
        c = y
        x = e * cy
        y = e + e - 1.
        y = (((sy * sy * 4. - 3.) * y * cz * d / 6. + x) *
              d / 4. - cz) * sy * d + tu

    b = cu * cy * cf - su * sy
    c = r * np.sqrt(sa * sa + b * b)
    d = su * cy + cu * sy * cf
    glat2 = (np.arctan2(d, c) + np.pi) % (2*np.pi) - np.pi
    c = cu * cy - su * sy * cf
    x = np.arctan2(sy * sf, c)
    c = ((-3. * c2a + 4.) * f + 4.) * c2a * f / 16.
    d = ((e * cy * c + cz) * sy * c + y) * sa
    glon2 = ((glon1 + x - (1. - c) * d * f + np.pi) % (2*np.pi)) - np.pi	

    baz = (np.arctan2(sa, b) + np.pi) % (2 * np.pi)

    glon2 *= 180./np.pi
    glat2 *= 180./np.pi
    baz *= 180./np.pi

    return (glon2, glat2, baz)

def great(m, startlon, startlat, azimuth,*args, **kwargs):
    glon1 = startlon
    glat1 = startlat
    glon2 = glon1
    glat2 = glat1

    step = 50

    glon2, glat2, baz = shoot(glon1, glat1, azimuth, step)
    if azimuth-180 >= 0:
        while glon2 <= startlon:
            m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,**kwargs)
            azimuth = baz + 180.
            glat1, glon1 = (glat2, glon2)

            glon2, glat2, baz = shoot(glon1, glat1, azimuth, step)
    elif azimuth-180 < 0:
        while glon2 >= startlon:
            m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,**kwargs)
            azimuth = baz + 180.
            glat1, glon1 = (glat2, glon2)

            glon2, glat2, baz = shoot(glon1, glat1, azimuth, step)

fig = plt.figure(figsize=(11.7,8.3))

#Custom adjust of the subplots

plt.subplots_adjust(left=0.05,right=0.95,top=0.90,bottom=0.05,wspace=0.15,hspace=0.05)

ax = plt.subplot(111)

#Let's create a basemap of the world

m = Basemap(resolution='l')
m.drawcountries()
m.drawcoastlines()

glon1 = 4.360515
glat1 = 50.79747

azimuth = 270.
maxdist = 9000
glon2, glat2, baz = shoot(glon1, glat1, azimuth, maxdist)
m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,color='red', lw=4.)

azimuth = 90.
maxdist = 12000
glon2, glat2, baz = shoot(glon1, glat1, azimuth, maxdist)
m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,color='blue', lw=4.)

azimuth = 225.
maxdist = 15000
glon2, glat2, baz = shoot(glon1, glat1, azimuth, maxdist)
m.drawgreatcircle(glon1, glat1, glon2, glat2,del_s=50,color='green', lw=4.)

startlon = 4.360515
startlat = 50.79747
azimuth = 165.
great(m,startlon, startlat, azimuth, color='orange', lw=2.0)
great(m,startlon, startlat, azimuth+180., color='orange', lw=2.0)

plt.savefig('tutorial08.png',dpi=300)

plt.show()