PyHarm#

This chapter demonstrates how to perform some basic tasks with PyHarm. For a detailed description of the classes, methods, functions, etc. to be used, see the PyHarm (Python API) chapter. The source codes of the examples can be found in the cookbook/python directory.

In your Python code, import PyHarm in double precision by:

>>> import pyharm as ph

To import PyHarm in single precision, use:

>>> import pyharmf as phf

Both precisions can be combined in a single Python file, too:

>>> import pyharm as ph
>>> import pyharmf as phf

The chapters that follow show how to work with PyHarm in double and single precision, respectively. Quadruple precision is not supported in PyHarm.

PyHarm in double precision#

PyHarm attempts to follow the rules of the object-oriented programming paradigm. Several classes are introduced to represent spherical harmonic coefficients, coordinates of evaluation points/cells and Fourier coefficients of Legendre functions. Most of the classes are accompanied by factory methods for instantiation, so users should always instantiate the classes using these methods. The factory methods are summarized within the class documentation. Alternatively, they can be found by searching for the classmethod flag that precedes the class method names in the documentation.

Integration with NumPy#

PyHarm supports NumPy arrays. All arrays consisting of floating point data must

be instances of the numpy.ndarray class,
have the dtype flag set to numpy.float64 (or to numpy.float32 if the library was compiled in single precision), and
have the flags.c_contiguous attribute set to True.

As an example, let’s create a 1D numpy array x that meets all of these conditions.

>>> import numpy as np
>>> x = np.zeros((5,))
>>> x.dtype
dtype('float64')
>>> x.flags.c_contiguous
True

In this code snippet, we did not specify the dtype and order flags when calling numpy.zeros, so numpy automatically used their default values, namely numpy.float64 and 'C'. In some cases, you might need to explicitly specify the two flags:

>>> import numpy as np
>>> x = np.zeros((5,), dtype=np.float64, order='C')
>>> x.dtype
dtype('float64')
>>> x.flags.c_contiguous
True

For further details, see the NumPy documentation.

If an array passed to PyHarm does not meet all the required criteria, PyHarm will throw an error.

Working with spherical harmonic coefficients#

This example shows how to read/write spherical harmonic coefficients from/to files and how to compute (difference) degree variances.

import pyharm as ph


# INPUTS
# =============================================================================
# Define the path to an input "gfc" text file with spherical harmonic
# coefficients.  For details on the structure of the "gfc" file, see the
# documentation.
shcs_in_file = '../../data/input/EGM96-degree10.gfc'


# Maximum harmonic degree to read the spherical harmonic coefficients
nmax = 10;


# Define the path to an output binary file with spherical harmonic
# coefficients.  For details on the structure of the output binary file, see
# the documentation of the "shc" module of "CHarm".
shcs_out_file = '../../data/output/EGM96-degree10.shcs'
# =============================================================================


# =============================================================================
# Read the spherical harmonic coefficients using the factory method "from_file"
# of the "Shc" class from the "ph.shc" module.
# -----------------------------------------------------------------------------
shcs = ph.shc.Shc.from_file('gfc', shcs_in_file, nmax)
# -----------------------------------------------------------------------------


# Now print some more or less randomly chosen spherical harmonic coefficients
# -----------------------------------------------------------------------------
# Let's start with zonal coefficients
n = 9
m = 0
cnm, snm = shcs.get_coeffs(n, m)
print(f'C({n}, {m}) = {cnm}')
print(f'S({n}, {m}) = {snm}')


# Now some tesseral coefficients
n = 9
m = 4
cnm, snm = shcs.get_coeffs(n, m)
print(f'C({n}, {m}) = {cnm}')
print(f'S({n}, {m}) = {snm}')


# And finally some sectorial coefficients
n = 9
m = 9
cnm, snm = shcs.get_coeffs(n, m)
print(f'C({n}, {m}) = {cnm}')
print(f'S({n}, {m}) = {snm}')
# -----------------------------------------------------------------------------


