Teleconnection Visualization Skill

Purpose

This skill provides best practices for visualizing climate teleconnections, correlation patterns, and atmospheric fields, specifically designed for reproducing figures from the Horel & Wallace (1981) study.

Core Visualization Types

1. Global Map Projections

Setup

python

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np

# Create figure with orthographic projection (like original paper)
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(central_longitude=180))

# Add features
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.OCEAN, facecolor='white')
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.gridlines(draw_labels=False, linewidth=0.5, alpha=0.5)

Plotting Spatial Patterns

python

def plot_spatial_pattern(data, ax=None, cmap='RdBu_r', vmin=None, vmax=None,
                         title='', contour_levels=None):
    """
    Plot spatial pattern on map.
    
    Args:
        data: xarray DataArray with lat/lon coordinates
        ax: Matplotlib axis with cartopy projection
        cmap: Colormap name
        vmin/vmax: Color scale limits
        title: Plot title
        contour_levels: Contour intervals (if None, uses contourf)
    """
    if ax is None:
        fig = plt.figure(figsize=(12, 8))
        ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    # Ensure longitude is in correct range
    lon = data.lon.values
    if lon.max() > 180:
        lon = ((lon + 180) % 360) - 180
        data = data.assign_coords(lon=lon).sortby('lon')
    
    # Set symmetric color scale for anomalies
    if vmin is None or vmax is None:
        vmax = np.abs(data).max().values
        vmin = -vmax
    
    # Plot
    if contour_levels is not None:
        im = ax.contourf(data.lon, data.lat, data,
                        levels=contour_levels,
                        transform=ccrs.PlateCarree(),
                        cmap=cmap, vmin=vmin, vmax=vmax,
                        extend='both')
        # Add contour lines
        cs = ax.contour(data.lon, data.lat, data,
                       levels=contour_levels,
                       colors='k', linewidths=0.5,
                       transform=ccrs.PlateCarree())
        ax.clabel(cs, inline=True, fontsize=8)
    else:
        im = data.plot(ax=ax, transform=ccrs.PlateCarree(),
                      cmap=cmap, vmin=vmin, vmax=vmax,
                      add_colorbar=False)
    
    # Add features
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
    ax.add_feature(cfeature.BORDERS, linewidth=0.3, linestyle=':')
    
    # Title
    ax.set_title(title, fontsize=12, fontweight='bold')
    
    # Colorbar
    cbar = plt.colorbar(im, ax=ax, orientation='horizontal',
                       pad=0.05, aspect=40, shrink=0.8)
    
    return im, cbar

# Example: Plot EOF pattern
fig = plt.figure(figsize=(14, 6))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=180))
plot_spatial_pattern(eof1, ax=ax, title='First EOF of SST Anomalies',
                    contour_levels=np.linspace(-0.3, 0.3, 13))
plt.savefig('figures/eof1_sst.png', dpi=300, bbox_inches='tight')

2. Correlation Maps

python

def plot_correlation_map(corr, pval, ax=None, significance_level=0.05,
                        title='Correlation Coefficients'):
    """
    Plot correlation coefficients with significance hatching.
    
    Args:
        corr: Correlation coefficient field (xarray DataArray)
        pval: P-value field
        ax: Matplotlib axis
        significance_level: Significance threshold
        title: Plot title
    """
    if ax is None:
        fig = plt.figure(figsize=(14, 8))
        ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    # Plot correlation
    im = corr.plot(ax=ax, transform=ccrs.PlateCarree(),
                  cmap='RdBu_r', vmin=-1, vmax=1,
                  add_colorbar=False)
    
    # Add significance stippling
    significant = pval < significance_level
    if significant.any():
        # Create stipple pattern
        lat_sig, lon_sig = np.meshgrid(corr.lat, corr.lon)
        ax.scatter(lon_sig[significant.T], lat_sig[significant.T],
                  s=1, c='k', alpha=0.3, transform=ccrs.PlateCarree(),
                  label=f'p < {significance_level}')
    
    # Add features
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
    ax.gridlines(draw_labels=True, linewidth=0.5, alpha=0.5)
    
    # Colorbar
    cbar = plt.colorbar(im, ax=ax, orientation='horizontal',
                       pad=0.05, label='Correlation Coefficient')
    
    ax.set_title(title, fontsize=12, fontweight='bold')
    
    return im

