.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/spherical_harmonics_basis.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_spherical_harmonics_basis.py: Understanding spherical harmonics ================================= Spherical harmonics are the spherical counterpart of Fourier modes. Just as Fourier modes provide a basis for representing periodic functions, spherical harmonics provide a basis for representing functions defined on the sphere. This tutorial introduces their main mathematical properties, visualises some basis functions, and shows how ``S2FFT`` decomposes a signal into spherical harmonic coefficients. .. image:: https://colab.research.google.com/assets/colab-badge.svg :align: center :alt: Open in Google Colab :target: https://colab.research.google.com/github/astro-informatics/s2fft/tree/gh-pages/_colab_notebooks/spherical_harmonics_basis.ipynb If you are working on this notebook in Google Colab, you will need to have Google Colab install ``s2fft``. You can do this by adding a cell to the top of the notebook with the following content: .. code-block:: bash !pip install cartopy s2fft &> /dev/null and then running that cell. .. GENERATED FROM PYTHON SOURCE LINES 32-46 Spherical coordinates --------------------- A point on the unit sphere is described by two angles: - :math:`\theta \in [0,\pi]` is the colatitude, measured down from the north pole. - :math:`\phi \in [0,2\pi)` is the longitude around the sphere, measured eastward from a reference `meridian `__. We write a signal on the sphere as :math:`f(\theta,\phi)`. Depending on the application, it could represent quantities such as temperature, elevation, or radiation intensity. .. GENERATED FROM PYTHON SOURCE LINES 48-72 Spherical harmonics ------------------- The spherical harmonic :math:`Y_{\ell m}(\theta,\phi)` is indexed by two integers: - The degree :math:`\ell=0,1,2,\dots` controls its overall angular scale. - The order :math:`m=-\ell,\dots,\ell` controls its variation with longitude. A spherical harmonic can be written as .. math:: Y_{\ell m}(\theta,\phi) = N_{\ell m} P_{\ell}^m(\cos\theta) \exp(i m\phi), where :math:`N_{\ell m}` is a normalisation constant, :math:`P_{\ell}^{m}` is an `associated Legendre polynomial `__ controlling variation with colatitude, and :math:`\exp(i m\phi)` controls variation with longitude. Consequently, spherical harmonics are generally complex-valued. .. GENERATED FROM PYTHON SOURCE LINES 74-114 Orthonormality -------------- The surface-area element on the unit sphere, also called the `solid-angle element `__, is .. math:: \mathrm{d}\Omega = \sin\theta \ \mathrm{d}\theta \mathrm{d}\phi. The factor :math:`\sin\theta` comes from expressing surface area in spherical coordinates. For two square-integrable functions :math:`f` and :math:`g`, the inner product on the sphere is .. math:: \langle f,g\rangle = \int_{\mathbb{S}^{2}} f(\theta,\phi) g^{*}(\theta,\phi) \ \mathrm{d}\Omega. Here, :math:`*` denotes complex conjugation. The spherical harmonics are orthonormal with respect to the inner product: .. math:: \langle Y_{\ell m},Y_{\ell' m'}\rangle = \int_{\mathbb{S}^{2}} Y_{\ell m}(\theta,\phi) Y^{*}_{\ell' m'}(\theta,\phi) \ \mathrm{d}\Omega = \delta_{\ell \ell'} \delta_{m m'}. The Kronecker delta :math:`\delta_{ab}` equals one when :math:`a=b` and zero otherwise. Consequently, distinct spherical harmonics have zero inner product, while the inner product of a spherical harmonic with itself is one. .. GENERATED FROM PYTHON SOURCE LINES 116-155 Expanding a signal in spherical harmonics ----------------------------------------- The spherical harmonics form a complete orthonormal basis for :math:`L^{2}(\mathbb{S}^{2})`, the space of `square-integrable functions `__ on the sphere. Any signal in this space can therefore be expanded as .. math:: f(\theta,\phi) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} f_{\ell m} Y_{\ell m}(\theta,\phi), where orthonormality allows each coefficient to be isolated by an inner product: .. math:: f_{\ell m} = \langle f,Y_{\ell m}\rangle = \int_{\mathbb{S}^{2}} f(\theta,\phi) Y^{*}_{\ell m}(\theta,\phi) \ \mathrm{d}\Omega. Computing the coefficients :math:`f_{\ell m}` is called the forward spherical harmonic transform, or spherical harmonic analysis. Reconstructing :math:`f` from these coefficients is called the inverse spherical harmonic transform, or spherical harmonic synthesis. For a real-valued signal, the coefficients satisfy the conjugate-symmetry relation .. math:: f_{\ell,-m} = (-1)^{m} f_{\ell m}^{*}. An individual spherical harmonic is generally complex-valued, so we set ``reality=False``. .. GENERATED FROM PYTHON SOURCE LINES 157-179 Band-limited signals -------------------- ``S2FFT`` computes transforms for signals with a finite harmonic band-limit :math:`L`. A signal is band-limited at :math:`L` if :math:`f_{\ell m}=0` for every :math:`\ell \geq L`. Its expansion then becomes the finite sum .. math:: f(\theta,\phi) = \sum_{\ell=0}^{L-1} \sum_{m=-\ell}^{\ell} f_{\ell m} Y_{\ell m}(\theta,\phi). The largest degree represented is therefore :math:`L-1`. Increasing :math:`L` allows finer angular structure to be represented, but also increases the number of coefficients and the computational cost of the transforms. A sampling theorem specifies a finite set of sample locations at which a band-limited signal can be represented and transformed exactly, up to numerical precision. Here we use the `McEwen-Wiaux (2012) `__ sampling scheme. .. GENERATED FROM PYTHON SOURCE LINES 181-195 .. code-block:: Python import jax jax.config.update("jax_enable_x64", True) import cartopy.crs as ccrs import numpy as np from matplotlib import pyplot as plt import s2fft L = 32 sampling = "mw" .. GENERATED FROM PYTHON SOURCE LINES 196-214 How S2FFT stores harmonic coefficients -------------------------------------- S2FFT stores coefficients in an array ``flm`` of shape ``(L, 2 * L - 1)``. The first array index corresponds directly to the degree :math:`\ell`. The order :math:`m` is shifted by :math:`L-1` so that its negative and positive values can be represented by standard non-negative array indices: .. math:: f_{\ell m} \quad \longleftrightarrow \quad \mathtt{flm[\ell,\ m+L-1]}. Thus, :math:`m=0` is stored in the central column with index :math:`L-1`. The rectangular array also contains entries for which :math:`|m|>\ell`. These do not correspond to valid spherical harmonic coefficients and are kept at zero. .. GENERATED FROM PYTHON SOURCE LINES 216-217 Choose a valid degree and order: 0 <= ell < L, and -ell <= m <= ell. .. GENERATED FROM PYTHON SOURCE LINES 217-228 .. code-block:: Python ell = 2 m = 1 # Use S2FFT's helper function to get the coefficient array shape described above. flm_shape = s2fft.sampling.s2_samples.flm_shape(L) # Set all coefficients to zero, then set f_{ell m} = 1. flm = np.zeros(flm_shape, dtype=np.complex128) flm[ell, m + L - 1] = 1.0 .. GENERATED FROM PYTHON SOURCE LINES 229-242 Since :math:`f_{\ell m}=1` is the only non-zero coefficient, the spherical harmonic expansion reduces to .. math:: f(\theta,\phi) = 1 \cdot Y_{\ell m}(\theta,\phi) = Y_{\ell m}(\theta,\phi). The inverse transform therefore evaluates the chosen basis function at the MW sample locations. An individual spherical harmonic is generally complex-valued, so we set ``reality=False``. .. GENERATED FROM PYTHON SOURCE LINES 244-253 .. code-block:: Python y_ell_m = s2fft.inverse( flm, L=L, sampling=sampling, method="jax", reality=False, ) .. GENERATED FROM PYTHON SOURCE LINES 254-256 Spherical harmonics are generally complex-valued, so we visualise the real part. .. GENERATED FROM PYTHON SOURCE LINES 258-274 .. code-block:: Python fig, ax = plt.subplots( figsize=(3, 1.5), subplot_kw={"projection": ccrs.Mollweide()}, ) ax.imshow( y_ell_m.real, transform=ccrs.PlateCarree(), cmap="viridis", ) ax.set_title(rf"$\ell={ell},\ m={m}$") fig.tight_layout() plt.show() .. image-sg:: /tutorials/images/sphx_glr_spherical_harmonics_basis_001.png :alt: $\ell=2,\ m=1$ :srcset: /tutorials/images/sphx_glr_spherical_harmonics_basis_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 275-281 Visualising the basis --------------------- We can repeat the same construction, setting one coefficient to one and all remaining coefficients to zero, for every valid degree and order :math:`(\ell,m)` up to a chosen maximum degree. .. GENERATED FROM PYTHON SOURCE LINES 283-325 .. code-block:: Python max_degree = 3 ell_values = np.arange(max_degree + 1) m_values = np.arange(-max_degree, max_degree + 1) fig, axes = plt.subplots( len(ell_values), len(m_values), figsize=(16, 8), subplot_kw={"projection": ccrs.Mollweide()}, ) for ell in ell_values: for m in m_values: ax = axes[ell, m + max_degree] # Hide positions that do not correspond to a valid order. if abs(m) > ell: ax.set_axis_off() continue flm = np.zeros(flm_shape, dtype=np.complex128) flm[ell, m + L - 1] = 1.0 y_ell_m = s2fft.inverse( flm, L=L, sampling=sampling, method="jax", reality=False, ) ax.imshow( y_ell_m.real, transform=ccrs.PlateCarree(), cmap="viridis", ) ax.set_title(rf"$\ell={ell},\ m={m}$", fontsize=9) fig.tight_layout() plt.show() .. image-sg:: /tutorials/images/sphx_glr_spherical_harmonics_basis_002.png :alt: $\ell=0,\ m=0$, $\ell=1,\ m=-1$, $\ell=1,\ m=0$, $\ell=1,\ m=1$, $\ell=2,\ m=-2$, $\ell=2,\ m=-1$, $\ell=2,\ m=0$, $\ell=2,\ m=1$, $\ell=2,\ m=2$, $\ell=3,\ m=-3$, $\ell=3,\ m=-2$, $\ell=3,\ m=-1$, $\ell=3,\ m=0$, $\ell=3,\ m=1$, $\ell=3,\ m=2$, $\ell=3,\ m=3$ :srcset: /tutorials/images/sphx_glr_spherical_harmonics_basis_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 326-330 Moving down the rows increases :math:`\ell`, producing finer angular structure. Moving across a row changes :math:`m`. Increasing :math:`|m|` produces more variation with longitude. When :math:`m=0`, the harmonic does not vary with longitude. .. GENERATED FROM PYTHON SOURCE LINES 332-346 Decomposing a signal into spherical harmonics --------------------------------------------- A band-limited signal is a weighted sum of spherical harmonic basis functions. To demonstrate this, we construct a signal with three non-zero coefficients: .. math:: f(\theta, \phi) = Y_{0,0}(\theta,\phi) + 0.8Y_{2,1}(\theta,\phi) - 0.4Y_{3,-2}(\theta,\phi). The inverse transform evaluates this weighted sum at the MW sampling nodes. .. GENERATED FROM PYTHON SOURCE LINES 348-362 .. code-block:: Python signal_flm = np.zeros(flm_shape, dtype=np.complex128) signal_flm[0, 0 + L - 1] = 1.0 signal_flm[2, 1 + L - 1] = 0.8 signal_flm[3, -2 + L - 1] = -0.4 signal = s2fft.inverse( signal_flm, L=L, sampling=sampling, method="jax", reality=False, ) .. GENERATED FROM PYTHON SOURCE LINES 363-365 These coefficients do not satisfy conjugate symmetry, so the resulting signal is complex-valued. We visualise its real part below. .. GENERATED FROM PYTHON SOURCE LINES 367-383 .. code-block:: Python fig, ax = plt.subplots( figsize=(3, 1.5), subplot_kw={"projection": ccrs.Mollweide()}, ) ax.imshow( signal.real, transform=ccrs.PlateCarree(), cmap="viridis", ) ax.set_title(r"$f(\theta,\phi)$") fig.tight_layout() plt.show() .. image-sg:: /tutorials/images/sphx_glr_spherical_harmonics_basis_003.png :alt: $f(\theta,\phi)$ :srcset: /tutorials/images/sphx_glr_spherical_harmonics_basis_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 384-390 Recovering the spherical harmonic coefficients ---------------------------------------------- Starting from the sampled signal, the forward transform recovers its spherical harmonic coefficients, which are the weights of the basis functions. .. GENERATED FROM PYTHON SOURCE LINES 392-401 .. code-block:: Python recovered_flm = s2fft.forward( signal, L=L, sampling=sampling, method="jax", reality=False, ) .. GENERATED FROM PYTHON SOURCE LINES 402-405 We display the magnitudes of the recovered coefficients up to degree three. The non-zero entries should occur at the same degree and order pairs used to construct the signal. .. GENERATED FROM PYTHON SOURCE LINES 407-443 .. code-block:: Python coefficient_magnitudes = np.abs( recovered_flm[: max_degree + 1, L - 1 - max_degree : L + max_degree] ) # Leave invalid coefficients blank coefficient_magnitudes[abs(m_values[None]) > ell_values[:, None]] = np.nan fig, ax = plt.subplots(figsize=(7, 4)) im = ax.pcolormesh( m_values, ell_values, coefficient_magnitudes, cmap="viridis", edgecolors="white", linewidth=2, ) # Label the non-zero coefficients non_zero = ~np.isclose(np.nan_to_num(coefficient_magnitudes), 0) for ell, m_index in np.argwhere(non_zero): ax.text( m_values[m_index], ell, f"{coefficient_magnitudes[ell, m_index]:.2g}", ha="center", va="center", ) ax.invert_yaxis() ax.set(aspect=1, xlabel=r"$m$", ylabel=r"$\ell$", xticks=m_values, yticks=ell_values) fig.colorbar(im, ax=ax, label=r"$|f_{\ell m}|$", shrink=0.8) fig.tight_layout() plt.show() .. image-sg:: /tutorials/images/sphx_glr_spherical_harmonics_basis_004.png :alt: spherical harmonics basis :srcset: /tutorials/images/sphx_glr_spherical_harmonics_basis_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 444-447 Finally, we compare the recovered coefficients with the original coefficients. For the MW sampling scheme, the error should be close to machine precision. .. GENERATED FROM PYTHON SOURCE LINES 449-451 .. code-block:: Python max_error = np.max(np.abs(recovered_flm - signal_flm)) print(f"Maximum coefficient error: {max_error:.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Maximum coefficient error: 7.08e-15 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 17.032 seconds) .. _sphx_glr_download_tutorials_spherical_harmonics_basis.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: spherical_harmonics_basis.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: spherical_harmonics_basis.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: spherical_harmonics_basis.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_