Frequency Decomposition of Zonal Winds - Spectral Analysis for Geophysical Data Cookbook

Overview¶

This notebook performs spectral analysis of zonal wind data at 850 hPa to decompose tropical atmospheric variability into different frequency bands using the Fast Fourier Transform (FFT).

Prerequisites¶

Concepts	Importance	Notes
Fourier analysis	Necessary

Time to learn: 20 minutes

Imports¶

import s3fs
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

Access data¶

We are accessing data stored on Project Pythia’s Jetstream2 storage.

URL = 'https://js2.jetstream-cloud.org:8001/' #Locate and read a file
fs = s3fs.S3FileSystem(anon=True, client_kwargs=dict(endpoint_url=URL))

uwind_ncep_ncar_store = s3fs.S3Map(
    root=f'pythia/uwind-ncep-ncar.zarr',
    s3=fs,
    check=False
)

uwind_ncep_ncar = xr.open_zarr(uwind_ncep_ncar_store)
uwind_ncep_ncar

Spectral analysis¶

Before carrying out the analysis, let’s take a look at the data on the first day:

projPC = ccrs.PlateCarree()
fig = plt.figure(figsize=(10, 5), layout='constrained')
ax = plt.subplot(111, projection=projPC)
uwind_ncep_ncar.uwnd.isel(level=0, time=0).plot(transform=projPC, cbar_kwargs={'orientation': 'horizontal'})
ax.scatter(90, 0, marker='x', color='k', s=50, transform=projPC, zorder=3)
ax.coastlines()
ax.gridlines(draw_labels=True, alpha=0.7)

<cartopy.mpl.gridliner.Gridliner at 0x7f0d7bebf770>

/home/runner/micromamba/envs/spectral-cookbook-dev/lib/python3.14/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

We will perform the FFT at one location first. Let’s select a location over the eastern Indian Ocean, indicated by the X in the above plot.

uwind_ei = uwind_ncep_ncar.uwnd.isel(level=0).sel(lat=0, lon=90)
uwind_ei.plot()

We need to center the data by removing its time mean.

centered_series = uwind_ei - uwind_ei.mean(dim='time')

centered_series.plot()

We now perform the FFT, also computing the explained variance in percent.

sampling_interval = 1  # Assuming daily data, so 1 day between samples

freqs = np.fft.fftfreq(len(centered_series), sampling_interval)
# Handle division by zero (set zero frequency to infinity)
periods = np.where(freqs != 0, 1 / np.abs(freqs), np.inf)

# Compute Fourier coefficients and power spectrum
fourier_coeffs = np.fft.fft(centered_series)
amplitude = np.abs(fourier_coeffs)
power = (amplitude ** 2)

# Normalize power spectrum by total power and scale by series variance
normalized_power = (power/np.sum(power)) * centered_series.var().compute().values
# Total explained variance of the centered series
exp_var = (normalized_power / centered_series.var().compute().values) * 100.

/tmp/ipykernel_4101/537958354.py:5: RuntimeWarning: divide by zero encountered in divide
  periods = np.where(freqs != 0, 1 / np.abs(freqs), np.inf)

plt.plot(periods,  exp_var, color='k')
plt.xscale('log')
plt.yscale('log')
plt.xlim(1e1, 1.5e4)
plt.ylim(1e-4, 3e0)
plt.axvline(365.25, color='r', linestyle='--', label='365 days', zorder=-1)
plt.axvline(365.25/2, color='b', linestyle='--', label=f'{365.25/2} days', zorder=-1)
plt.axvline(365.25/3, color='y', linestyle='--', label=f'{365.25/3} days', zorder=-1)
plt.axvline(365.25/4, color='g', linestyle='--', label=f'{365.25/4} days', zorder=-1)
plt.xlabel('Days')
plt.ylabel('Explained variance')
plt.legend()

There are noticeable spectral peaks at fractions of one year.

Now, we will perform the FFT at every location in the dataset. Then, we will compute the explained variance over various spectral bands associated with different kinds of variability.

Frequency Bands:

Annual Cycle (AC): 364-366 days
Intra-Seasonal (IS): 20-90 days (includes MJO timescales)
Inter-Annual (IA): 2-7 years (ENSO and longer-term variability)
Background (BG): >100 days (low-frequency variability)
Subseasonal (SS): <80 days (high-frequency variability)

In detail, the method is:

For each grid point, extract time series and center by removing the mean
Apply FFT to compute frequency spectrum and convert to period domain
Calculate normalized power spectrum scaled by variance
Apply frequency masks to isolate different bands
Compute variance and explained variance (%) for each band

As an example from the single location, we can compute the explained variance associated with the annual cycle, giving us about 12%:

# mask to extract different bands
ac_period = (364, 366)  # days
mask_AC = np.where((np.abs(periods) >= ac_period[0]) &
                    (periods <= ac_period[1]))

# Compute explained variances for each band
exp_var_AC = np.nansum(exp_var[mask_AC])
exp_var_AC

np.float32(0.1229185)

uwind_ncep_ncar.uwnd.shape

(17167, 2, 25, 144)

Set up variables and empty arrays:

data_array = uwind_ncep_ncar.isel(level=0).uwnd.values

# Get the shape of the data array
dtime, dlat, dlon = data_array.shape

# Reshape data to 2D array: (time, space) for easier processing
data_reshaped = data_array.reshape(dtime, dlat * dlon)

# Variance for different frequency bands
var_AC = np.zeros(dlat * dlon) * np.nan  # Annual Cycle
var_IS = np.zeros(dlat * dlon) * np.nan  # Intra-Seasonal Cycle
var_IA = np.zeros(dlat * dlon) * np.nan  # Inter-Annual Cycle
var_BG = np.zeros(dlat * dlon) * np.nan  # Background (longer than 100 days)
var_SS = np.zeros(dlat * dlon) * np.nan  # Subseasonal (shorter than 80 days)


# Explained variance for different frequency bands
exp_var_AC = np.zeros(dlat * dlon) * np.nan  # Annual Cycle
exp_var_IS = np.zeros(dlat * dlon) * np.nan  # Intra-Seasonal
exp_var_IA = np.zeros(dlat * dlon) * np.nan  # Inter-Annual Cycle

# Background (longer than 100 days)
exp_var_BG = np.zeros(dlat * dlon) * np.nan
# Subseasonal (shorter than 80 days)
exp_var_SS = np.zeros(dlat * dlon) * np.nan

# Annual, Intra-Seasonal, Inter-annual periods in days
ac_period = (364, 366)  # days
is_period = (20, 90)  # days
ia_period = (2*365.25, 7*365.25)  # days

# Background and Subseasonal cutoff periods
bg_cutoff = 100  # days
ss_cutoff = 80  # days

Here, we loop through all grid points and do the same analysis as before, storing the results in the empty arrays we just created.

for i in range(dlat * dlon):
    # Extract time series for current gridpoint
    series = data_reshaped[:, i]

    # Define sampling interval (daily data)
    sampling_interval = 1

    # Calculate mean and center the series
    mean_value = np.mean(series)
    centered_series = series - mean_value

    # Compute FFT frequencies and convert to periods
    freqs = np.fft.fftfreq(len(centered_series), sampling_interval)
    # Handle division by zero (set zero frequency to infinity)
    periods = np.where(freqs != 0, 1 / np.abs(freqs), np.inf)

    # Compute Fourier coefficients and power spectrum
    fourier_coeffs = np.fft.fft(centered_series)
    amplitud = np.abs(fourier_coeffs)
    power = (amplitud ** 2)

    # Normalize power spectrum by total power and scale by series variance
    normalized_power = (power/np.sum(power)) * np.var(centered_series)
    # Total explained variance of the centered series
    exp_var = (normalized_power / np.var(centered_series)) * 100.

    # mask to extract different bands
    mask_AC = np.where((np.abs(periods) >= ac_period[0]) &
                       (periods <= ac_period[1]))
    mask_IS = np.where((np.abs(periods) >= is_period[0]) &
                       (periods <= is_period[1]))
    mask_IA = np.where((np.abs(periods) >= ia_period[0]) &
                       (periods <= ia_period[1]))

    mask_BG = np.where(np.abs(periods) >= bg_cutoff)
    mask_SS = np.where(np.abs(periods) <= ss_cutoff)

    # Compute variances for each band
    var_AC[i] = np.nansum(normalized_power[mask_AC])
    var_IS[i] = np.nansum(normalized_power[mask_IS])
    var_IA[i] = np.nansum(normalized_power[mask_IA])
    var_BG[i] = np.nansum(normalized_power[mask_BG])
    var_SS[i] = np.nansum(normalized_power[mask_SS])

    # Compute explained variances for each band
    exp_var_AC[i] = np.nansum(exp_var[mask_AC])
    exp_var_IS[i] = np.nansum(exp_var[mask_IS])
    exp_var_IA[i] = np.nansum(exp_var[mask_IA])
    exp_var_BG[i] = np.nansum(exp_var[mask_BG])
    exp_var_SS[i] = np.nansum(exp_var[mask_SS])