# Example
fig, ax = plt.subplots(1, 1, figsize=(14, 8),
                       subplot_kw={'projection': ccrs.PlateCarree()})
plot_correlation_map(corr_z700_sst, pval_z700_sst, ax=ax,
                    title='700 mb Height vs SST Index')

3. Time Series Plots

python

def plot_timeseries(data_dict, ylabel='', title='', figsize=(12, 4),
                   colors=None, linestyles=None, show_zero_line=True):
    """
    Plot multiple time series on same axes.
    
    Args:
        data_dict: Dictionary of {label: data_array}
        ylabel: Y-axis label
        title: Plot title
        figsize: Figure size
        colors: List of colors for each series
        linestyles: List of line styles
        show_zero_line: Add horizontal line at y=0
    """
    fig, ax = plt.subplots(1, 1, figsize=figsize)
    
    if colors is None:
        colors = plt.cm.tab10(np.linspace(0, 1, len(data_dict)))
    if linestyles is None:
        linestyles = ['-'] * len(data_dict)
    
    for (label, data), color, ls in zip(data_dict.items(), colors, linestyles):
        if hasattr(data, 'time'):
            ax.plot(data.time, data, label=label, color=color,
                   linestyle=ls, linewidth=1.5)
        else:
            ax.plot(data, label=label, color=color,
                   linestyle=ls, linewidth=1.5)
    
    if show_zero_line:
        ax.axhline(0, color='k', linestyle='--', linewidth=0.5, alpha=0.5)
    
    ax.set_ylabel(ylabel, fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.legend(frameon=False, fontsize=10)
    ax.grid(True, alpha=0.3)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    
    return fig, ax

# Example: Plot multiple indices
indices = {
    'SST Index': sst_index,
    'SOI': soi,
    '200 mb Index': z200_index
}
plot_timeseries(indices, ylabel='Normalized Anomaly',
               title='Climate Indices (1951-1978)')
plt.savefig('figures/climate_indices_timeseries.png', dpi=300, bbox_inches='tight')

4. Composite Difference Maps

python

def plot_composite_difference(pos_comp, neg_comp, difference, 
                              figsize=(18, 5), title_prefix=''):
    """
    Create three-panel plot: positive composite, negative composite, difference.
    
    Args:
        pos_comp: Positive episode composite
        neg_comp: Negative episode composite  
        difference: Difference (positive - negative)
        figsize: Figure size
        title_prefix: Prefix for panel titles
    """
    fig = plt.figure(figsize=figsize)
    
    # Positive composite
    ax1 = fig.add_subplot(1, 3, 1, projection=ccrs.PlateCarree())
    plot_spatial_pattern(pos_comp, ax=ax1, 
                        title=f'{title_prefix} Positive Episodes')
    
    # Negative composite
    ax2 = fig.add_subplot(1, 3, 2, projection=ccrs.PlateCarree())
    plot_spatial_pattern(neg_comp, ax=ax2,
                        title=f'{title_prefix} Negative Episodes')
    
    # Difference
    ax3 = fig.add_subplot(1, 3, 3, projection=ccrs.PlateCarree())
    plot_spatial_pattern(difference, ax=ax3,
                        title=f'{title_prefix} Difference')
    
    plt.tight_layout()
    return fig

# Example
fig = plot_composite_difference(pos_200mb, neg_200mb, diff_200mb,
                                title_prefix='200 mb Height')
plt.savefig('figures/composite_200mb.png', dpi=300, bbox_inches='tight')

5. Station Location Maps

python

def plot_station_locations(stations, ax=None, marker_size=100, 
                          add_labels=True, label_offset=(1, 1)):
    """
    Plot station locations on map.
    
    Args:
        stations: Dictionary of {name: (lat, lon)}
        ax: Matplotlib axis with cartopy projection
        marker_size: Size of station markers
        add_labels: Whether to add station name labels
        label_offset: (dx, dy) offset for labels in degrees
    """
    if ax is None:
        fig = plt.figure(figsize=(14, 8))
        ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    # Plot stations
    lats = [loc[0] for loc in stations.values()]
    lons = [loc[1] for loc in stations.values()]
    
    ax.scatter(lons, lats, s=marker_size, c='red', marker='*',
              edgecolors='black', linewidths=0.5,
              transform=ccrs.PlateCarree(), zorder=5)
    
    # Add labels
    if add_labels:
        for name, (lat, lon) in stations.items():
            ax.text(lon + label_offset[0], lat + label_offset[1], name,
                   fontsize=9, fontweight='bold',
                   transform=ccrs.PlateCarree(),
                   bbox=dict(boxstyle='round', facecolor='white', alpha=0.7))
    
    # Add features
    ax.add_feature(cfeature.LAND, facecolor='lightgray')
    ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
    ax.gridlines(draw_labels=True, linewidth=0.5, alpha=0.5)
    
    return ax

# Example
stations = {
    'Hilo': (19.7, -155.1),
    'Tahiti': (-17.5, -149.5),
    'Darwin': (-12.4, 130.9),
    'Nairobi': (-1.3, 36.8),
    'San Juan': (18.4, -66.0)
}

fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=180))
plot_station_locations(stations, ax=ax)
ax.set_title('Primary 200 mb Height Stations', fontsize=14, fontweight='bold')
plt.savefig('figures/station_locations.png', dpi=300, bbox_inches='tight')

