Estimate Red-noise Spectrum

Note

Go to the end to download the full example code.

Estimate Red-noise Spectrum#

Before proceeding with all the steps, first import some necessary libraries and packages

import pooch
import xarray as xr
import numpy as np
import pandas as pd
import easyclimate as ecl
import matplotlib.pyplot as plt
from matplotlib import ticker

You can download following datasets here:

pooch.retrieve(
    "https://psl.noaa.gov/data/correlation/nina34.anom.data",
    known_hash=None,
    fname = "nina34.anom.data",
    path = ".",
)

Now we begin to read the txt file

# Read data from txt file
file_path = "nina34.anom.data"
with open(file_path, 'r') as file:
    lines = file.readlines()

# Parse data
years = []
values = []
for line in lines:
    # Skip empty lines or comment lines
    if not line.strip() or line.startswith('#'):
        continue
    parts = line.split()

    # Check if the first column is a year (integer)
    try:
        year = int(parts[0])  # Attempt to convert the first column to an integer
    except ValueError:
        print(f"Skipping invalid line (first column is not a year): {line.strip()}")
        continue

    # Check if there are 12 monthly values
    if len(parts[1:]) != 12:
        print(f"Skipping invalid line (missing 12 monthly values): {line.strip()}")
        continue

    # Convert the 12 monthly values to floats
    try:
        monthly_values = list(map(float, parts[1:]))
    except ValueError:
        print(f"Skipping invalid line (contains non-numeric data): {line.strip()}")
        continue

    years.append(year)
    values.append(monthly_values)

# Create time index
time = pd.date_range(start=f'{years[0]}-01', periods=len(years) * 12, freq='ME')

# Flatten the data into a 1D array
flat_values = np.array(values).flatten()

# Create xarray.DataArray
nino34 = xr.DataArray(
    flat_values,
    dims=['time'],
    coords={'time': time},
    attrs={'description': 'Nino 3.4 Index', 'units': 'Celsius'}
)

# Replace -99.99 with NaN
nino34 = nino34.where(nino34 != -99.99)

Skipping invalid line (missing 12 monthly values): 1948        2025
Skipping invalid line (first column is not a year): -99.99
Skipping invalid line (first column is not a year): Nino Anom 3.4 Index  using ersstv5 from CPC
Skipping invalid line (first column is not a year): https://psl.noaa.gov/data/climateindices/list/for info

Filter needed time range

nino34 = nino34.isel(time = slice(24,-12))
nino34

<xarray.DataArray (time: 900)> Size: 7kB
array([-1.99, -1.69, -1.42, -1.54, -1.75, -1.27, -1.01, -0.97, -0.98,
       -1.03, -1.23, -1.31, -1.3 , -1.04, -0.38, -0.23, -0.01,  0.  ,
        0.3 ,  0.17,  0.51,  0.49,  0.55,  0.31,  0.13, -0.01, -0.11,
       -0.02, -0.14, -0.54, -0.76, -0.56, -0.36, -0.46, -0.78, -0.39,
        0.2 ,  0.24,  0.29,  0.22,  0.35,  0.39,  0.14,  0.09,  0.29,
        0.16,  0.18,  0.41,  0.43,  0.28, -0.38, -1.17, -0.81, -0.93,
       -1.18, -1.43, -1.6 , -1.48, -1.13, -1.33, -0.93, -0.95, -1.06,
       -1.22, -1.28, -1.18, -1.15, -1.35, -1.43, -2.31, -2.45, -2.03,
       -1.2 , -0.99, -0.83, -0.97, -0.8 , -0.92, -1.06, -1.17, -0.98,
       -0.97, -1.14, -0.89, -0.51, -0.22,  0.18,  0.41,  0.62,  0.63,
        0.87,  0.84,  0.72,  0.7 ,  0.92,  1.3 ,  1.78,  1.49,  0.98,
        0.46,  0.38,  0.26,  0.03, -0.01, -0.32, -0.27,  0.05,  0.02,
        0.53,  0.43,  0.19,  0.06, -0.23, -0.35, -0.86, -0.77, -0.8 ,
       -0.48, -0.66, -0.42, -0.28, -0.46, -0.3 , -0.33, -0.25, -0.49,
       -0.41, -0.15, -0.28, -0.49, -0.45, -0.38, -0.31, -0.19, -0.35,
       -0.46, -0.18, -0.06, -0.41, -0.67, -0.94, -1.01, -0.63, -0.63,
       -0.59, -0.57, -0.48, -0.69, -0.89, -0.65, -0.53, -0.53, -0.78,
       -0.75, -0.95, -0.93, -0.78, -0.53, -0.1 , -0.03, -0.3 , -0.11,
        0.49,  0.62,  0.68,  0.64,  0.77,  1.02,  0.79,  0.38, -0.26,
       -0.87, -1.12, -1.14, -0.96, -1.26, -1.4 , -1.35, -1.44, -1.37,
...
        1.52,  1.25,  0.9 ,  0.38, -0.22, -0.69, -1.07, -1.39, -1.6 ,
       -1.69, -1.64, -1.6 , -1.54, -1.11, -0.93, -0.77, -0.52, -0.38,
       -0.43, -0.65, -0.8 , -1.05, -1.19, -1.06, -0.87, -0.67, -0.61,
       -0.5 , -0.32,  0.02,  0.25,  0.47,  0.38,  0.26,  0.16, -0.25,
       -0.53, -0.52, -0.25, -0.25, -0.4 , -0.42, -0.39, -0.38, -0.18,
       -0.2 , -0.14, -0.17, -0.49, -0.62, -0.28,  0.08,  0.32,  0.23,
       -0.06, -0.03,  0.29,  0.44,  0.75,  0.71,  0.51,  0.42,  0.47,
        0.7 ,  0.92,  1.18,  1.46,  1.93,  2.21,  2.36,  2.72,  2.66,
        2.57,  2.26,  1.62,  0.91,  0.3 , -0.03, -0.48, -0.58, -0.58,
       -0.74, -0.76, -0.5 , -0.43, -0.08,  0.03,  0.22,  0.37,  0.34,
        0.25, -0.16, -0.43, -0.56, -0.97, -0.98, -0.98, -0.78, -0.8 ,
       -0.51, -0.2 ,  0.04,  0.12,  0.09,  0.47,  0.9 ,  0.9 ,  0.89,
        0.65,  0.71,  0.81,  0.62,  0.55,  0.45,  0.35,  0.04,  0.03,
        0.48,  0.52,  0.52,  0.6 ,  0.37,  0.48,  0.36, -0.27, -0.34,
       -0.3 , -0.59, -0.83, -1.26, -1.42, -1.15, -1.  , -1.  , -0.8 ,
       -0.72, -0.46, -0.28, -0.39, -0.53, -0.55, -0.94, -0.94, -1.06,
       -0.94, -0.89, -0.97, -1.11, -1.11, -0.75, -0.69, -0.97, -1.07,
       -0.99, -0.9 , -0.85, -0.72, -0.46, -0.11,  0.14,  0.46,  0.84,
        1.02,  1.35,  1.6 ,  1.72,  2.02,  2.03,  1.81,  1.52,  1.13,
        0.77,  0.23,  0.17,  0.04, -0.12, -0.26, -0.27, -0.25, -0.6 ])
