<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introduction to pySHTOOLS
pyshtools is a Python module that provides functions for working with spherical harmonics. These routines can be used to perform spherical harmonic transforms and reconstructions, rotations of data expressed in spherical harmonics, and multitaper spectral analyses on the sphere. The base functions are fast routines written in Fortran 95 from the SHTOOLS package. pyshtools provides easy access to these routines by use of Python-wrapper functions and a compact class interface.
To get started, import the standard matplotlib library for graphics, numpy for mathematical extensions to Python, and pyshtools
Step1:  The pyshtools module contains the three classes SHCoeffs, SHGrid, and SHWindow, a submodule shtools that contains all of the Python-wrapped Fortran functions, and a series of submodules that groups these functions (constant,  legendre, expand, shio, spectralanalysis, localizedspectralanalysis, rotate, gravmag, and other). SHCoeffs, SHGrid, and SHWindow bundle the vast majority of pyshtools functionality. In this notebook we will demonstrate how to use the classes related to grids and coefficients.
Grid and coefficient classes
Let's start by first creating a power spectrum that follows a power law with exponent -2, up to and including degree 100. To avoid a division by zero, we will set the degree 0 term to zero
Step2:  Next, we will create a random realization of spherical harmonic coefficients whose expected power is given by the spectrum we just created
Step3:  This creates a new class instance of SHCoeffs that contains several attributes and methods. By default, pyshtools assumes that the coefficients are real and that they are normalized using the '4pi' convention exlcuding the Condon-Shortley phase factor. This is the standard normalization in Geodesy and many fields of geophysics and spectral analysis. Other normalization can be specified explicitly by specifying the optional parameter normalization, which can be either '4pi', 'ortho', or 'schmidt'. The Condon-Shortley phase can be included by setting the optional parameter csphase to -1, and if you wanted complex coefficients, you could set kind='complex'.
from_random() is just one way to create a set of spherical harmonic coefficients. The other constructor methods are from_file() to read the coefficients from a file using the SHRead routine, from_zeros() if you just want all the coefficients to be set to zero, and from_array() if you already have a numpy array of the coefficients.
Next, let's calculate the power spectrum (i.e., the power per spherical harmonic degree) and plot it. pyshtools provides a built in plotting function to do this, and as we will see below, the power spectrum can also be returned as a numpy array using the get_powerperdegree() method.
Step4:  The spherical harmonic coefficients of a given SHCoeffs class instance can be converted easily to a different normalization using the return_coeffs() method. Here, we will convert the coeffients to the orthonormalized convention using the Condon-Shortley phase, which is common in many fields of physics and seismology. Also, let's just return the first 50 degrees of the function
Step5:  If you ever forget how your data are normalized, you can use the built in info() method to remind you
Step6:  We will next calculate the power spectrum of our two functions and plot them along with our theoretical input spectrum
Step7:  As you can see, the power spectrum of the random realization follows closely the input spectrum. Furthermore, the power spectrum is displayed in units of '4pi' spherical harmonics and therefore independend of the employed normalization. You will also note that we used the method get_degrees() to return a numpy list of the degrees from 0 up to the maximum value of the coefficients.
Next, let's expand the data onto a grid and plot it. We first use the expand() method, which returns a new instance of the class SHGrid, and then plot the data using the built in method plot_rawdata()
Step8:  When expanding data onto a grid, pyshtools provides three options. The default is to use grid='DH', which is an equally-sampled grid in latitude and longitude that conforms to Driscoll and Healy's (1994) sampling theorem. If you would like a grid that is oversampled in longitude by a factor or two, such that the number of longitude bands is twice that as the number of latitude bands, use grid='DH2' instead.
The third grid is constructed explicitly for use with Gauss-Legendre quadrature integration techniques. This grid contains about half as many latitude bands as an equivalent DH grid, but the latitudes are unequally sampled at the zeros of the Legendre polynomial of the maximum degree of the coefficients. The following commands show how to expand the spherical harmonic coefficents onto a GLQ grid, plot it, and output lists that contain the latitudes and longitudes (in degrees) for each row and column of the grid
Step9:  Once again, if you ever forget how your grid was constructed, the info() method provides you with everything you need to know
Step10:  Sometimes you need to set individual spherical harmonic coefficients to a specified value. For example, for planetary topography, the degree 2 order 0 term can be atypically large as a result of the planet's rotation. Let's set this term equal to zero and replot it using an over-sampled DH grid. Though we will modify only a single coefficient here, you could also provide input arrays for l, m, and value. When using the set_coeffs() method, you set the cosine terms using positive values for the orders, and the sine terms using negative orders
Step11:  If you want to extract the grid as a numpy array, this can be done using the get_grid() method. If instead you want to extract the coefficients as a numpy array, this can be done using the get_coeffs() method. Here, we will just return the first few degrees of the spherical harmonic coefficients and verify that the coefficient was indeed set to zero
Step12:  It is also easy to rotate either the physical body or the underlying coordinate frame of the function expressed in spherical harmonics. This is accomplished using three Euler angles, alpha, beta and gamma. There are many different conventions for these angles, so please read the documentation before doing this blindly! Here, we will rotate the point at the north pole to 60 degrees north latitude and 180 degrees east longitude
Step13:  By default, this routine expects the angles to be in degrees, which we specified redundantly by the provided optional parameter. If the angles were in radians, you would instead set this to False.
You can also perform simply arithmetic operations on SHCoeffs and SHGrid class instances. Let's create a second set of random coefficients whose power is 10 times less than our first set and expand this onto a grid.
Step14:  We can easily add, subtract multiply, and divide sets of coefficients and grids. If you need to make an exact copy of either, you can use the copy() method.
Step15:  Finally, pyshtools provides many physical constants that are useful when working with planetary data. As an example, here is the GM of Earth, along with information on it's source
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from __future__ import print_function # only necessary if using Python 2.x
import matplotlib.pyplot as plt
import numpy as np
import pyshtools
Explanation: Introduction to pySHTOOLS
pyshtools is a Python module that provides functions for working with spherical harmonics. These routines can be used to perform spherical harmonic transforms and reconstructions, rotations of data expressed in spherical harmonics, and multitaper spectral analyses on the sphere. The base functions are fast routines written in Fortran 95 from the SHTOOLS package. pyshtools provides easy access to these routines by use of Python-wrapper functions and a compact class interface.
To get started, import the standard matplotlib library for graphics, numpy for mathematical extensions to Python, and pyshtools:
End of explanation
degrees = np.arange(101, dtype=float)
degrees[0] = np.inf
power = degrees**(-2)
Explanation: The pyshtools module contains the three classes SHCoeffs, SHGrid, and SHWindow, a submodule shtools that contains all of the Python-wrapped Fortran functions, and a series of submodules that groups these functions (constant,  legendre, expand, shio, spectralanalysis, localizedspectralanalysis, rotate, gravmag, and other). SHCoeffs, SHGrid, and SHWindow bundle the vast majority of pyshtools functionality. In this notebook we will demonstrate how to use the classes related to grids and coefficients.
Grid and coefficient classes
Let's start by first creating a power spectrum that follows a power law with exponent -2, up to and including degree 100. To avoid a division by zero, we will set the degree 0 term to zero:
End of explanation
clm = pyshtools.SHCoeffs.from_random(power)
Explanation: Next, we will create a random realization of spherical harmonic coefficients whose expected power is given by the spectrum we just created:
End of explanation
fig, ax = clm.plot_powerperdegree()
Explanation: This creates a new class instance of SHCoeffs that contains several attributes and methods. By default, pyshtools assumes that the coefficients are real and that they are normalized using the '4pi' convention exlcuding the Condon-Shortley phase factor. This is the standard normalization in Geodesy and many fields of geophysics and spectral analysis. Other normalization can be specified explicitly by specifying the optional parameter normalization, which can be either '4pi', 'ortho', or 'schmidt'. The Condon-Shortley phase can be included by setting the optional parameter csphase to -1, and if you wanted complex coefficients, you could set kind='complex'.
from_random() is just one way to create a set of spherical harmonic coefficients. The other constructor methods are from_file() to read the coefficients from a file using the SHRead routine, from_zeros() if you just want all the coefficients to be set to zero, and from_array() if you already have a numpy array of the coefficients.
Next, let's calculate the power spectrum (i.e., the power per spherical harmonic degree) and plot it. pyshtools provides a built in plotting function to do this, and as we will see below, the power spectrum can also be returned as a numpy array using the get_powerperdegree() method.
End of explanation
clm_ortho = clm.return_coeffs(normalization='ortho', csphase=-1, lmax=50)
Explanation: The spherical harmonic coefficients of a given SHCoeffs class instance can be converted easily to a different normalization using the return_coeffs() method. Here, we will convert the coeffients to the orthonormalized convention using the Condon-Shortley phase, which is common in many fields of physics and seismology. Also, let's just return the first 50 degrees of the function:
End of explanation
clm_ortho.info()
Explanation: If you ever forget how your data are normalized, you can use the built in info() method to remind you:
End of explanation
fix, ax = plt.subplots(1,1)
ax.plot(clm.get_degrees(), power, '-k', 
        clm.get_degrees(), clm.get_powerperdegree(), '-r', 
        clm_ortho.get_degrees(), clm_ortho.get_powerperdegree(), '-b')
