Drawing spherical harmonic density plots on the surface of a sphere in tikz/pgfplots

No idea if this could be done "natively" with asymptote, pstricks or TikZ, or even by calculating all the data in an external program and then plotting it with pgfplots. I went for doing everything in python and simply including the resulting image.

This adapts the python code written by Alex J. DeCaria, taken from the ipython notebook linked on this page.

On the python side this requires numpy, scipy, matplotlib and mpl_tookits.basemap, and on the (pdf) latex side it must be compiled with -shell-escape. Most (but not all) options are parametrized with keys that can be called from latex.

\documentclass[tikz,border=5]{standalone}
\usepackage{filecontents}
\begin{filecontents*}{shpl.py}
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
def plot(filename, m, n, longitude=0, latitude=0, inches=(1,1), 
         cmap='RdYlBu', points=100):

    figure, ax = plt.subplots(1,1)
    figure.set_size_inches(*inches)

    lon = np.linspace(0, 2*np.pi, points)
    lat = np.linspace(-np.pi / 2, np.pi / 2, points)
    colat = lat + np.pi / 2
    d = np.zeros((len(lon), len(colat)), dtype=np.complex64)

    meshed_grid = np.meshgrid(lon, lat)
    lat_grid = meshed_grid[1]
    lon_grid = meshed_grid[0]

    mp = Basemap(projection='ortho', lat_0=latitude, lon_0=longitude, ax=ax)
    mp.drawmapboundary()
    mp.drawmeridians(np.arange(0, 360, 30))
    mp.drawparallels(np.arange(-90, 90, 30))

    for j, yy in enumerate(colat):
        for i, xx in enumerate(lon):
            d[i,j] = sp.sph_harm(m, n, xx, yy)

    drm = np.round(np.real(d) / np.max(np.real(d)), 2)
    x, y = mp(np.degrees(lon_grid), np.degrees(lat_grid))
    mp.pcolor(x, y, np.transpose(drm), cmap=cmap)

    figure.savefig(filename, transparent=True)
\end{filecontents*}

\newif\ifshpoverwrite
\tikzset{%
  spherical harmonics/.cd,
    overwrite/.is if=shpoverwrite,
    file/.store in=\shpfilename,
    m/.store in=\shpm,
    n/.store in=\shpn,
    longitude/.store in=\shplongitude,
    latitude/.store in=\shplatitude,
    cmap/.store in=\shpcmap,
    points/.store in=\shppoints,
    inches/.store in=\shpinches,
    longitude=0, latitude=0,
    cmap=RdYlBu,  points=100, inches={(1,1)}
}
\def\sphericalharmonicplot#1{%
  \tikzset{spherical harmonics/.cd,#1}%
  \edef\pythoncommand{python -c "import shpl; 
    shpl.plot('\shpfilename', \shpm, \shpn,
              latitude=\shplatitude, longitude=\shplongitude,
              cmap='\shpcmap', points=\shppoints, inches=\shpinches)"}%
  \ifshpoverwrite
    \immediate\write18{\pythoncommand}%
  \else
    \IfFileExists{\shpfilename}{}{\immediate\write18{\pythoncommand}}%
  \fi%
  \includegraphics{\shpfilename}%
}
\begin{document} 
\begin{tikzpicture}[x=1in,y=1in]
\foreach \m/\n [count=\i from 0] in {0/1, 0/2, 0/3, 1/1, 1/2, 1/3, 
  2/2, 2/3, 3/6, 4/5, 5/7, 6/10}
\node [label=270:{$m=\m,\,n=\n$}] at ({floor(\i/3)*1.5}, {-mod(\i,3)*1.5})
  {\sphericalharmonicplot{file=sph\i.png, m=\m, n=\n,
    longitude=-100, latitude=30}};
\end{tikzpicture}
\end{document}

I could not replicate @Mark solution in Python 3 due to the missing "basemap" in the current Matplotlib module.

I finally found a neat solution which outputs 3 types of plots

Please check the Python code here