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

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()```

## 5 thoughts on “Matplotlib Basemap tutorial 08: Shooting Great Circles”

1. Ricardo Gama says:

Dear Thomas,
This tutorial is great and very helpful. Thank you.
Best regards,
Ricardo

2. Juan Fernandez says:

Hi!

Nice tutorials you have for Basemap.

By the way. I am trying to draw some great circles that cross longitude 180 (or -180, as you want to see it). So my map has the borders at longitudes -220 an -45 and when I try to draw a great circle from west of longitude 180 to west of longitude 180, the great circles are not great circles but strange lines. Any idea of how to fix this?

Thanx!

Juan

3. Juan Fernandez says:

sorry I mean from west of longitude -180 to east of longitude -180…