ax.set(yscale='log', xlabel='degree', ylabel='power')
Explanation: We will next calculate the power spectrum of our two functions and plot them along with our theoretical input spectrum:
End of explanation
grid = clm.expand()
fig, ax = grid.plot_rawdata()
Explanation: As you can see, the power spectrum of the random realization follows closely the input spectrum. Furthermore, the power spectrum is displayed in units of '4pi' spherical harmonics and therefore independend of the employed normalization. You will also note that we used the method get_degrees() to return a numpy list of the degrees from 0 up to the maximum value of the coefficients.
Next, let's expand the data onto a grid and plot it. We first use the expand() method, which returns a new instance of the class SHGrid, and then plot the data using the built in method plot_rawdata():
End of explanation
grid_glq = clm.expand(grid='GLQ')
grid_glq.plot_rawdata()
lats = grid_glq.get_lats()
lons = grid_glq.get_lons()
Explanation: When expanding data onto a grid, pyshtools provides three options. The default is to use grid='DH', which is an equally-sampled grid in latitude and longitude that conforms to Driscoll and Healy's (1994) sampling theorem. If you would like a grid that is oversampled in longitude by a factor or two, such that the number of longitude bands is twice that as the number of latitude bands, use grid='DH2' instead.
The third grid is constructed explicitly for use with Gauss-Legendre quadrature integration techniques. This grid contains about half as many latitude bands as an equivalent DH grid, but the latitudes are unequally sampled at the zeros of the Legendre polynomial of the maximum degree of the coefficients. The following commands show how to expand the spherical harmonic coefficents onto a GLQ grid, plot it, and output lists that contain the latitudes and longitudes (in degrees) for each row and column of the grid:
End of explanation
grid_glq.info()
grid.info()
Explanation: Once again, if you ever forget how your grid was constructed, the info() method provides you with everything you need to know:
End of explanation
l = 2
m = 0
value = 0.
clm.set_coeffs(value, l, m)
grid_dh2 = clm.expand(grid='DH2')
fig, ax = grid_dh2.plot_rawdata()
grid_dh2.info()
Explanation: Sometimes you need to set individual spherical harmonic coefficients to a specified value. For example, for planetary topography, the degree 2 order 0 term can be atypically large as a result of the planet's rotation. Let's set this term equal to zero and replot it using an over-sampled DH grid. Though we will modify only a single coefficient here, you could also provide input arrays for l, m, and value. When using the set_coeffs() method, you set the cosine terms using positive values for the orders, and the sine terms using negative orders:
End of explanation
coeffs = clm.get_coeffs(lmax=4)
print(coeffs, end='\n\n')
print('c20 = {:f}'.format(coeffs[0,2,0]))
Explanation: If you want to extract the grid as a numpy array, this can be done using the get_grid() method. If instead you want to extract the coefficients as a numpy array, this can be done using the get_coeffs() method. Here, we will just return the first few degrees of the spherical harmonic coefficients and verify that the coefficient was indeed set to zero:
End of explanation
clat = 60.
clon = 180.
alpha = 0.
beta = -(90.-clat)
gamma = -clon
clm_rotated = clm.rotate(alpha, beta, gamma, degrees=True)
grid_rotated = clm_rotated.expand()
fig, ax = grid_rotated.plot_rawdata()
Explanation: It is also easy to rotate either the physical body or the underlying coordinate frame of the function expressed in spherical harmonics. This is accomplished using three Euler angles, alpha, beta and gamma. There are many different conventions for these angles, so please read the documentation before doing this blindly! Here, we will rotate the point at the north pole to 60 degrees north latitude and 180 degrees east longitude:
End of explanation
clm2 = pyshtools.SHCoeffs.from_random(power/10.)
grid2 = clm2.expand()
fig, ax = grid2.plot_rawdata()
Explanation: By default, this routine expects the angles to be in degrees, which we specified redundantly by the provided optional parameter. If the angles were in radians, you would instead set this to False.
You can also perform simply arithmetic operations on SHCoeffs and SHGrid class instances. Let's create a second set of random coefficients whose power is 10 times less than our first set and expand this onto a grid.
End of explanation
grid3 = 100 + 2*grid2 + grid**2 + grid/2 -100 - grid2*2 -grid*grid -grid/2
grid4 = grid3.copy()
Explanation: We can easily add, subtract multiply, and divide sets of coefficients and grids. If you need to make an exact copy of either, you can use the copy() method.
End of explanation
print(pyshtools.constant.gm_earth, end='\n\n')
pyshtools.constant.gm_earth.info()
Explanation: Finally, pyshtools provides many physical constants that are useful when working with planetary data. As an example, here is the GM of Earth, along with information on it's source:
End of explanation
<END_TASK>