# Matplotlib Basemap tutorial 09: Drawing circles

In the previous tutorial, I defined a “shoot” method to compute the landing point of a shoot from one point, to a given azimuth and distance. Using this logic, it’s possible to find the points situated at a given distance from a “centre” point, a circle.

The goal:

Drawing circles of a given radius around any point on earth.

The process:

We call the “shoot” method 360 times, one time per azimuth. We store the obtained locations in an array, that we plot. The “equi” method is defined as :

```def equi(m, centerlon, centerlat, radius, *args, **kwargs):
glon1 = centerlon
glat1 = centerlat
X = []
Y = []
for azimuth in range(0, 360):
glon2, glat2, baz = shoot(glon1, glat1, azimuth, radius)
X.append(glon2)
Y.append(glat2)
X.append(X[0])
Y.append(Y[0])

#m.plot(X,Y,**kwargs) #Should work, but doesn't...
X,Y = m(X,Y)
plt.plot(X,Y,**kwargs)```

In the image above, we called the method in a loop:

```radii = [500,1000,2000,3000,4000]

# Set number 1:
centerlon = 4.360515
centerlat = 50.79747

Fairly easy, isn’t it ?

Note : I’ve changed the projection type in this example, which led to a small bug: calling m.plot(X,Y) didn’t work, I have to convert X,Y using m(X,Y) …

The code is after this break:

```#
# BaseMap example by geophysique.be
# tutorial 09
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 equi(m, centerlon, centerlat, radius, *args, **kwargs):
glon1 = centerlon
glat1 = centerlat
X = []
Y = []
for azimuth in range(0, 360):
glon2, glat2, baz = shoot(glon1, glat1, azimuth, radius)
X.append(glon2)
Y.append(glat2)
X.append(X[0])
Y.append(Y[0])

#~ m.plot(X,Y,**kwargs) #Should work, but doesn't...
X,Y = m(X,Y)
plt.plot(X,Y,**kwargs)

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

ax = plt.subplot(111)

#Let's create a basemap of the world

m = Basemap(resolution='l',projection='robin',lon_0=0)

m.drawcountries()

m.drawcoastlines()
m.fillcontinents(color='grey',lake_color='white')
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary(fill_color='white')

# Set number 1:
centerlon = 4.360515
centerlat = 50.79747

# Set number 2:
centerlon = -64.360515
centerlat = -30.79747

# Set number 3:
centerlon = 104.360515
centerlat = -40.79747

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

plt.show()```

## 2 thoughts on “Matplotlib Basemap tutorial 09: Drawing circles”

1. Matthew Brown says:

What units must the radii you define be in order for the plotting to be accurate?

2. Thomas Lecocq says:

Matthew,

I think it’s the same units as your map, thus degrees, meters, it depends.