6. Multi-Panel Time Series

python

def plot_station_timeseries(station_data, layout='vertical', figsize=None):
    """
    Create multi-panel plot of station time series.
    
    Args:
        station_data: Dictionary of {station_name: data_array}
        layout: 'vertical' or 'grid'
        figsize: Figure size
    """
    n_stations = len(station_data)
    
    if layout == 'vertical':
        nrows, ncols = n_stations, 1
        if figsize is None:
            figsize = (12, 2*n_stations)
    else:  # grid
        ncols = 2
        nrows = (n_stations + 1) // 2
        if figsize is None:
            figsize = (14, 3*nrows)
    
    fig, axes = plt.subplots(nrows, ncols, figsize=figsize, sharex=True)
    if nrows == 1 and ncols == 1:
        axes = [axes]
    else:
        axes = axes.flatten()
    
    for i, (station, data) in enumerate(station_data.items()):
        ax = axes[i]
        
        # Plot time series
        if hasattr(data, 'time'):
            ax.plot(data.time, data, 'k-', linewidth=1)
        else:
            ax.plot(data, 'k-', linewidth=1)
        
        # Zero line
        ax.axhline(0, color='r', linestyle='--', linewidth=0.5, alpha=0.5)
        
        # Formatting
        ax.set_ylabel('Normalized\nAnomaly', fontsize=9)
        ax.set_title(station, fontsize=10, fontweight='bold', loc='left')
        ax.grid(True, alpha=0.3)
        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
    
    # Remove extra axes
    for i in range(n_stations, len(axes)):
        axes[i].set_visible(False)
    
    # Common x-label
    axes[-1].set_xlabel('Year', fontsize=10)
    
    plt.tight_layout()
    return fig, axes

# Example
station_z200 = {
    'Nairobi': z200_nairobi,
    'Darwin': z200_darwin,
    'Hilo': z200_hilo,
    'Tahiti': z200_tahiti
}
fig, axes = plot_station_timeseries(station_z200, layout='grid')
plt.savefig('figures/station_200mb_timeseries.png', dpi=300, bbox_inches='tight')