# Now let's save the coefficients to a binary file and then read them back to
# a new "Shc" class instance.  Just for the fun...
# -----------------------------------------------------------------------------
shcs.to_file('bin', nmax, shcs_out_file)


# Read back the coefficients to a new class instance called "shcs2"
shcs2 = ph.shc.Shc.from_file('bin', shcs_out_file, nmax)
# -----------------------------------------------------------------------------


# Compute degree variances from "shcs"
# -----------------------------------------------------------------------------
dv = shcs.dv()


n = 0
print('\n')
print(f'Degree variance for harmonic degree {n} = {dv[n]}')
n = 4
print(f'Degree variance for harmonic degree {n} = {dv[n]}')
n = 10
print(f'Degree variance for harmonic degree {n} = {dv[n]}')
# -----------------------------------------------------------------------------


# Now check whether "shcs" and "shcs2" contain the same coefficients by
# computing difference degree variances
# -----------------------------------------------------------------------------
ddv = shcs.ddv(shcs2)


n = 0
print('\n')
print(f'Difference degree variance for harmonic degree {n} = {ddv[n]}')
n = 4
print(f'Difference degree variance for harmonic degree {n} = {ddv[n]}')
n = 10
print(f'Difference degree variance for harmonic degree {n} = {ddv[n]}')
# -----------------------------------------------------------------------------


print('Great, all done!')
# =============================================================================

Spherical harmonic synthesis and analysis#

This section shows several examples.

A closed-loop test of analysis and synthesis with point data values. First, spherical harmonic coefficients are loaded from a text file. Then, they are used to synthesize a signal, from which a new set of coefficients is finally computed by surface spherical harmonic analysis. The two coefficients sets are then finally compared (should be equal, that is, the difference degree amplitudes should be negligibly small).
A solid spherical harmonic synthesis of point values at scattered points.
A solid spherical harmonic synthesis of point values at a custom grid of points.
A solid spherical harmonic synthesis of block-mean values at scattered cells.
A solid spherical harmonic synthesis of block mean values at a grid of cells.
Surface spherical harmonic analysis with block-mean values in cells.

import numpy as np
import pyharm as ph


# INPUTS
# =============================================================================
# Define the path to an input text file with spherical harmonic coefficients.
# For details on the structure of the text file, see the documentation.
shcs_in_file = '../../data/input/EGM96-degree10-mtx.txt'


# Maximum harmonic degree of coefficients to read.  The same degree is used
# later in the harmonic synthesis and harmonic analysis.
nmax = 10
# =============================================================================


# =============================================================================
print('Closed-loop experiment -- Spherical harmonic synthesis and analysis')
print('===========================')


# Read spherical harmonic coefficients from the input text file.
shcs = ph.shc.Shc.from_file('mtx', shcs_in_file, nmax)


# Now let's say we do not want to use the zero-degree term.  This can be
# achieved simply by setting the respective "C00" coefficient to zero.
shcs.set_coeffs(n=0, m=0, c=0.0)


# Compute the Gauss--Legendre grid for a given "nmax" on a sphere that is
# "1000.0" metres above the reference sphere of spherical harmonic
# coefficients.  PyHarm represents the Gauss--Legendre grid via the
# "PointGridGL" class.  For the Driscoll--Healy grids, PyHarm has the
# "PointGridDH1" and "PointGridDH2" classes (see the documentation).
grd_pnt = ph.crd.PointGridGL(nmax, shcs.r + 1000.0)


# Perform the synthesis.  Since the Gauss--Legendre grid resides "1000.0"
# metres above the reference sphere of the coefficients, performed is
# solid spherical harmonic synthesis.
f = ph.shs.point(grd_pnt, shcs, nmax)


# Note on parallelization
# -----------------------
#
# If you compiled the library with OpenMP parallelization enabled, you can
# control the number of threads being used either by an environmental variable
# or in your Python code, for instance, like this:
#
# from threadpoolctl import threadpool_limits
# with threadpool_limits(limits=2):
#     f = ph.shs.point(grd_pnt, shcs, nmax)


# Print some synthesized values
print('Print some synthesized values of the signal...')
i = 0
j = 0
print(f'f({i}, {j}) = {f[i, j]}')