Coordinates:
  * time     (time) datetime64[ns] 7kB 1950-01-31 1950-02-28 ... 2024-12-31
Attributes:
    description:  Nino 3.4 Index
    units:        Celsius


Calculate red noise

result_redfit = ecl.filter.calc_redfit(nino34)
result_redfit = result_redfit.assign_coords({"period_month": (result_redfit.period)/12})
result_redfit

<xarray.Dataset> Size: 22kB
Dimensions:       (freq: 451)
Coordinates:
  * freq          (freq) float32 2kB 0.0 0.001111 0.002222 ... 0.4978 0.4989 0.5
    period        (freq) float32 2kB inf 900.0 450.0 300.0 ... 2.009 2.004 2.0
    period_month  (freq) float32 2kB inf 75.0 37.5 25.0 ... 0.1674 0.167 0.1667
Data variables:
    gxx           (freq) float32 2kB 1.603e-13 2.266 5.423 ... 0.01003 0.008833
    gxx_corr      (freq) float32 2kB 4.819 3.171 5.445 ... 0.008368 0.007121
    gred_th       (freq) float32 2kB 66.46 64.73 60.06 ... 0.03038 0.03038
    gred          (freq) float32 2kB 2.211e-12 46.25 59.82 ... 0.03643 0.03769
    corrFac       (freq) float32 2kB 3.326e-14 0.7145 0.9961 ... 1.199 1.24
    chi2_80       (freq) float32 2kB 107.0 104.2 96.68 ... 0.04891 0.04891
    chi2_90       (freq) float32 2kB 153.0 149.0 138.3 ... 0.06995 0.06995
    chi2_95       (freq) float32 2kB 199.1 194.0 180.0 ... 0.09105 0.09104
    chi2_99       (freq) float32 2kB 306.0 298.1 276.5 ... 0.1399 0.1399 0.1399
Attributes:
    Description:                                Estimating red-noise spectra ...
    Input:                                      OFAC = 1.0, HIFAC = 1.0, n50 ...
    Initial values:                             idum = -8.5992514e+32, Data v...
    Results:                                    Avg. autocorr. coeff., rho = ...
    Equality of theoretical and data spectrum:  90-% acceptance region: rcrit...
    Elapsed time:                               0 [s]
    About:                                      Michael Schulz, Manfred Mudel...
    Reference:                                  Schulz, M. and Mudelsee, M. (...
    Python platform:                            Shen yulu => https://github.c...


Draw the red noise graph

fig, ax = plt.subplots()

result_redfit.gxx.plot(ax = ax, x = 'period_month', color = 'black')
result_redfit.chi2_95.plot(ax = ax, ls = '--', x = 'period_month', label = '95% CI', color = 'red')
result_redfit.chi2_90.plot(ax = ax, ls = '--', x = 'period_month', label = '90% CI', color = 'blue')
ax.legend()

ax.set_xscale('log', base = 2, subs = None)
ax.xaxis.set_major_formatter(ticker.ScalarFormatter())

ax.set_xlim(2, 16)

ax.set_xlabel('Period (years)')
ax.set_ylabel('Spectral Amplitude')
plot redfit
Text(38.347222222222214, 0.5, 'Spectral Amplitude')

Total running time of the script: (0 minutes 0.633 seconds)