/tmp/ipykernel_4101/2317136713.py:15: RuntimeWarning: divide by zero encountered in divide
  periods = np.where(freqs != 0, 1 / np.abs(freqs), np.inf)

Reshape the results back to (latitude, longitude) for plotting:

var_AC = var_AC.reshape(dlat, dlon)
var_IS = var_IS.reshape(dlat, dlon)
var_IA = var_IA.reshape(dlat, dlon)
var_BG = var_BG.reshape(dlat, dlon)
var_SS = var_SS.reshape(dlat, dlon)

exp_var_AC = exp_var_AC.reshape(dlat, dlon)
exp_var_IS = exp_var_IS.reshape(dlat, dlon)
exp_var_IA = exp_var_IA.reshape(dlat, dlon)
exp_var_BG = exp_var_BG.reshape(dlat, dlon)
exp_var_SS = exp_var_SS.reshape(dlat, dlon)

# Compute standard deviations
total_std = np.std(data_array, axis=0)
std_AC = np.sqrt(var_AC)
std_IS = np.sqrt(var_IS)
std_IA = np.sqrt(var_IA)
std_BG = np.sqrt(var_BG)
std_SS = np.sqrt(var_SS)

Plots¶

fig = plt.figure(figsize=(10, 5), layout='constrained')
ax = plt.subplot(111, projection=projPC)
plot = ax.pcolormesh(uwind_ncep_ncar.lon, uwind_ncep_ncar.lat, exp_var_AC, vmin=0, vmax=90)
ax.coastlines()
ax.gridlines(draw_labels=True, alpha=0.7)
plt.colorbar(plot, orientation='horizontal', label='Explained variance (%)')
plt.title('Annual Cycle')

fig = plt.figure(figsize=(10, 5), layout='constrained')
ax = plt.subplot(111, projection=projPC)
plot = ax.pcolormesh(uwind_ncep_ncar.lon, uwind_ncep_ncar.lat, exp_var_IS, vmin=0, vmax=90)
ax.coastlines()
ax.gridlines(draw_labels=True, alpha=0.7)
plt.colorbar(plot, orientation='horizontal', label='Explained variance (%)')
plt.title('Intra-Seasonal Variability [20-90 days]')

fig = plt.figure(figsize=(10, 5), layout='constrained')
ax = plt.subplot(111, projection=projPC)
plot = ax.pcolormesh(uwind_ncep_ncar.lon, uwind_ncep_ncar.lat, exp_var_IA, vmin=0, vmax=90)
ax.coastlines()
ax.gridlines(draw_labels=True, alpha=0.7)
plt.colorbar(plot, orientation='horizontal', label='Explained variance (%)')
plt.title('Interannual Variability [2-7 years]')

fig = plt.figure(figsize=(10, 5), layout='constrained')
ax = plt.subplot(111, projection=projPC)
plot = ax.pcolormesh(uwind_ncep_ncar.lon, uwind_ncep_ncar.lat, exp_var_BG, vmin=0, vmax=90)
ax.coastlines()
ax.gridlines(draw_labels=True, alpha=0.7)
plt.colorbar(plot, orientation='horizontal', label='Explained variance (%)')
plt.title('Background [>100 days]')

fig = plt.figure(figsize=(10, 5), layout='constrained')
ax = plt.subplot(111, projection=projPC)
plot = ax.pcolormesh(uwind_ncep_ncar.lon, uwind_ncep_ncar.lat, exp_var_SS, vmin=0, vmax=90)
ax.coastlines()
ax.gridlines(draw_labels=True, alpha=0.7)
plt.colorbar(plot, orientation='horizontal', label='Explained variance (%)')
plt.title('Subseasonal Variability [<80 days]')

Summary¶

In this notebook, we decomposed daily 850-hPa zonal winds into different frequency bands using the Fourier transform.

What’s next?¶

TBD

Resources and references¶

TBD