Rotated Maximum Covariance Analysis

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

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

from xeofs.models import MCA, MCARotator

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 = MCA(n_modes=20, standardize=False, use_coslat=True)
mca.fit(da1, da2, dim="time")

<xeofs.models.mca.MCA object at 0x7f016048d2d0>

Apply Varimax-rotation to MCA solution

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

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)>
array([[[5.91798633e-081, 4.09037030e-064, 6.23425561e-058, ...,
         4.45216591e-046, 7.99246088e-107, 6.95533368e-234],
        [9.61231974e-109, 4.30909052e-079, 2.41956104e-054, ...,
         1.13316213e-079, 1.17033263e-086, 1.58776474e-105],
        [3.65945612e-169, 1.03636735e-150, 5.85938600e-088, ...,
         2.66091078e-037, 1.39445455e-071, 1.89947470e-183],
        ...,
        [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, ...,
         2.15049395e-115, 1.25690777e-037, 3.79422747e-010],
        [0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ...,
         3.68434950e-006, 4.54442588e-001, 7.47361402e-001],
        [5.82970051e-304, 1.04446785e-286, 1.67982320e-322, ...,
         3.12415723e-031, 2.97910617e-009, 3.83983102e-013],
...
        [3.81951811e-033, 1.24766309e-035, 9.83759118e-037, ...,
         1.08208863e-007, 2.09241618e-007, 9.75026228e-007],
        [7.63507975e-018, 1.84998265e-018, 1.85005250e-018, ...,
         1.47665784e-002, 2.23166837e-002, 3.74986857e-002],
        [5.23820545e-013, 5.59326448e-014, 1.46539526e-014, ...,
         7.79688798e-001, 7.97072146e-001, 8.38057312e-001]],

       [[1.33009336e-001, 3.88106643e-002, 2.39386496e-003, ...,
         9.33561796e-001, 6.41357062e-001, 5.12874473e-001],
        [2.47971905e-002, 1.49401803e-002, 1.01803200e-003, ...,
         9.48611093e-001, 3.16366009e-001, 6.29374852e-002],
        [1.00284803e-002, 8.50825810e-002, 1.51720312e-001, ...,
         2.62222025e-001, 1.27842313e-002, 1.43419078e-004],
        ...,
        [6.95159651e-010, 1.10849674e-010, 2.79734189e-011, ...,
         2.20338645e-004, 2.58483821e-003, 2.44608586e-002],
        [2.08872724e-002, 1.93770600e-002, 2.08168263e-002, ...,
         3.74564462e-004, 2.56337651e-003, 1.63483691e-002],
        [9.06082861e-001, 8.64987130e-001, 8.39176635e-001, ...,
         5.47872690e-003, 1.49439193e-002, 3.60352909e-002]]])
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:
    model:             Rotated MCA
    n_modes:           10
    power:             1
    max_iter:          1000
    rtol:              1e-08
    squared_loadings:  False
    software:          xeofs
    version:           1.2.0
    date:              2023-10-29 13:29:48

xarray.DataArray

'pvalues_of_left_homogeneous_patterns'

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

• 5.918e-81 4.09e-64 6.234e-58 5.974e-49 ... 0.005479 0.01494 0.03604

  array([[[5.91798633e-081, 4.09037030e-064, 6.23425561e-058, ...,
           4.45216591e-046, 7.99246088e-107, 6.95533368e-234],
          [9.61231974e-109, 4.30909052e-079, 2.41956104e-054, ...,
           1.13316213e-079, 1.17033263e-086, 1.58776474e-105],
          [3.65945612e-169, 1.03636735e-150, 5.85938600e-088, ...,
           2.66091078e-037, 1.39445455e-071, 1.89947470e-183],
          ...,
          [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, ...,
           2.15049395e-115, 1.25690777e-037, 3.79422747e-010],
          [0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ...,
           3.68434950e-006, 4.54442588e-001, 7.47361402e-001],
          [5.82970051e-304, 1.04446785e-286, 1.67982320e-322, ...,
           3.12415723e-031, 2.97910617e-009, 3.83983102e-013],
  ...
          [3.81951811e-033, 1.24766309e-035, 9.83759118e-037, ...,
           1.08208863e-007, 2.09241618e-007, 9.75026228e-007],
          [7.63507975e-018, 1.84998265e-018, 1.85005250e-018, ...,
           1.47665784e-002, 2.23166837e-002, 3.74986857e-002],
          [5.23820545e-013, 5.59326448e-014, 1.46539526e-014, ...,
           7.79688798e-001, 7.97072146e-001, 8.38057312e-001]],
  
         [[1.33009336e-001, 3.88106643e-002, 2.39386496e-003, ...,
           9.33561796e-001, 6.41357062e-001, 5.12874473e-001],
          [2.47971905e-002, 1.49401803e-002, 1.01803200e-003, ...,
           9.48611093e-001, 3.16366009e-001, 6.29374852e-002],
          [1.00284803e-002, 8.50825810e-002, 1.51720312e-001, ...,
           2.62222025e-001, 1.27842313e-002, 1.43419078e-004],
          ...,
          [6.95159651e-010, 1.10849674e-010, 2.79734189e-011, ...,
           2.20338645e-004, 2.58483821e-003, 2.44608586e-002],
          [2.08872724e-002, 1.93770600e-002, 2.08168263e-002, ...,
           3.74564462e-004, 2.56337651e-003, 1.63483691e-002],
          [9.06082861e-001, 8.64987130e-001, 8.39176635e-001, ...,
           5.47872690e-003, 1.49439193e-002, 3.60352909e-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

  n_modes :

  10

  power :

  1

  max_iter :

  1000

  rtol :

  1e-08

  squared_loadings :

  False

  software :

  xeofs

  version :

  1.2.0

  date :

  2023-10-29 13:29:48

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")