In [1]:
''' Generates a sampling of spherical harmonics for illustration of
their properties'''
# Written by Alex DeCaria for submission to the John Hunter Excellence in
# Plotting contest at SciPy 2014.
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import scipy.special as sp
# A collection of indices (order, degree) for the spherical harmonics
indices = ((0,1),(0,2),(0,3),(1,1),(1,2),(1,3),\
(2,2),(2,3),(3,6),(4,5),(5,7),(6,10))
#Setup figure and plot parameters
pcolor = True # If True, plots pcolor image instead of filled contour
cmap = cm.get_cmap('BrBG_r') # Set color map
figure, ax = plt.subplots(4,3)
figure.set_size_inches(6, 8)
figure.set_facecolor('#ffde99')
# Create latitude and longitude arrays
lon = np.linspace(0,2*np.pi,100)
lat = np.linspace(-np.pi/2,np.pi/2,100)
colat = lat+np.pi/2 # colatitude array
# create a 2-D array to hold spherical harmonic values
d = np.zeros((len(lon),len(colat)),dtype = np.complex64)
# Generate 2-D lat/lon grids
meshed_grid = np.meshgrid(lon, lat)
lat_grid = meshed_grid[1]
lon_grid = meshed_grid[0]
for i, a in enumerate(ax.flat):
# set up orthographic map projection
mp = Basemap(projection='ortho', lat_0 = 30, lon_0 = -100,ax=a)
# draw the edge of the map projection region (the projection limb)
mp.drawmapboundary()
# draw lat/lon grid lines every 30 degrees.
mp.drawmeridians(np.arange(0, 360, 30))
mp.drawparallels(np.arange(-90, 90, 30))
# Generate spherical harmonic data
m, n = indices[i]
print('Creating plot for m = {0:d}, n = {1:d}'.format(m,n))
for j, yy in enumerate(colat):
for i, xx in enumerate(lon):
d[i,j] = sp.sph_harm(m,n,xx,yy)
# Create normalized data for real component
drm = np.round(np.real(d)/np.max(np.real(d)), 2)
# Create projection coordinates
x, y = mp(np.degrees(lon_grid), np.degrees(lat_grid))
#Plot data
if pcolor:
CS = mp.pcolor(x,y,np.transpose(drm),cmap=cmap)
else:
CS = mp.contourf(x,y,np.transpose(drm),20,cmap=cmap)
a.annotate('{0:d},{1:d}'.format(m,n),(-0.05,0.95),xycoords = 'axes fraction')
# Create and plot colorbar
cax = figure.add_axes([0.15,0.03,0.7,0.03])
cb = plt.colorbar(CS, cax=cax,ticks = np.linspace(-1.0, 1.0, 5),
orientation = 'horizontal')
cax.set_xlabel('Normalized Amplitude',size = 'x-large', labelpad = -55)
#figure.savefig('Spherical-harmonics.pdf', facecolor = figure.get_facecolor(), edgecolor='none')
plt.show()
In []: