PredictingThe Onset of the Indian Monsoon

Joint Math Meetings, Seattle, 2025

Sushovan Majhi

George Washington University, Washington D.C., USA

Today’s Talk

Background on the Indian Monsoon
Planet earth and TDA
Topology of a time-series
Predicting the onset and withdrawal of monsoon
Future work

The Indian Monsoon

Rain in India

The climate of India predominately depends on rainfall
Average: \(899\) mm with variation \(\pm 20\%\)
Monsoon is the typical rainy season
- Onset: May-July
- Withdrawal: Sep-Oct
What is the Monsoon?

The Monsoon System

Definition according to the Indian Meteorological Department (IMD)

The seasonal reversal of the direction of winds along the shores of the Indian Ocean, especially in the Arabian Sea, which blow from the southwest for half of the year and from the northeast for the other half.

Only partly understood
The most anticipated weather phenomenon
The biggest financial bet
- Affecting the Indian flora and fauna, economy, and agriculture
Notoriously difficult to predict

Monsoon is Coming!

The main objective of the project is develop an early warning system using topological data analysis (TDA).

Climate Science & TDA

Nuisances of Climate Data

high dimensionality
complexity of realistic models
presence of noise
missing values

Successful Applications of TDA

European Topsoil
- Savic, Toth, and Duponchel (2017)
Wildfire
- Kim and Vogel (2019)
Atmospheric River Patters
- Muszynski et al. (2019)
Weather Regimes
- Strommen et al. (2023)

Topology of a Time-series

Chaos

A chaotic system is deterministic but sensitive to initial conditions and manifests unpredictable patterns.

Lorenz System

\[ \begin{aligned} \dot{x} &= \sigma(y-x)\\ \dot{y} &= x(\rho-z)-y\\ \dot{z} &= xy-\beta z \end{aligned} \] Parameters: \(\sigma = 10\), \(\rho = 28\), \(\beta = 8 / 3\).

Chaos is still discernible when projected on the \(xy\)-plane.
How much of it still detectable using only a single signal, e.g. \(x\)?

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 18 },
  marginTop: 50,
  marginBottom: 50,
  height: 1000,
  width: width,
  x: { domain: [-20, 20], ticks: 15 },
  y: { domain: [-23, 30], ticks: 15 },
  marks: [
    Plot.line(lorenzData.slice(0, n), {
      x: "x",
      y: "y",
      strokeWidth: 2,
      stroke: "red"
    }),
    Plot.ruleX([-20]),
    Plot.ruleY([-23])
  ]
});

mutable play = false;
viewof playPause = Inputs.button(" ⏯️", {
  value: false,
  reduce: () => { if(mutable n >= 3999) mutable n = 0; mutable play = !play; } 
});

Time-Delay Embedding

Let’s now consider only one signal: \(x\)

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 18 },
  height: 250,
  width: width,
  x: {label: "t →"},
  y: { ticks: 4 },
  marks: [
    Plot.line(lorenzData.slice(500,600), {
      x: "n",
      y: "x"
    }),
    Plot.ruleX([w, w+delay], { stroke: "red" }),
    Plot.text([[w+0.2, -5], [w+delay+0.6, -5]], {text: [`d₁`, `d₂`]}),
    Plot.ruleX([500]),   
    Plot.ruleY([-12]),   
  ]
});

viewof w = Inputs.range([500, 600-delay-1], { value: 510, step: 1 });

The topology of the time-series is shown below:

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 18 },
  width: width,
  height: 420,
  x: {label: 'd₁', ticks: 15 },
  y: {label: 'd₂', ticks: 10 },
  marks: [
    Plot.ruleX([22.5]),
    Plot.ruleY([22]),
    Plot.dot(
      delayTS(
        lorenzData.slice(500,600).map((d) => d.z),
        delay
      ),
      {
        stroke: "white",
        r: (d, id) => id == w-500 ? 100 : 20,
        fill: (d, id) => id == w-500 ? "red" : "none"
      }
    )
  ]
});

viewof delay = Inputs.range([1, 50], { value: 1, label: tex`\text{Delay }(\tau)=`, step: 1 });

Sliding Window

We still consider only one signal: \(x\)

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 18 },
  height: 300,
  width: width,
  x: { ticks: 20, label: 't →' },
  y: {ticks: 7 },
  marks: [
    Plot.ruleX([500, 1300]),
    Plot.ruleY([-20]),
    Plot.line(lorenzData.slice(500), {
      x: "n",
      y: "x",
      stroke: (d) => d.n <= 1300 ? "green" : "red"
    })
  ]
});

