Rotated Maximum Covariance Analysis

Rotated Maximum Covariance Analysis (MCA) between two data sets.

# Load packages and data:
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from cartopy.crs import Orthographic, PlateCarree
from cartopy.feature import LAND
from matplotlib.gridspec import GridSpec

import xeofs as xe

Create 2 different DataArrays

t2m = xr.tutorial.load_dataset("air_temperature")["air"]
da1 = t2m.isel(lon=slice(0, 26))
da2 = t2m.isel(lon=slice(27, None))

Perform MCA

mca = xe.cross.MCA(n_modes=20, standardize=False, use_coslat=True)
mca.fit(da1, da2, dim="time")

<xeofs.cross.mca.MCA object at 0x7e858f640e50>

Apply Varimax-rotation to MCA solution

rot = xe.cross.MCARotator(n_modes=10)
rot.fit(mca)

<xeofs.cross.mca_rotator.MCARotator object at 0x7e858f643bd0>

Get rotated singular vectors, projections (PCs), homogeneous and heterogeneous patterns:

singular_vectors = rot.components()
scores = rot.scores()
hom_pats, pvals_hom = rot.homogeneous_patterns()
het_pats, pvals_het = rot.heterogeneous_patterns()

When two fields are expected, the output of the above methods is a list of length 2, with the first and second entry containing the relevant object for X and Y. For example, the p-values obtained from the two-sided t-test for the homogeneous patterns of X are:

pvals_hom[0]

<xarray.DataArray 'pvalues_of_left_homogeneous_patterns' (mode: 10, lat: 25,
                                                          lon: 26)>
1.562e-89 1.552e-68 7.788e-62 2.969e-52 ... 0.002105 0.005704 0.01526 0.03629
Coordinates:
  * lat      (lat) float32 15.0 17.5 20.0 22.5 25.0 ... 65.0 67.5 70.0 72.5 75.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 255.0 257.5 260.0 262.5
  * mode     (mode) int64 1 2 3 4 5 6 7 8 9 10
Attributes: (9)

xarray.DataArray

'pvalues_of_left_homogeneous_patterns'

• mode: 10
• lat: 25
• lon: 26

• 1.562e-89 1.552e-68 7.788e-62 2.969e-52 ... 0.005704 0.01526 0.03629

  array([[[1.56208256e-089, 1.55239986e-068, 7.78794124e-062, ...,
           2.23630076e-050, 8.37546668e-120, 6.12486970e-251],
          [2.13377869e-114, 9.51856072e-082, 3.11078791e-056, ...,
           6.63535708e-085, 1.64027702e-088, 1.04268512e-106],
          [7.81870839e-175, 1.21813874e-154, 4.64846387e-091, ...,
           1.38573825e-037, 4.80721743e-072, 1.33333492e-185],
          ...,
          [0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ...,
           0.00000000e+000, 0.00000000e+000, 0.00000000e+000],
          [0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ...,
           0.00000000e+000, 0.00000000e+000, 0.00000000e+000],
          [0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ...,
           0.00000000e+000, 0.00000000e+000, 0.00000000e+000]],
  
         [[0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ...,
           4.93958967e-127, 8.45513761e-042, 1.03197344e-010],
          [0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ...,
           1.73897067e-006, 4.63761210e-001, 7.64196234e-001],
          [2.34844968e-315, 5.10677495e-295, 0.00000000e+000, ...,
           3.29483746e-031, 3.24285441e-009, 2.47165248e-013],
  ...
          [1.08085510e-031, 4.58045008e-034, 3.60794214e-035, ...,
           7.71988815e-008, 1.69368171e-007, 8.68890142e-007],
          [1.50775155e-017, 3.35814264e-018, 3.12179389e-018, ...,
           1.32467333e-002, 2.10227389e-002, 3.66950212e-002],
          [1.51966734e-012, 1.69550861e-013, 4.50853230e-014, ...,
           7.46662192e-001, 7.71075026e-001, 8.17296768e-001]],
  
         [[1.23875224e-001, 3.38571130e-002, 1.78864650e-003, ...,
           6.87399789e-001, 7.72005690e-001, 5.05935071e-001],
          [2.49910224e-002, 1.60226695e-002, 9.86692885e-004, ...,
           8.65706271e-001, 2.93450126e-001, 5.20229461e-002],
          [9.82710346e-003, 8.71321224e-002, 1.49840305e-001, ...,
           2.52859669e-001, 9.57841679e-003, 8.40102447e-005],
          ...,
          [3.06865333e-009, 5.47748409e-010, 1.42061389e-010, ...,
           2.00371635e-004, 2.25804878e-003, 2.13782959e-002],
          [2.51266033e-002, 2.26673739e-002, 2.39401727e-002, ...,
           3.52103529e-004, 2.33547169e-003, 1.47653066e-002],
          [9.69410778e-001, 9.88052997e-001, 9.62623629e-001, ...,
           5.70368551e-003, 1.52586025e-002, 3.62908323e-002]]])

