amescap.Spectral_utils

Spectral_utils contains wave analysis routines. Note the dependencies on scipy.signal.

Third-party Requirements:

  • numpy

  • amescap.Script_utils

  • scipy.signal

  • ImportError

  • Exception

  • pyshtools (optional)

Module Contents

Functions

diurn_extract(VAR, N, tod, lon)

Extract the diurnal component of a field. Original code by John

extract_diurnal_harmonics(kmx, tmx, varIN, tod, lon)

Obtain west and east propagating waves. This is a Python

import_pyshtools()

Attempt to import pyshtools and set the global variable

reconstruct_diurn(amp, phas, tod, lon[, sumList])

Reconstructs a field wave based on its diurnal harmonics

zeroPhi_filter(VAR, btype, low_highcut, fs[, axis, ...])

A temporal filter that uses a forward and backward pass to prevent

zonal_construct(COEFFS_flat, VAR_shape[, btype, ...])

Recomposition into spherical harmonics

zonal_decomposition(VAR)

Decomposition into spherical harmonics. [A. Kling, 2020]
diurn_extract(VAR, N, tod, lon)

Extract the diurnal component of a field. Original code by John Wilson. Adapted by Alex Kling April, 2021

Parameters:

  • VAR (1D or ND array) – field to process. Time of day dimension must be first (e.g., [tod, time, lat, lon] or [tod]

  • N (int) – number of harmonics to extract (N=1 for diurnal, N=2 for diurnal AND semi diurnal, etc.)

  • tod (1D array) – universal time of day in sols (0-1.) If provided in hours (0-24), it will be normalized.

  • lon (1D array or float) – longitudes 0-360

Returns:

the amplitudes & phases for the Nth first harmonics, (e.g., size [Nh, time, lat, lon])

Return type:

ND arrays

extract_diurnal_harmonics(kmx, tmx, varIN, tod, lon)
Obtain west and east propagating waves. This is a Python

implementation of John Wilson’s space_time routine. Richard Urata 2025.

Parameters:

  • lon (1D array) – longitude [°] (0-360)

  • tod (1D array) – time_of_day [hour] (e.g., a 24 hour sol sampled every hour could be [0.5, 1.5, 2.5, ..., 23.5])

  • varIN (array) – variable for the Fourier analysis, expected to be scaled by zonal mean. First axis must be time and last axis must be lon (e.g., varIN[time, tod, lon], varIN[time, tod, lat, lon], or varIN[time, tod, lat, lon, lev])

  • kmx (int) – the number of longitudinal wavenumbers to extract (max = nlon/2)

  • tmx (int) – the number of tidal harmonics to extract (max = nsamples/2)

Returns:

(ampe) East propagating wave amplitude [same unit as varIN]; (ampw) West propagating wave amplitude [same unit as varIN]; (phasee) East propagating phase [°]; (phasew) West propagating phase [°]

Note

1. ampe, ampw, phasee, and phasew have dimensions [kmx, tmx] or [kmx, tmx, lat] or [kmx, tmx, lat, lev] etc.

2. The x and y axes may be constructed as follows, which will display the eastern and western modes:

klon = np.arange(0, kmx) # [wavenumber] [cycle/sol]
ktime = np.append(-np.arange(tmx, 0, -1), np.arange(0, tmx))
KTIME, KLON = np.meshgrid(ktime, klon)
amplitude = np.concatenate((ampw[:, ::-1], ampe), axis = 1)
phase = np.concatenate((phasew[:, ::-1], phasee), axis = 1)
import_pyshtools()

Attempt to import pyshtools and set the global variable PYSHTOOLS_AVAILABLE accordingly. Check if pyshtools is available and return its availability status

Returns:

True if pyshtools is available, False otherwise

Return type:

bool

reconstruct_diurn(amp, phas, tod, lon, sumList=[])

Reconstructs a field wave based on its diurnal harmonics

Parameters:

  • amp (array) – amplitude of the signal. Harmonics dimension FIRST (e.g., [N, time, lat, lon])

  • phas (array) – phase of the signal [hr]. Harmonics dimension FIRST

  • tod (1D array) – time of day in universal time [hr]

  • lon (1D array or float) – longitude for converting universal -> local time

  • sumList (list, optional) – the harmonics to include when reconstructing the wave (e.g., sumN=[1, 2, 4]), defaults to []

Returns:

a variable with reconstructed harmonics with N dimension FIRST and time of day SECOND ([N, tod, time, lat, lon]). If sumList is provided, the wave output harmonics will be aggregated (i.e., size = [tod, time, lat, lon])

Return type:

_type_

zeroPhi_filter(VAR, btype, low_highcut, fs, axis=0, order=4, add_trend=False)

A temporal filter that uses a forward and backward pass to prevent phase shift. Alex Kling 2020.

Parameters:

  • VAR (array) – values for filtering 1D or ND array. Filtered dimension must be FIRST. Adjusts axis as necessary.

  • btype (str) – filter type (i.e., “low”, “high” or “band”)

  • low_high_cut (int) – low, high, or [low, high] cutoff frequency depending on the filter [Hz or m-1]

  • fs (int) – sampling frequency [Hz or m-1]

  • axis (int) – if data is an ND array, this identifies the filtering dimension

  • order (int) – order for the filter

  • add_trend (bool) – if True, return the filtered output. If false, return the trend and filtered output.

Returns:

the filtered data

Note

Wn=[low, high] are expressed as a function of the Nyquist frequency

zonal_construct(COEFFS_flat, VAR_shape, btype=None, low_highcut=None)

Recomposition into spherical harmonics

Parameters:

  • COEFFS_flat (array) – Spherical harmonic coefficients as a flattened array, (e.g., [time, lat, lon] or [time x lev, 2, lat, lon])

  • VAR_shape (tuple) – shape of the original variable

  • btype (str or None) – filter type: “low”, “high”, or “band”. If None, returns reconstructed array using all zonal wavenumbers

  • low_high_cut (int or list) – low, high or [low, high] zonal wavenumber

Returns:

detrended variable reconstructed to original size (e.g., [time, lev, lat, lon])

Note

The minimum and maximum wavelenghts in [km] are computed:: dx = 2*np.pi * 3400 L_min = (1./kmax) * dx L_max = 1./max(kmin, 1.e-20) * dx if L_max > 1.e20:

L_max = np.inf

print(“(kmin,kmax) = ({kmin}, {kmax})

-> dx min = {L_min} km, dx max = {L_max} km”)

zonal_decomposition(VAR)

Decomposition into spherical harmonics. [A. Kling, 2020]

Parameters:

VAR – Detrend variable for decomposition. Lat is SECOND to LAST dimension and lon is LAST (e.g., [time,lat,lon] or [time,lev,lat,lon])

Returns:

(COEFFS) coefficient for harmonic decomposion, shape is flattened (e.g., [time, 2, lat/2, lat/2] [time x lev, 2, lat/2, lat/2]); (power_per_l) power spectral density, shape is re-organized (e.g., [time, lat/2] or [time, lev, lat/2])

Note

Output size is ([...,lat/2, lat/2]) as lat is the smallest dimension. This matches the Nyquist frequency.