Below are the time-delay embeddings of the above windows:

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 22 },
  width: width,
  height: 800,
  x: { ticks: 15, label: 'd₁'},
  y: { ticks: 15, label: 'd₂'},
  marks: [
    Plot.ruleX([-15]),
    Plot.ruleY([-15]),
    Plot.dot(
      delayTS(
        lorenzData.map((d) => d.x),
        delay1
      ).slice(500, 1300),
      {
        stroke: "green",
        r: 8
      }
    )
  ]
});

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 22 },
  width: width,
  height: 800,
  x: { ticks: 15, label: 'd₁' },
  y: { ticks: 15, label: 'd₂'},
  marks: [
    Plot.ruleX([-18]),
    Plot.ruleY([-18]),
    Plot.dot(
      delayTS(
        lorenzData.map((d) => d.x),
        delay1
      ).slice(1300),
      {
        stroke: "red",
        r: 8
      }
    )
  ]
});

viewof delay1 = Inputs.range([1, 50], { value: 1, label: tex`\text{Delay }(\tau)=`, step: 1 });

TDA Pipeline

We use overlapping windows to detect topological changes, e.g. transition to chaos.

take the input signal (time-series)
set a window size \(w\), delay \(\tau\), and delay coordinates \(k\)
position the window at the beginning of the time-series, and do the following:
- construct a point-cloud using time-delay embedding
- compute the persistence diagram (PD) of the point-cloud
- additionally, compute topological summaries like persistence landscapes, \(L^p\)-norms, etc
slide the window one-step forward, and repeat STEP 3 until the end of the time-series is reached.
In our case, we accumulate the landscape \(2\)-norms to output another time-series.

Predicting the Monsoon

The Monsoon Index

The daily average difference of rainfall across \(7\) different weather stations across India.

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 18 },
  width: width,
  height: 600,
  marginBottom: 50,
  y: { ticks: 10 },
  marks: [
    Plot.ruleY([-14]),
    Plot.line(
      data,
      {
        x: "Date",
        y: "rain",
      }
    )
  ]
});

data = FileAttachment("idx_1972.csv").csv({typed: true});
norms = FileAttachment("idx_1972_norms.csv").csv({typed: true});

Q: Can you guess the official onset from the signal?
A: June 18 was the official date.

Inference using the Sliding Window

Window size \(w=30\) and embedding dimension \(k=2\).

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 12 },
  width: width,
  height: 200,
  marks: [
    Plot.ruleY([-14]),
    Plot.dot(
      data,
      {
        x: "Date",
        y: "rain",
        stroke: (d, index) => index >= w1 && index <= w1 + w_size ? "red" : "white"
      }
    )
  ]
});

w_size = 30
viewof w1 = Inputs.range([0, data.length-w_size-1], { value: 0, step: 1 });

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 10 },
  width: width,
  height: 500,
  x: { domain: [-14,21], label: 'd₁'},
  y: { domain: [-14,21], label: 'd₂'},
  marks: [
    Plot.ruleX([-14]),
    Plot.ruleY([-14]),
    Plot.dot(
      delayTS(
        data.slice(w1, w1 + w_size).map((d) => d.rain),
        delay3
      ),
      {
        stroke: "red",
        r: 4,
      }
    )
  ]
});

viewof delay3 = Inputs.range([1, 20], { value: 1, label: tex`Delay`, step: 1 });

Persistence Pipeline

Go to embedding dimension: \(k=7\)
Window size: \(w=30\) (delay: \(\tau=1\))

The time-series of \(L^2\) norms of persistence landscapes:

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 13 },
  width: width,
  height: 500,
  marks: [
    Plot.ruleX([Date.parse('1972-01-01')]),
    Plot.ruleX([Date.parse('1972-12-31')], { stroke: 'red' }),
    Plot.ruleY([0]),
    dualAxisY(data, {y: "rain", anchor: "left", color: "white"}),
    dualAxisY(norms, {y: "norm", anchor: "right", color: "red"}),
    Plot.lineY(
      data,
      Plot.normalizeY("extent", {x: "Date", y: "rain", stroke: "white"})
    ),
    Plot.lineY(
      norms,
      Plot.normalizeY("extent", {x: "Date", y: "norm", stroke: "red"})
    )
  ]
});

The Onset and Withdrawal of Monsoon

