Case Study 2 — Climate Data as Tensors
The Pacific Climate Research Consortium's Temperature Analysis
Context
The Pacific Climate Research Consortium (PCRC) maintains gridded temperature records covering the Pacific basin. The data is naturally represented as a 3-dimensional tensor: latitude $\times$ longitude $\times$ time. Each entry $\mathcal{T}_{i,j,t}$ records the monthly mean temperature at grid point $(i, j)$ for month $t$.
Traditional analysis (plotting time series at individual locations, computing spatial averages) discards structure. Tensor decomposition methods extract the dominant spatial patterns and their temporal evolution simultaneously — revealing phenomena like El Nino/Southern Oscillation (ENSO) patterns, seasonal cycles, and long-term warming trends.
This case study constructs a synthetic but physically realistic climate tensor, applies PCA (via matrix unfolding) and CP decomposition, and interprets the extracted patterns.
Constructing the Climate Tensor
We simulate a temperature field with three embedded physical modes: a seasonal cycle, an ENSO-like pattern, and a long-term warming trend.
import numpy as np
def generate_climate_tensor(
n_lat: int = 30,
n_lon: int = 60,
n_months: int = 360,
seed: int = 42,
) -> dict:
"""Generate a synthetic climate tensor with known spatial-temporal patterns.
Three physical modes are embedded:
1. Seasonal cycle: varies with latitude, annual period
2. ENSO-like pattern: concentrated in equatorial Pacific, ~4-year period
3. Warming trend: global, with polar amplification
Args:
n_lat: Number of latitude grid points.
n_lon: Number of longitude grid points.
n_months: Number of monthly time steps.
seed: Random seed.
Returns:
Dictionary with the tensor and component information.
"""
rng = np.random.default_rng(seed)
# Spatial grids (Pacific basin: 30S-30N, 120E-80W)
lat = np.linspace(-30, 30, n_lat)
lon = np.linspace(120, 280, n_lon)
lat_grid, lon_grid = np.meshgrid(lat, lon, indexing="ij")
time_months = np.arange(n_months)
# --- Mode 1: Seasonal cycle ---
# Spatial: depends on latitude (stronger at higher latitudes)
spatial_seasonal = np.cos(np.radians(lat_grid)) * 0.5 + 0.5
spatial_seasonal = spatial_seasonal * np.abs(np.sin(np.radians(lat_grid)))
# Temporal: annual cycle
temporal_seasonal = 10.0 * np.cos(2 * np.pi * time_months / 12)
mode_seasonal = np.einsum("ij,t->ijt",
spatial_seasonal / np.linalg.norm(spatial_seasonal),
temporal_seasonal)
# --- Mode 2: ENSO-like pattern ---
# Spatial: concentrated near equator, strongest in central-eastern Pacific
spatial_enso = np.exp(-(lat_grid**2) / (2 * 8**2)) # Equatorial
spatial_enso *= np.exp(-((lon_grid - 220)**2) / (2 * 40**2)) # Central Pacific
# Temporal: irregular oscillation with ~4-year period
temporal_enso = 5.0 * np.sin(2 * np.pi * time_months / 48
+ 0.3 * np.sin(2 * np.pi * time_months / 120))
mode_enso = np.einsum("ij,t->ijt",
spatial_enso / np.linalg.norm(spatial_enso),
temporal_enso)
# --- Mode 3: Warming trend ---
# Spatial: global but with polar amplification
spatial_warming = 1.0 + 0.5 * np.abs(np.sin(np.radians(lat_grid)))
# Temporal: linear trend
temporal_warming = 0.05 * time_months / 12 # ~0.05°C per year
mode_warming = np.einsum("ij,t->ijt",
spatial_warming / np.linalg.norm(spatial_warming),
temporal_warming)
# --- Combined tensor ---
# Climatological mean
mean_temp = 25.0 - 15.0 * np.abs(np.sin(np.radians(lat_grid)))
T_mean = np.repeat(mean_temp[:, :, np.newaxis], n_months, axis=2)
# Full tensor: mean + modes + noise
noise = 0.5 * rng.standard_normal((n_lat, n_lon, n_months))
T = T_mean + mode_seasonal + mode_enso + mode_warming + noise
return {
"tensor": T,
"lat": lat,
"lon": lon,
"time_months": time_months,
"true_modes": {
"seasonal": {"spatial": spatial_seasonal, "temporal": temporal_seasonal},
"enso": {"spatial": spatial_enso, "temporal": temporal_enso},
"warming": {"spatial": spatial_warming, "temporal": temporal_warming},
},
"shape": T.shape,
}
climate = generate_climate_tensor()
T = climate["tensor"]
print(f"Climate tensor shape: {T.shape}")
print(f" {T.shape[0]} latitude points × {T.shape[1]} longitude points × {T.shape[2]} months")
print(f" Total elements: {T.size:,}")
print(f" Temperature range: [{T.min():.1f}, {T.max():.1f}] °C")
Tensor Unfolding and PCA
A tensor can be "unfolded" (matricized) along each mode to produce a matrix. PCA on the spatial unfolding reveals the dominant spatial patterns; PCA on the temporal unfolding reveals the dominant temporal patterns.
def unfold_tensor(T: np.ndarray, mode: int) -> np.ndarray:
"""Unfold (matricize) a 3-tensor along the specified mode.
Mode 0: (n_lat, n_lon * n_time) — each row is a latitude's full history
Mode 1: (n_lon, n_lat * n_time) — each row is a longitude's full history
Mode 2: (n_time, n_lat * n_lon) — each row is a spatial field at one time
Args:
T: 3-tensor of shape (n_lat, n_lon, n_time).
mode: Mode along which to unfold (0, 1, or 2).
Returns:
Unfolded matrix.
"""
return np.reshape(np.moveaxis(T, mode, 0), (T.shape[mode], -1))
def spatial_pca(T: np.ndarray, n_components: int = 5) -> dict:
"""Perform PCA on the temporal unfolding to extract spatial patterns.
Unfolding along mode 2 gives a matrix where each row is a flattened
spatial field at one time step. The principal components are spatial
patterns (EOFs in climate science terminology).
Args:
T: Climate tensor (n_lat, n_lon, n_time).
n_components: Number of components to extract.
Returns:
Dictionary with spatial patterns (EOFs) and temporal coefficients (PCs).
"""
n_lat, n_lon, n_time = T.shape
# Unfold: (n_time, n_lat * n_lon)
X = unfold_tensor(T, mode=2)
# Center (remove temporal mean at each grid point)
X_mean = X.mean(axis=0)
X_centered = X - X_mean
# SVD
U, sigma, Vt = np.linalg.svd(X_centered, full_matrices=False)
# Spatial patterns (EOFs): reshape rows of Vt back to spatial grids
eofs = Vt[:n_components, :].reshape(n_components, n_lat, n_lon)
# Temporal coefficients (PCs): columns of U scaled by singular values
pcs = U[:, :n_components] * sigma[:n_components]
# Explained variance
total_var = np.sum(sigma**2)
explained_ratio = sigma[:n_components]**2 / total_var
return {
"eofs": eofs,
"pcs": pcs,
"explained_variance_ratio": explained_ratio,
"cumulative_variance": np.cumsum(explained_ratio),
"singular_values": sigma[:n_components],
"mean_field": X_mean.reshape(n_lat, n_lon),
}
result = spatial_pca(T, n_components=5)
print("Spatial PCA Results")
print("=" * 40)
for i in range(5):
print(f"EOF {i+1}: {result['explained_variance_ratio'][i]*100:.1f}% "
f"(cumulative: {result['cumulative_variance'][i]*100:.1f}%)")
In climate science, these spatial patterns are called Empirical Orthogonal Functions (EOFs) and the temporal coefficients are called Principal Components (PCs). The terminology differs from ML convention, but the mathematics is identical: EOF analysis is PCA on the temporal unfolding of the data tensor.
Interpreting the Extracted Modes
def interpret_modes(
result: dict,
climate: dict,
n_modes: int = 3,
) -> None:
"""Compare extracted EOFs/PCs to the known true modes.
This is possible because we generated the data with known structure.
In practice, interpretation relies on domain knowledge.
Args:
result: Output of spatial_pca.
climate: Output of generate_climate_tensor.
n_modes: Number of modes to interpret.
"""
true_modes = climate["true_modes"]
for i in range(n_modes):
eof = result["eofs"][i]
pc = result["pcs"][:, i]
# Correlate with known spatial patterns
correlations = {}
for name, mode in true_modes.items():
spatial = mode["spatial"]
# Normalize for correlation
s_norm = spatial / np.linalg.norm(spatial)
e_norm = eof / np.linalg.norm(eof)
correlations[name] = abs(np.sum(s_norm * e_norm))
best_match = max(correlations, key=correlations.get)
print(f"\nEOF {i+1} (explains {result['explained_variance_ratio'][i]*100:.1f}%)")
print(f" Best match: {best_match} "
f"(spatial correlation: {correlations[best_match]:.3f})")
print(f" All correlations: "
+ ", ".join(f"{k}: {v:.3f}" for k, v in correlations.items()))
# Temporal correlation with best match
true_temporal = true_modes[best_match]["temporal"]
t_corr = abs(np.corrcoef(pc, true_temporal[:len(pc)])[0, 1])
print(f" Temporal correlation: {t_corr:.3f}")
interpret_modes(result, climate)
With the synthetic data, the extracted modes should closely match the embedded patterns:
- EOF 1 corresponds to the seasonal cycle (largest signal, latitude-dependent)
- EOF 2 corresponds to the ENSO-like pattern (equatorial Pacific concentration)
- EOF 3 corresponds to the warming trend (global, with polar amplification)
CP Decomposition
The CANDECOMP/PARAFAC (CP) decomposition represents a tensor as a sum of rank-1 tensors:
$$\mathcal{T} \approx \sum_{r=1}^{R} \lambda_r \, \mathbf{a}_r \circ \mathbf{b}_r \circ \mathbf{c}_r$$
where $\circ$ denotes the outer product. Each component is defined by three vectors: one for each mode (latitude, longitude, time). Unlike matrix SVD, the CP decomposition does not have an exact closed-form solution — it is computed iteratively via alternating least squares (ALS).
def cp_als(
T: np.ndarray,
rank: int,
max_iter: int = 100,
tol: float = 1e-6,
seed: int = 42,
) -> dict:
"""Compute CP decomposition via alternating least squares.
Args:
T: 3-tensor of shape (I, J, K).
rank: Number of components.
max_iter: Maximum ALS iterations.
tol: Convergence tolerance (relative change in reconstruction error).
seed: Random seed for initialization.
Returns:
Dictionary with factor matrices and reconstruction info.
"""
rng = np.random.default_rng(seed)
I, J, K = T.shape
# Initialize factor matrices randomly
A = rng.standard_normal((I, rank)) # Latitude factors
B = rng.standard_normal((J, rank)) # Longitude factors
C = rng.standard_normal((K, rank)) # Time factors
# Normalize
for r in range(rank):
A[:, r] /= np.linalg.norm(A[:, r])
B[:, r] /= np.linalg.norm(B[:, r])
C[:, r] /= np.linalg.norm(C[:, r])
lambdas = np.ones(rank)
prev_error = np.inf
for iteration in range(max_iter):
# Update A: unfold mode 0, solve A @ (C ⊙ B)^T ≈ T_(0)
# where ⊙ is the Khatri-Rao (columnwise Kronecker) product
CB = np.zeros((J * K, rank))
for r in range(rank):
CB[:, r] = np.kron(C[:, r], B[:, r])
T_0 = unfold_tensor(T, 0) # (I, J*K)
A = T_0 @ CB @ np.linalg.inv(CB.T @ CB + 1e-10 * np.eye(rank))
# Normalize columns of A
for r in range(rank):
norm = np.linalg.norm(A[:, r])
if norm > 0:
lambdas[r] = norm
A[:, r] /= norm
# Update B: unfold mode 1
CA = np.zeros((I * K, rank))
for r in range(rank):
CA[:, r] = np.kron(C[:, r], A[:, r])
T_1 = unfold_tensor(T, 1) # (J, I*K)
B = T_1 @ CA @ np.linalg.inv(CA.T @ CA + 1e-10 * np.eye(rank))
for r in range(rank):
norm = np.linalg.norm(B[:, r])
if norm > 0:
lambdas[r] *= norm
B[:, r] /= norm
# Update C: unfold mode 2
AB = np.zeros((I * J, rank))
for r in range(rank):
AB[:, r] = np.kron(B[:, r], A[:, r])
T_2 = unfold_tensor(T, 2) # (K, I*J)
C = T_2 @ AB @ np.linalg.inv(AB.T @ AB + 1e-10 * np.eye(rank))
for r in range(rank):
norm = np.linalg.norm(C[:, r])
if norm > 0:
lambdas[r] *= norm
C[:, r] /= norm
# Reconstruction error
T_approx = np.zeros_like(T)
for r in range(rank):
T_approx += lambdas[r] * np.einsum("i,j,k->ijk", A[:, r], B[:, r], C[:, r])
error = np.linalg.norm(T - T_approx) / np.linalg.norm(T)
if abs(prev_error - error) / max(prev_error, 1e-10) < tol:
print(f"Converged at iteration {iteration + 1}, "
f"relative error: {error:.6f}")
break
prev_error = error
# Sort by importance
order = np.argsort(-np.abs(lambdas))
lambdas = lambdas[order]
A = A[:, order]
B = B[:, order]
C = C[:, order]
return {
"lambdas": lambdas,
"A": A, # Latitude factors
"B": B, # Longitude factors
"C": C, # Time factors
"n_iter": iteration + 1,
"relative_error": error,
}
# Remove the mean field before decomposition
T_anomaly = T - T.mean(axis=2, keepdims=True)
cp_result = cp_als(T_anomaly, rank=5, max_iter=200)
print(f"\nCP Decomposition Results")
print(f"Relative reconstruction error: {cp_result['relative_error']:.4f}")
print(f"Component weights: {cp_result['lambdas'].round(1)}")
Connecting CP Decomposition to PCA
The CP decomposition and spatial PCA are related but not identical:
- PCA (via matrix unfolding) treats the tensor as a matrix by collapsing two modes. The spatial EOFs are orthogonal by construction, but they may mix physical phenomena that have correlated spatial patterns.
- CP decomposition maintains the three-mode structure. Each component has separate latitude, longitude, and time factors. The components are not orthogonal, but they can better separate physically distinct phenomena.
def compare_methods(
pca_result: dict,
cp_result: dict,
climate: dict,
) -> None:
"""Compare PCA and CP decomposition results.
Args:
pca_result: Output of spatial_pca.
cp_result: Output of cp_als.
climate: Original climate data with true modes.
"""
n_lat, n_lon = climate["tensor"].shape[:2]
print("\n=== Method Comparison ===\n")
print("PCA explained variance (top 5):")
for i in range(min(5, len(pca_result['explained_variance_ratio']))):
print(f" EOF {i+1}: {pca_result['explained_variance_ratio'][i]*100:.1f}%")
print(f"\nCP relative reconstruction error: {cp_result['relative_error']:.4f}")
print(f"CP component weights: {cp_result['lambdas'][:5].round(1)}")
# Reconstruct spatial patterns from CP factors
print("\nCP spatial patterns (A ⊗ B for each component):")
for r in range(min(3, len(cp_result['lambdas']))):
spatial_cp = np.outer(cp_result['A'][:, r], cp_result['B'][:, r])
# Compare to known modes
for name, mode in climate['true_modes'].items():
sp = mode['spatial']
corr = abs(np.sum(
(spatial_cp / np.linalg.norm(spatial_cp)) *
(sp / np.linalg.norm(sp))
))
if corr > 0.3:
print(f" Component {r+1}: matches '{name}' "
f"(correlation: {corr:.3f})")
compare_methods(result, cp_result, climate)
Key Takeaways
-
Tensors are the natural representation for spatiotemporal data. Collapsing a 3D climate field into a 2D table discards structural information. Tensor methods preserve the multimodal nature of the data.
-
Matrix unfolding + PCA is the simplest approach. It reduces the tensor to a matrix problem (solved by SVD) and produces orthogonal spatial patterns. Climate scientists have used this approach (Empirical Orthogonal Functions) since the 1950s.
-
CP decomposition separates modes along each dimension. Instead of a single spatial pattern, each component has separate latitude and longitude factors. This can better capture phenomena with separable spatial structure (e.g., ENSO, which varies smoothly in both latitude and longitude).
-
The connection to linear algebra is direct. The spatial modes extracted by PCA are eigenvectors of the spatial covariance matrix $\mathbf{C} = \frac{1}{T}\mathbf{X}^\top \mathbf{X}$, where $\mathbf{X}$ is the temporal unfolding. The singular values quantify each mode's importance. The first few modes capture physically meaningful phenomena; the remaining modes capture noise.
-
This analysis extends to deep learning. In Part II, we will use convolutional neural networks to learn spatial patterns from climate data (Chapter 8), and the latent representations learned by these networks serve a similar role to the EOFs derived here — but with the capacity to capture nonlinear spatial patterns that PCA cannot.
Case Study 2 for Chapter 1 of Advanced Data Science