i = grd_pnt.lat.size // 2
j = grd_pnt.lat.size // 2
print(f'f({i}, {j}) = {f[i, j]}')
print('')


# Now use the synthesized signal and compute back its spherical harmonic
# coefficients by harmonic analysis.  The output coefficients in "shcs2" should
# be the same as the input ones in "shcs".  Note that the signal "f" resides on
# a sphere that is "1000.0" metres above the reference sphere "shcs.r".  The
# last input parameter therefore allows user to define the reference radius of
# the output coefficients.  In this case, this value is set to "shcs.r" to
# ensure that "shcs2.r" will be the same as "shcs.r".
shcs2 = ph.sha.point(grd_pnt, f, nmax, shcs.mu, shcs.r)


# Now check whether "shcs" and "shcs2" are the same by computing their
# difference degree amplitudes
dda = shcs.dda(shcs2)


# Print some difference degree amplitudes
n = 0
print('Now print the difference degree amplitudes.  These should be'
      'very small, say, at the order of 1e-18 or less...')
print(f'Difference degree amplitude for harmonic degree {n} = {dda[n]}')
n = 4
print(f'Difference degree amplitude for harmonic degree {n} = {dda[n]}')
n = 10
print(f'Difference degree amplitude for harmonic degree {n} = {dda[n]}')


# We will not need some of the variables anymore, so let's free the memory by
# deleting them
del grd_pnt, shcs2, f, dda


print('===========================')
# =============================================================================


# =============================================================================
print('')
print('Solid spherical harmonic synthesis at scattered points')
print('===========================')


# Create a "PointSctr" class instance to store three scattered points.  The
# instance will be created from three numpy arrays that hold the coordinates of
# the scattered points.
lat      = np.array([0.1, 0.436231, -0.9651], dtype=np.float64)
lon      = np.array([0.0, 1.53434, 4.2316], dtype=np.float64)
r        = np.array([1000.0, 2000.0, 3000.0], dtype=np.float64) + shcs.r
sctr_pnt = ph.crd.PointSctr.from_arrays(lat, lon, r)


# Do the synthesis
f = ph.shs.point(sctr_pnt, shcs, nmax)


# It's really this easy!


# Print the synthesized values
print(f'{f}')


del sctr_pnt, f
print('===========================')
# =============================================================================


# =============================================================================
print('')
print('Solid spherical harmonic synthesis at a custom grid of points')
print('===========================')

# Define some grid points as numpy arrays.  Then create a "PointGrid" class
# instance that is designed to store custom point grids.
nlat    = 5
nlon    = 10
lat     = np.linspace(np.pi / 2.0, -np.pi / 2.0, nlat)
lon     = np.linspace(0.0, 2.0 * np.pi, nlon)
r       = np.zeros((nlat,), dtype=np.float64) + shcs.r + np.arange(nlat)
grd_pnt = ph.crd.PointGrid.from_arrays(lat, lon, r)


# Do the synthesis
f = ph.shs.point(grd_pnt, shcs, nmax)


# Print the synthesized values
print(f'{f}')


del grd_pnt, f
print('===========================')
# =============================================================================


# =============================================================================
print('')
print('Solid spherical harmonic synthesis at scattered cells')
print('===========================')


# Now we define evaluation cells to synthesize area mean values.  To this end,
# we specify numpy arrays of minimum and maximum latitudes and longitudes and
# an array of spherical radii.  Using these arrays, we create a "CellSctr"
# class instance.
latmax    = np.array([0.323413435, -0.90234320952, 0.0])
latmin    = latmax - np.array([0.234323, 0.4456, np.pi / 2.0])
lonmin    = np.array([0.123456789, 4.3445234, 0.0])
lonmax    = lonmin + np.array([1.3235, 0.01, 2.0 * np.pi])
r         = np.zeros((latmin.size,), dtype=np.float64) + shcs.r + \
            np.arange(latmax.size) + 1000.0
cell_sctr = ph.crd.CellSctr.from_arrays(latmin, latmax, lonmin, lonmax, r)