7. Schematic Diagrams

python

def plot_teleconnection_schematic(centers_of_action, title=''):
    """
    Create schematic diagram showing teleconnection pattern.
    
    Args:
        centers_of_action: List of (lat, lon, sign, label) tuples
        title: Diagram title
    """
    fig = plt.figure(figsize=(14, 8))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.Orthographic(
        central_longitude=-100, central_latitude=45))
    
    # Add map features
    ax.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.5)
    ax.add_feature(cfeature.OCEAN, facecolor='lightblue', alpha=0.3)
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
    ax.gridlines(linewidth=0.5, alpha=0.3)
    
    # Plot centers of action
    for lat, lon, sign, label in centers_of_action:
        color = 'red' if sign > 0 else 'blue'
        marker = '^' if sign > 0 else 'v'
        
        ax.plot(lon, lat, marker=marker, markersize=15, color=color,
               markeredgecolor='black', markeredgewidth=1.5,
               transform=ccrs.PlateCarree(), zorder=5)
        
        # Add label
        ax.text(lon, lat-5, label, fontsize=10, fontweight='bold',
               ha='center', transform=ccrs.PlateCarree(),
               bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    
    # Draw arrows connecting centers
    for i in range(len(centers_of_action)-1):
        lat1, lon1 = centers_of_action[i][:2]
        lat2, lon2 = centers_of_action[i+1][:2]
        
        ax.plot([lon1, lon2], [lat1, lat2], 'k-', linewidth=2,
               alpha=0.5, transform=ccrs.Geodetic())
    
    ax.set_title(title, fontsize=14, fontweight='bold')
    
    return fig, ax

# Example: PNA pattern
pna_centers = [
    (20, -160, 1, 'H'),   # Subtropical Pacific high
    (45, -165, -1, 'L'),  # North Pacific low
    (55, -115, 1, 'H'),   # Western Canada high
    (30, -85, -1, 'L')    # SE US low
]
fig, ax = plot_teleconnection_schematic(pna_centers,
                                       title='Pacific/North American (PNA) Pattern')
plt.savefig('figures/pna_schematic.png', dpi=300, bbox_inches='tight')

Styling Best Practices

Color Schemes

python

# Diverging colormaps for anomalies (recommended)
diverging_cmaps = ['RdBu_r', 'BrBG', 'PiYG', 'PRGn', 'RdYlBu_r']

# Sequential colormaps for positive-definite fields
sequential_cmaps = ['YlOrRd', 'YlGnBu', 'Purples', 'Blues']

# Custom colormap for paper-like appearance
from matplotlib.colors import LinearSegmentedColormap

colors_neg = ['darkblue', 'blue', 'lightblue', 'white']
colors_pos = ['white', 'lightyellow', 'orange', 'darkred']
n_bins = 100

cmap_neg = LinearSegmentedColormap.from_list('neg', colors_neg, N=n_bins)
cmap_pos = LinearSegmentedColormap.from_list('pos', colors_pos, N=n_bins)

Professional Figure Formatting

python

# Set publication-quality defaults
plt.rcParams.update({
    'font.size': 10,
    'font.family': 'sans-serif',
    'font.sans-serif': ['Arial', 'DejaVu Sans'],
    'axes.labelsize': 11,
    'axes.titlesize': 12,
    'xtick.labelsize': 9,
    'ytick.labelsize': 9,
    'legend.fontsize': 9,
    'figure.dpi': 100,
    'savefig.dpi': 300,
    'savefig.bbox': 'tight',
    'axes.linewidth': 0.8,
    'grid.linewidth': 0.5,
    'lines.linewidth': 1.5,
})

# Contour levels matching original paper style
contour_levels = np.linspace(-1.0, 1.0, 21)  # 21 levels from -1 to 1

Annotation Best Practices

python

def add_figure_label(ax, label, loc='upper left'):
    """Add figure panel label (a), (b), etc."""
    ax.text(0.02, 0.98, label, transform=ax.transAxes,
           fontsize=12, fontweight='bold', va='top', ha='left',
           bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

def add_statistical_info(ax, corr_value, n_samples):
    """Add correlation and sample size information."""
    text = f'r = {corr_value:.2f}, n = {n_samples}'
    ax.text(0.98, 0.02, text, transform=ax.transAxes,
           fontsize=9, va='bottom', ha='right',
           bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

Creating Paper-Quality Figures

Figure 1: EOF Analysis

python

def create_figure1(eof_pattern, pc_timeseries, var_explained):
    """
    Reproduce Figure 1 from Horel & Wallace (1981):
    (a) First eigenvector spatial pattern
    (b) Time series of principal component
    """
    fig = plt.figure(figsize=(14, 10))
    
    # Panel (a): Spatial pattern
    ax1 = fig.add_subplot(2, 1, 1, projection=ccrs.PlateCarree(central_longitude=180))
    plot_spatial_pattern(eof_pattern, ax=ax1,
                        title=f'(a) First Eigenvector of Pacific SST ({var_explained:.1f}% variance)',
                        contour_levels=np.linspace(-0.05, 0.05, 11))
    add_figure_label(ax1, '(a)')
    
    # Panel (b): Time series
    ax2 = fig.add_subplot(2, 1, 2)
    ax2.plot(pc_timeseries.time, pc_timeseries, 'k-', linewidth=1.5)
    ax2.axhline(0, color='k', linestyle='--', linewidth=0.5)
    ax2.set_ylabel('Normalized Amplitude', fontsize=11)
    ax2.set_xlabel('Year', fontsize=11)
    ax2.set_title('(b) Time Series', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    ax2.spines['top'].set_visible(False)
    ax2.spines['right'].set_visible(False)
    add_figure_label(ax2, '(b)')
    
    plt.tight_layout()
    return fig

Figure 9: Correlation Panels

python

def create_figure9(corr_data_dict, pval_data_dict):
    """
    Reproduce Figure 9: Four-panel correlation maps.
    (a) 700 mb vs SST Index
    (b) 700 mb vs Fanning rainfall
    (c) 700 mb vs Southern Oscillation Index
    (d) 700 mb vs 200 mb Index
    """
    fig = plt.figure(figsize=(16, 12))
    
    panels = ['(a) SST Index', '(b) Fanning Rainfall',
              '(c) Southern Oscillation Index', '(d) 200 mb Index']
    
    for i, (title, corr_data, pval_data) in enumerate(
        zip(panels, corr_data_dict.values(), pval_data_dict.values())):
        
        ax = fig.add_subplot(2, 2, i+1, projection=ccrs.PlateCarree())
        plot_correlation_map(corr_data, pval_data, ax=ax, title=title)
        add_figure_label(ax, panels[i][:3])
    
    plt.tight_layout()
    return fig

Common Pitfalls

•Projection issues: Always specify transform=ccrs.PlateCarree() when plotting lat/lon data
•Color scale symmetry: For anomalies, always use symmetric color scales
•Missing data: Handle NaNs properly in contour plots
•Longitude wrap: Ensure longitude spans -180 to 180 or 0 to 360 consistently
•Resolution: Save figures at 300 DPI for publication quality

Export Guidelines

python

# Save in multiple formats
fig.savefig('figure1.png', dpi=300, bbox_inches='tight')
fig.savefig('figure1.pdf', bbox_inches='tight')  # Vector format
fig.savefig('figure1.eps', bbox_inches='tight')  # For LaTeX

Resources

•Cartopy documentation: https://scitools.org.uk/cartopy/
•Matplotlib gallery: https://matplotlib.org/stable/gallery/
•NCL graphics: https://www.ncl.ucar.edu/Applications/