• Coordinates: (3)
  • lat
    (lat)
    float32
    15.0 17.5 20.0 ... 70.0 72.5 75.0

    array([15. , 17.5, 20. , 22.5, 25. , 27.5, 30. , 32.5, 35. , 37.5, 40. , 42.5,
           45. , 47.5, 50. , 52.5, 55. , 57.5, 60. , 62.5, 65. , 67.5, 70. , 72.5,
           75. ], dtype=float32)

  • lon
    (lon)
    float32
    200.0 202.5 205.0 ... 260.0 262.5

    array([200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
           225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
           250. , 252.5, 255. , 257.5, 260. , 262.5], dtype=float32)

  • mode
    (mode)
    int64
    1 2 3 4 5 6 7 8 9 10

    array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])

• Indexes: (3)
  • mode
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype='int64', name='mode'))

  • lat
    PandasIndex

    PandasIndex(Index([15.0, 17.5, 20.0, 22.5, 25.0, 27.5, 30.0, 32.5, 35.0, 37.5, 40.0, 42.5,
           45.0, 47.5, 50.0, 52.5, 55.0, 57.5, 60.0, 62.5, 65.0, 67.5, 70.0, 72.5,
           75.0],
          dtype='float32', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Index([200.0, 202.5, 205.0, 207.5, 210.0, 212.5, 215.0, 217.5, 220.0, 222.5,
           225.0, 227.5, 230.0, 232.5, 235.0, 237.5, 240.0, 242.5, 245.0, 247.5,
           250.0, 252.5, 255.0, 257.5, 260.0, 262.5],
          dtype='float32', name='lon'))

• Attributes: (9)

  model :

  Rotated MCA

  software :

  xeofs

  version :

  1.2.0

  date :

  2024-09-04 15:06:39

  n_modes :

  10

  power :

  1

  max_iter :

  1000

  rtol :

  1e-08

  compute :

  True

Create a mask to identifiy where p-values are below 0.05

hom_mask = [values < 0.05 for values in pvals_hom]
het_mask = [values < 0.05 for values in pvals_het]

Plot some relevant quantities of mode 2.

lonlats = [
    np.meshgrid(pvals_hom[0].lon.values, pvals_hom[0].lat.values),
    np.meshgrid(pvals_hom[1].lon.values, pvals_hom[1].lat.values),
]
proj = [
    Orthographic(central_latitude=30, central_longitude=-120),
    Orthographic(central_latitude=30, central_longitude=-60),
]
kwargs1 = {"cmap": "BrBG", "vmin": -0.05, "vmax": 0.05, "transform": PlateCarree()}
kwargs2 = {"cmap": "RdBu", "vmin": -1, "vmax": 1, "transform": PlateCarree()}

mode = 2

fig = plt.figure(figsize=(7, 14))
gs = GridSpec(5, 2)
ax1 = [fig.add_subplot(gs[0, i], projection=proj[i]) for i in range(2)]
ax2 = [fig.add_subplot(gs[1, i], projection=proj[i]) for i in range(2)]
ax3 = [fig.add_subplot(gs[2, i], projection=proj[i]) for i in range(2)]
ax4 = [fig.add_subplot(gs[3, i]) for i in range(2)]

for i, a in enumerate(ax1):
    singular_vectors[i].sel(mode=mode).plot(ax=a, **kwargs1)

for i, a in enumerate(ax2):
    hom_pats[i].sel(mode=mode).plot(ax=a, **kwargs2)
    a.scatter(
        lonlats[i][0],
        lonlats[i][1],
        hom_mask[i].sel(mode=mode).values * 0.5,
        color="k",
        alpha=0.5,
        transform=PlateCarree(),
    )
for i, a in enumerate(ax3):
    het_pats[i].sel(mode=mode).plot(ax=a, **kwargs2)
    a.scatter(
        lonlats[i][0],
        lonlats[i][1],
        het_mask[i].sel(mode=mode).values * 0.5,
        color="k",
        alpha=0.5,
        transform=PlateCarree(),
    )

for i, a in enumerate(ax4):
    scores[i].sel(mode=mode).plot(ax=a)
    a.set_xlabel("")


for a in np.ravel([ax1, ax2, ax3]):
    a.coastlines(color=".5")
    a.add_feature(LAND)

plt.tight_layout()
plt.savefig("rotated_mca.jpg")