# Do the synthesis
f = ph.shs.cell(cell_sctr, shcs, nmax)


# Print the synthesized values
print(f'{f}')


del cell_sctr, f
print('===========================')
# =============================================================================


# =============================================================================
print('')
print('Solid spherical harmonic synthesis at a grid of cells')
print('===========================')


# Define the grid cells and create a "CellGrid" class instance from numpy
# arrays.
nlat     = 15
nlon     = 30
latmax   = np.linspace(np.pi / 2.0, -np.pi / 2.0, nlat, False)
latmin   = np.linspace(latmax[1], -np.pi / 2.0, nlat)
lonmin   = np.linspace(0.0, 2.0 * np.pi, nlon, False)
lonmax   = np.linspace(lonmin[1], 2.0 * np.pi, nlon)
r        = np.zeros((nlat,), dtype=np.float64) + shcs.r
grd_cell = ph.crd.CellGrid.from_arrays(latmin, latmax, lonmin, lonmax, r)


# Do the synthesis
f = ph.shs.cell(grd_cell, shcs, nmax)


# Print some of the synthesized values
# -----------------------------------------------------------------------------
i = 0
j = 10
print(f'f({i}, {j}) = {f[i, j]}')


i = 4
j = 3
print(f'f({i}, {j}) = {f[i, j]}')
# -----------------------------------------------------------------------------

print('===========================')
# =============================================================================


# =============================================================================
print('')
print('Spherical harmonic analysis of block-mean values in cells')
print('===========================')


# In this example, we use the "grd_cell" instance and the signal "f" from the
# previous example


# Do the analysis using the approximate quadrature
shcs2 = ph.sha.cell(grd_cell, f, nmax, ph.sha.CELL_AQ, shcs.mu, shcs.r)


# Print some of the computed coefficients.  Note that the harmonic analysis
# with block-mean values in cells is *not* exact, hence the coefficients will
# not be equal to the input ones.
n = [2, 9, 9, 9]  # List of harmonic degrees
m = [0, 0, 4, 9]  # List of harmonic orders
cnm, snm = shcs2.get_coeffs(n, m)
for i, (n1, m1) in enumerate(zip(n, m)):
    print(f'C({n1}, {m1}) = {cnm[i]}')
    print(f'S({n1}, {m1}) = {snm[i]}')


del grd_cell, shcs2, f, shcs


print('===========================')
# =============================================================================

First- and second-order gradients (gravitational vector and tensor)#

This example shows how to compute the full gravitational vector and the full gravitational tensor in the Local north-oriented reference frame.

import pyharm as ph
import numpy as np


# INPUTS
# =============================================================================
# Path to the input file with spherical harmonic coefficients in the 'gfc'
# format
shcs_file = '../../data/input/EGM96-degree10.gfc'
# =============================================================================


# =============================================================================
# First, let's get the maximum harmonic degree stored in "shcs_file".  This is
# especially useful to read all coefficients in "shcs_file" without knowing its
# maximum harmonic degree a priori.
nmax = ph.shc.Shc.nmax_from_file('gfc', shcs_file)


# Now read all coefficients in "shcs_file"
shcs = ph.shc.Shc.from_file('gfc', shcs_file, nmax)


# Compute the Gauss--Legendre point grid for an "nmax" given by the maximum
# degree stored in "shcs" and for a radius equal to the reference sphere,
# to which the spherical harmonic coefficients are scaled to.  We intentionally
# use here the Gauss--Legendre grid instead of the Driscoll--Healy grids in
# order to avoid the inaccuracies due to the singularities at the poles (see
# the documentation).
grd = ph.crd.PointGridGL(shcs.nmax, shcs.r)


# Compute the full first-order gradient (vector)
fx, fy, fz = ph.shs.point_grad1(grd, shcs, shcs.nmax)


# Compute the magnitude of the gravitational acceleration (no contribution due
# to the centrifugal force is considered here)
g = np.sqrt(fx**2 + fy**2 + fz**2)


