How to use pyshtools - 10 common examples
To help you get started, we’ve selected a few pyshtools examples, based on popular ways it is used in public projects.
Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
marcomusy / vtkplotter / examples / other / spherical_harmonics1.py View on Github 
for th in np.linspace(0, np.pi, N, endpoint=False):
lats = []
for ph in np.linspace(0, 2 * np.pi, N, endpoint=False):
p = spher2cart(rmax, th, ph)
intersections = shape.intersectWithLine([0, 0, 0], p) #
if len(intersections):
value = mag(intersections[0])
lats.append(value)
pts.append(intersections[0])
else:
lats.append(rmax)
pts.append(p)
agrid.append(lats)
agrid = np.array(agrid)
grid = pyshtools.SHGrid.from_array(agrid)
clm = grid.expand()
grid_reco = clm.expand(lmax=lmax) # cut "high frequency" components
#############################################################
# interpolate to a finer grid
agrid_reco = grid_reco.to_array()
# adding 1 column
dim = (lmax * 2 + 2) + 1
gridTmp = np.zeros((dim - 1, dim))
gridTmp[:, 0:dim - 1] = agrid_reco
col0 = gridTmp[:, 0]
gridTmp[:, -1] = col0
agrid_reco = gridTmpdef sphrot_shtools(f, x, lmax=None, latcols=True):
""" Rotate function on sphere f by Euler angles x (Z-Y-Z?) """
if 'pyshtools' not in sys.modules:
raise ImportError('pyshtools not available!')
if latcols:
f = f.T
c = SHGrid.from_array(f).expand()
c_r = c.rotate(*x, degrees=False)
f_r = c_r.expand(lmax=lmax).to_array()
if latcols:
f_r = f_r.T
return f_rshape1 = Sphere(alpha=0.2)
shape2 = vp.load(datadir + "icosahedron.vtk").normalize().lineWidth(1)
agrid1, actorpts1 = makeGrid(shape1, N)
vp.show(shape1, actorpts1, at=0)
agrid2, actorpts2 = makeGrid(shape2, N)
vp.show(shape2, actorpts2, at=1)
vp.camera.Zoom(1.2)
vp.interactive = False
clm1 = pyshtools.SHGrid.from_array(agrid1).expand()
clm2 = pyshtools.SHGrid.from_array(agrid2).expand()
# clm1.plot_spectrum2d() # plot the value of the sph harm. coefficients
# clm2.plot_spectrum2d()
for t in arange(0, 1, 0.005):
act21 = Points(morph(clm2, clm1, t, lmax), c="r", r=4)
act12 = Points(morph(clm1, clm2, t, lmax), c="g", r=4)
vp.show(act21, at=2, resetcam=0)
vp.show(act12, at=3)
vp.camera.Azimuth(2)
vp.show(interactive=1)vp = Plotter(shape=[2, 2], verbose=0, axes=3, interactive=0)
shape1 = Sphere(alpha=0.2)
shape2 = vp.load(datadir + "icosahedron.vtk").normalize().lineWidth(1)
agrid1, actorpts1 = makeGrid(shape1, N)
vp.show(shape1, actorpts1, at=0)
agrid2, actorpts2 = makeGrid(shape2, N)
vp.show(shape2, actorpts2, at=1)
vp.camera.Zoom(1.2)
vp.interactive = False
clm1 = pyshtools.SHGrid.from_array(agrid1).expand()
clm2 = pyshtools.SHGrid.from_array(agrid2).expand()
# clm1.plot_spectrum2d() # plot the value of the sph harm. coefficients
# clm2.plot_spectrum2d()
for t in arange(0, 1, 0.005):
act21 = Points(morph(clm2, clm1, t, lmax), c="r", r=4)
act12 = Points(morph(clm1, clm2, t, lmax), c="g", r=4)
vp.show(act21, at=2, resetcam=0)
vp.show(act12, at=3)
vp.camera.Azimuth(2)
vp.show(interactive=1)we need to use x convention,
alpha_x = -pi (contour clock-wise rotate along z by pi)
beta_x = -pi/2 (contour clock-wise rotate along new x by pi/2)
gamma_x = 0
then y convention is:
alpha_y = alpha_x - pi/2 = 0
beta_y = beta_x = -pi/2
gamma_y = gamma_x + pi/2 = pi/2
reference: https://shtools.oca.eu/shtools/pyshrotaterealcoef.html
'''
new_lighting = np.zeros(lighting.shape)
n = lighting.shape[0]
for i in range(n):
shMatrix = shtools_sh2matrix(lighting[i,:], self.SH_DEGREE)
# rotate coordinate
shMatrix = SHRotateRealCoef(shMatrix, np.array([0, -np.pi/2, np.pi/2]), self.dj)
# rotate object
#shMatrix = SHRotateRealCoef(shMatrix, np.array([-np.pi/2, np.pi/2, -np.pi/2]), self.dj)
new_lighting[i,:] = shtools_matrix2vec(shMatrix)
return new_lighting
we need to use x convention,
alpha_x = pi/2 (clock-wise rotate along z axis by pi/2)
beta_x = -pi/2 (contour clock-wise rotate along new x by pi/2)
gamma_x = 0
then y convention is:
alpha_y = alpha_x - pi/2 = 0
beta_y = beta_x = -pi/2
gamma_y = gamma_x + pi/2 = pi/2
reference: https://shtools.oca.eu/shtools/pyshrotaterealcoef.html
'''
new_lighting = np.zeros(lighting.shape)
n = lighting.shape[0]
for i in range(n):
shMatrix = shtools_sh2matrix(lighting[i,:], self.SH_DEGREE)
# rotate coordinate
shMatrix = SHRotateRealCoef(shMatrix, np.array([0, -np.pi/2, np.pi/2]), self.dj)
# rotate object
#shMatrix = SHRotateRealCoef(shMatrix, np.array([np.pi/2, -np.pi/2, 0]), self.dj)
new_lighting[i,:] = shtools_matrix2vec(shMatrix)
return new_lighting
def __init__(self, coeffs, gm, r0, omega=None, errors=None, lmax=None,
copy=False):
self._coeffs = _SHCoeffs.from_array(coeffs, lmax=lmax, copy=copy)
self.gm = gm
self.r0 = r0
self.omega = omega
self.errors = errors
if self.omega is not None:
self.centrifugal = _Centrifugal(omega=self.omega)def __init__(self, coeffs, gm, r0, errors=None, ell=None, omega=None):
if ell is not None:
self._ell = ell
else:
self._ell = _LevelEllipsoid()
if omega is None:
omega = self._ell.omega
self._coeffs = _SHGravCoeffs.from_array(coeffs=coeffs, gm=gm, r0=r0,
errors=errors, omega=omega, copy=True)def __init__(self):
self.SH_DEGREE = 2
self.dj = djpi2(self.SH_DEGREE)Args:
n (int): dimension is n x n
mode (str): sampling mode; DH or GLQ
Returns:
theta, phi (1D arrays): polar and azimuthal angles
"""
assert n % 2 == 0
j = np.arange(0, n)
if mode == 'DH':
return j*np.pi/n, j*2*np.pi/n
elif mode == 'ours':
return (2*j+1)*np.pi/2/n, j*2*np.pi/n
elif mode == 'GLQ':
from pyshtools.shtools import GLQGridCoord
phi, theta = GLQGridCoord(n-1)
# convert latitude to [0, np.pi/2]
return np.radians(phi+90), np.radians(theta)
elif mode == 'naive':
# repeat first and last points; useful for plotting
return np.linspace(0, np.pi, n), np.linspace(0, 2*np.pi, n)