Onset: first occurrence of norm larger than \(\mu+c_1*\sigma\)
Withdrawal: last occurrence of norm larger than \(\mu+c_2*\sigma\)

\(c_1=2\), \(c_2=1\):

Plot.plot({
  style: { fontFamily: 'monospace', fontSize: 13 },
  width: width,
  height: 500,
  marks: [
    Plot.text([[Date.parse('1972-06-15'), 1], [Date.parse('1972-09-26'), 1]], { text: ['June 15', 'Oct 26'] }),
    Plot.ruleY([0]),
    Plot.ruleX([Date.parse('1972-01-01')]),
    Plot.ruleX([Date.parse('1972-12-31')], { stroke: 'red'}),
    Plot.ruleX([Date.parse('1972-06-15')], { stroke: 'green', strokeWidth: 3 }),
    Plot.ruleX([Date.parse('1972-09-26')], { stroke: 'green', strokeWidth: 3}),
    dualAxisY(data, {y: "rain", anchor: "left", color: "white"}),
    dualAxisY(norms, {y: "norm", anchor: "right", color: "red"}),
    Plot.lineY(
      data,
      Plot.normalizeY("extent", {x: "Date", y: "rain", stroke: "white"})
    ),
    Plot.lineY(
      norms,
      Plot.normalizeY("extent", {x: "Date", y: "norm", stroke: "red"})
    )
  ]
});

Transition to Chaos

Comparing \(2\)-dimensional persistence diagrams

Future Directions

For a working early-warning system, we need a better resolution around the transition period
Enhance interpretability
Tune the parameters (e.g. window size \(w\), embedding dimension \(k\), delay \(k\)) to minimize false-negatives and false-positives
Instead of relying on time-delay embedding, use more signals
PCA to project the high-dimensional reconstruction of the phase space to gain insights into the attractors
Extend the study to other tropical monsoon systems.

Any Questions?

References

delayTS = function (ts, delay = 1, dim = 2) {
  return ts.map((d, i) => {
    const out = [];
    for (let j = 0; j < dim; j++) out.push(ts[i + j * delay]);
    return out;
  });
};
lorenz = function (
  [x, y, z] = [0, 0, 0],
  n = 100,
  { sigma = 10, r = 28, b = 8 / 3 } = {}
) {
  let F = function (t, [x, y, z]) {
    return [sigma * (y - x), r * x - x * z - y, x * y - b * z];
  };

  let s = new odex.Solver(3);
  let flow = [];

  for (let i = 0; i < n; i++) {
    [x, y, z] = s.solve(F, 0, [x, y, z], 0.01).y;
    flow.push({ n: i, x: x, y: y, z: z });
  }
  return flow;
};
lorenzData = lorenz([1, 1, 1], 4000);
odex = import("https://cdn.skypack.dev/odex@v2.0.4");
mutable n = 0;
{
  let i = 0;
  while (play && i < 4000) {
    yield Promises.delay(3, mutable n++);
  }
};

function dualAxisY(data, {y, ticks = 10, tickFormat, ...options} = {}) {
  const [y1, y2] = d3.extent(Plot.valueof(data, y));
  const scale = d3.scaleLinear().domain([y1, y2]);
  return Plot.axisY(d3.ticks(y1, y2, ticks), {...options, y: scale, tickFormat: scale.tickFormat(ticks, tickFormat)});
}

Kim, Hannah, and Christian Vogel. 2019. “Deciphering Active Wildfires in the Southwestern USA Using Topological Data Analysis.” Climate 7 (12): 135. https://doi.org/10.3390/cli7120135.

Muszynski, Grzegorz, Karthik Kashinath, Vitaliy Kurlin, Michael Wehner, and Prabhat. 2019. “Topological Data Analysis and Machine Learning for Recognizing Atmospheric River Patterns in Large Climate Datasets.” Geoscientific Model Development 12 (2): 613–28. https://doi.org/10.5194/gmd-12-613-2019.

Savic, Aleksandar, Gergely Toth, and Ludovic Duponchel. 2017. “Topological Data Analysis (TDA) Applied to Reveal Pedogenetic Principles of European Topsoil System.” Science of The Total Environment 586 (May): 1091–1100. https://doi.org/10.1016/j.scitotenv.2017.02.095.

Strommen, Kristian, Matthew Chantry, Joshua Dorrington, and Nina Otter. 2023. “A Topological Perspective on Weather Regimes.” Climate Dynamics 60 (5-6): 1415–45. https://doi.org/10.1007/s00382-022-06395-x.