# Now compute the full second-order gradient (tensor)
fxx, fxy, fxz, fyy, fyz, fzz = ph.shs.point_grad2(grd, shcs, shcs.nmax)


# Check whether "fxx + fyy + fzz" is zero within numerical errors
print(f'The largest error of the gravitational tensor trace is '
      f'{np.abs(fxx + fyy + fzz).max()} s**-2')


print('Great, all done!')
# =============================================================================

Fourier coefficients of Legendre functions#

import pyharm as ph


# Maximum harmonic degree to compute the Fourier coefficients of Legendre
# functions
nmax = 500


# Create a "ph.leg.Pnmj" class instance using the factory method called
# "from_zeros"
pnmj = ph.leg.Pnmj.from_zeros(nmax, ph.leg.PMNJ)


# Compute the Fourier coefficients
pnmj.coeffs()


# Print some Fourier coefficients
# Harmonic degree
n = 123
# Harmonic order
m = 23
# Wave-number-related variable
j = 12
k = pnmj.j2k(n, j)
c = pnmj.get_coeff(n, m, j)
print(f'Fourier coefficients for degree {n}, order {m} and wave-number {k} = '
      f'{c}')


n = 360
m = 358
j = 101
k = pnmj.j2k(j, k)
c = pnmj.get_coeff(n, m, j)
print(f'Fourier coefficients for degree {n}, order {m} and wave-number {k} = '
      f'{c}')


print('Great, all done!')

Integrals#

Finally, let’s compute an integral of a product of two spherical harmonics over a restricted domain on the unit sphere. Then, we do the same, but with the product of Legendre functions only.

import pyharm as ph


# Type of the first spherical harmonic function, its harmonic degree and order,
# respectively
i1 = 0
n1 = 287
m1 = 122


# Type of the second spherical harmonic function, its harmonic degree and
# order, respectively
i2 = 1
n2 = 34
m2 = 9


# Minimum and maximum co-latitudes of the integration domain
cltmin = 0.1
cltmax = 0.9


# Minimum and maximum longitudes of the intration domain
lonmin = 0.5
lonmax = 0.6


# Compute the Fourier coefficients
nmax = max(n1, n2)
pnmj = ph.leg.fourier_coeffs(nmax)


# Compute the integral of a product of two spherical harmonics
iy = ph.integ.yi1n1m1yi2n2m2(cltmin, cltmax, lonmin, lonmax,
                             i1, n1, m1, i2, n2, m2, pnmj)


# Print the value of the integral
print(f'The integral of the product of two spherical harmonics '
      f'i1 = {i1}, n1 = {n1}, m1 = {m1}, i2 = {i2}, n2 = {n2}, m2 = {m2} '
      f'is {iy}')


# Now compute the integral of a product of two Legendre functions over
# a restricted domain
ip = ph.integ.pn1m1pn2m2(cltmin, cltmax, n1, m1, n2, m2, pnmj)


print(f'The integral of the product of two Legendre functions '
      f'n1 = {n1}, m1 = {m1}, n2 = {n2}, m2 = {m2} is {ip}')

PyHarm in single precision#

A few simple rules need to be obeyed when working with PyHarm in single precision.

Make sure you have compiled the library in single precision (see Installing).
In Python, import the PyHarm package as
```
>>> import pyharmf as phf
```

All floating point constants must be of the numpy.float32 data type. For instance, while in double precision you can write this:

>>> import pyharm as ph
>>> shcs = ph.shc.Shc.from_zeros(10, 1.0, 1.0)

in single precision, you have to write:

>>> import numpy as np
>>> import pyharmf as phf
>>> shcs = phf.shc.Shc.from_zeros(10, np.float32(1.0), np.float32(1.0))

All floating point numpy arrays must be of the numpy.float32 data type. So when creating a numpy array to be passed to PyHarm, you must set the dtype flag to numpy.float32:
```
>>> import numpy as np
>>> import pyharmf as phf
>>> x = np.zeros((5,), dtype=np.float32)
>>> phf.crd.PointSctr.from_arrays(x, x, x)
```

PyHarm

Contents