Step by step Tutorial¶

In this tutorial we will go throught the entire implementation of spy4cast MCA methodology

On the first section we will do everything using xarray and on the second one we will introduce the abstractions that spy4cast implements to enhace productivity

Using Xarray

Configuration
Preprocessing
MCA
Plot results

Using Spy4Cast

Configuration
Preprocessing
MCA
Plot results
Crossvalidation

Using xarray¶

In this case of study we will go through the relation between sea surface temperature in the Atlantic ocean and the clorophyl concentration in the shouth coast of the Iberian Peninsula. It is the code used to replicate the results from Lopez-Parages, Jorge, et al. 2022.

[1]:

import xarray as xr
import os
import numpy as np
import scipy
from scipy import signal, sparse, stats
import matplotlib.pyplot as plt
from matplotlib import gridspec
import cartopy.crs as ccrs

/Users/pabloduran/opt/anaconda3/envs/Spy4Cast/lib/python3.9/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

3. Perform MCA¶

[22]:

nm = 3  # Calculate up to 3 modes
alpha = 0.1  # Significance level for Pearson correlation

Calculate the convariance matrix

[23]:

c = np.dot(y[~y_land_mask], np.transpose(z[~z_land_mask]))

Singular Value Decomposition

[24]:

r, _d, q = sparse.linalg.svds(c, k=nm, which='LM')  # Which LM = Large magnitude

# Modes are reversed so we reverse them in r, d and q
r = r[:, ::-1]
_d = _d[::-1]
q = q[::-1, :]

svdvals = scipy.linalg.svdvals(c)
scf = svdvals[:nm] / np.sum(svdvals)

[25]:

# Time series of covariability
u = np.dot(np.transpose(y[~y_land_mask]), r)
v = np.dot(np.transpose(z[~z_land_mask]), q.transpose())

Us = ((u - u.mean(0)) / u.std(0)).transpose()  # Standarized anom across axis 0
Vs = ((v - v.mean(0)) / v.std(0)).transpose()

Perform index regression to get the modes of variability

[26]:

ny, nz = y.shape[0], z.shape[0]
nt = y.shape[1]  # Must be the same as z.shape[1]

RUY = np.zeros([ny, nm], dtype=np.float32)
SUY = np.zeros([ny, nm], dtype=np.float32)
pvalruy = np.zeros([ny, nm], dtype=np.float32)

RUZ = np.zeros([nz, nm], dtype=np.float32)
SUZ = np.zeros([nz, nm], dtype=np.float32)
pvalruz = np.zeros([nz, nm], dtype=np.float32)

RUY[y_land_mask, :] = np.nan
pvalruy[y_land_mask, :] = np.nan

RUZ[z_land_mask, :] = np.nan
pvalruz[z_land_mask, :] = np.nan

for i in range(nm):
    SUY[:, i] = y.dot(Us[i, :]) / nt
    result = np.apply_along_axis(stats.pearsonr, 1, y[~y_land_mask], Us[i, :])
    RUY[~y_land_mask, i], pvalruy[~y_land_mask, i] = result[:, 0], result[:, 1]

    SUZ[:, i] = z.dot(Us[i, :]) / nt
    result = np.apply_along_axis(stats.pearsonr, 1, z[~z_land_mask], Us[i, :])
    RUZ[~z_land_mask, i], pvalruz[~z_land_mask, i] = result[:, 0], result[:, 1]

RUY_sig = RUY.copy()
SUY_sig = SUY.copy()
RUY_sig[pvalruy > alpha] = np.nan
SUY_sig[pvalruy > alpha] = np.nan

RUZ_sig = RUZ.copy()
SUZ_sig = SUZ.copy()
RUZ_sig[pvalruz > alpha] = np.nan
SUZ_sig[pvalruz > alpha] = np.nan

4. Plot MCA¶

[27]:

fig = plt.figure(figsize=(14, 8))

<Figure size 1400x800 with 0 Axes>

[28]:

gs = gridspec.GridSpec(3, 3)

[29]:

for i in range(nm):
    ax = fig.add_subplot(gs[i, 0])
    ax.plot(sst["year"].values, Us[i, :], color="green", label="Us")
    ax.plot(sst["year"].values, Vs[i, :], color="blue", label="Vs")
    if i == 0:
        ax.legend()

for i in range(nm):
    ax = fig.add_subplot(gs[i, 1], projection=ccrs.PlateCarree())

    ru, ru_sig, lats, lons = RUY, RUY_sig, sst["latitude"], sst["longitude"]
    t = ru[:, i].transpose().reshape((len(lats), len(lons)))
    th = ru_sig[:, i].transpose().reshape((len(lats), len(lons)))

    im = ax.contourf(
        lons, lats, t, cmap="bwr",
        extend='both', transform=ccrs.PlateCarree(),
    )
    ax.contourf(
        lons, lats, th, colors='none', hatches='..', extend='both',
        transform=ccrs.PlateCarree()
    )
    ax.coastlines()
    fig.colorbar(im, ax=ax, orientation="horizontal")

for i in range(nm):
    ax = fig.add_subplot(gs[i, 2], projection=ccrs.PlateCarree())

    ru, ru_sig, lats, lons = RUZ, RUZ_sig, chl["lat"], chl["lon"]
    t = ru[:, i].transpose().reshape((len(lats), len(lons)))
    th = ru_sig[:, i].transpose().reshape((len(lats), len(lons)))

    im = ax.contourf(
        lons, lats, t, cmap="viridis",
        extend='both', transform=ccrs.PlateCarree(),
    )
    ax.contourf(
        lons, lats, th, colors='none', hatches='..', extend='both',
        transform=ccrs.PlateCarree()
    )
    ax.coastlines()
    fig.colorbar(im, ax=ax, orientation="horizontal")

fig.suptitle("MCA")

[29]:

Text(0.5, 0.98, 'MCA')

[30]:

fig

[30]:

Using Spy4Cast¶

[31]:

from spy4cast import Dataset, Region, Month
from spy4cast.spy4cast import Preprocess, MCA, Crossvalidation

1. Configuration¶

Configuration is srtraight forward with the abstractions spy4cast provides.

This enhances productivity and makes it easier to change the parameters of the run. It is also quite straight-forward to declare seasons that span through a year like NDJ, which requires more logic with xarray

Also it provides a skip option to reduce the dimension of the dataset. This makes ir faster for the first runs and you can decrease it to 0 for the final output

[55]:

predictor = Dataset(PREDICTOR_NAME, DATASET_DIR).open("sst")
sst_region = Region(
    lat0=5, latf=25,
    lon0=-75, lonf=-20,
    month0=Month.AUG, monthf=Month.SEP,
    year0=1997, yearf=2019,
)  # PREDICTOR: Y
predictor.slice(sst_region, skip=0)

predictand = Dataset(PREDICTAND_NAME, DATASET_DIR).open("CHL")
chl_region = Region(
    lat0=36, latf=37,
    lon0=-5.3, lonf=-2,
    month0=Month.MAR, monthf=Month.APR,
    year0=1998, yearf=2020,
)  # PRECITAND: Z
predictand.slice(chl_region, skip=0)

[55]:

<spy4cast.dataset.Dataset at 0x7f9ec9321640>

Climatology and anomaly plots can be created easily with meteo.Clim and meteo.Anom as well as time series

[56]:

from spy4cast.meteo import Clim

Clim(predictand, "map").plot(levels=[0, .2, .4, .5, .6, .8, 1, 2], cmap="viridis", figsize=(10, 4))

[56]:

((<Figure size 1000x400 with 2 Axes>,), (<GeoAxesSubplot:>,))

2. Preprocessing¶

Preprocessing only takes one line, as everything we did before is included.

It also performs a Butterwirth filter if you provide the period and order keyword arguments. To do so it uses scipy.signal.butter along with scipy.signal.filtfilt functions along the time axis

[59]:

predictor_preprocessed = Preprocess(predictor)
predictand_preprocessed = Preprocess(predictand)

[INFO] Preprocessing data for variable sst took: 0.031 seconds
[INFO] Preprocessing data for variable CHL took: 0.017 seconds

To visalize easily the region of both variables we can use .plot method.

[60]:

predictor_preprocessed.plot(selected_year=1999, figsize=(8, 4), levels=20)
predictand_preprocessed.plot(selected_year=2000, figsize=(8, 4), cmap="viridis", levels=[0, .1, .2, .4, .5, .6, .7, .8, .9, 1, 2])

[60]:

((<Figure size 800x400 with 2 Axes>,),
 (<GeoAxesSubplot:title={'center':'Year 2000'}>,))

To save the calculated matrices, spy4cast includes a convenient method in case there is an intention to share the results

[61]:

DATA_FOLDER = 'step-by-step-tutorial/data/'
PREDICTOR_PREPROCESSED_PREFIX = 'predictor_'
PREDICTAND_PREPROCESSED_PREFIX = 'predictand_'

predictor_preprocessed.save(PREDICTOR_PREPROCESSED_PREFIX, DATA_FOLDER)
predictand_preprocessed.save(PREDICTAND_PREPROCESSED_PREFIX, DATA_FOLDER)

[INFO] Saving Preprocess data in `step-by-step-tutorial/data/predictor_*.npy`
[INFO] Saving Preprocess data in `step-by-step-tutorial/data/predictand_*.npy`

They will be saved in .npy format. Eventhough this format is not the best it promotes compatibility and it is flexible so that it can be used with any other tools. In the future, if there is a demand from the users, it will be easy to move to any format it is suggested. To load the different matrices that are saved only the np.load method is needed.

However, spy4cast has a more convenient function to get the same object you have saved from the previous call

This is useful to write code in the following manner:

# Once these have been run we set to False and the program is faster
RUN_PREPROCESS_Y = True
RUN_PREPROCESS_Z = True

if RUN_PREPROCESS_Y:
    y = Preprocess(sst)
    y.save("y_", "data")
else:
    y = Preprocess.load("y_", "data")

if RUN_PREPROCESS_Z:
    z = Preprocess(chl)
    z.save("z_", "data")
else:
    z = Preprocess.load("z_", "data")

mca = MCA(y, z, 3, .1)

### 3. MCA

MCA is built into spy4cast. It performs the exact same methodology as before and returns a variable contaning all the matrices as before

This is very convenient for user experience as all the date is stored together

[62]:

nm = 3
alpha = .1

mca = MCA(predictor_preprocessed, predictand_preprocessed, nm, alpha)

[INFO] Applying MCA
    Shapes: Z(24750, 23)
            Y(1100, 23)
    Regions: Z MA (36.00ºN, 37.00ºN - 5.30ºW, 2.00ºW)
            Y AS (5.00ºN, 25.00ºN - 75.00ºW, 20.00ºW)
       Took: 4.556 seconds

[63]:

mca.var_names

[63]:

('psi',
 'RUY',
 'RUY_sig',
 'SUY',
 'SUY_sig',
 'RUZ',
 'RUZ_sig',
 'SUZ',
 'SUZ_sig',
 'pvalruz',
 'pvalruy',
 'Us',
 'Vs',
 'scf',
 'alpha',
 'nm')

As in Preprocess, there is a method to save the data into .npy format. This is useful when the user needs to change the plots and wants to use the results without running the methodology again

As in Preprocess, it can be used the following way:

RUN_PREPROCESS_MCA = True

if RUN_PREPROCESS_Y:
    mca = MCA(y, z, nm, alpha)
    mca.save("mca_", "data")
else:
    mca = MCA.load("mca_", "data", dsy=y, dsz=z)

[64]:

MCA_PREFIX = 'mca_'

mca.save(MCA_PREFIX, DATA_FOLDER)
# mca = MCA.load(MCA_PREFIX, DATA_FOLDER, dsy=predictor_preprocessed, dsz=predictand_preprocessed)

[INFO] Saving MCA data in `step-by-step-tutorial/data/mca_*.npy`

4. Plot MCA¶

Plotting takes one command also and it is very easy to save the plots and to change the visulization

These plots must be used as a quick-diagnosis tool and not as a final result. We encourage users to write their own plot functions.

[65]:

PLOTS_FOLDER = 'step-by-step-tutorial/plots'

MCA_PLOT_NAME = 'mca.png'

mca.plot(save_fig=True, cmap='viridis', name=MCA_PLOT_NAME,
         figsize=(14, 8),
         width_ratios=[1, 1, 1],
         folder=PLOTS_FOLDER,)

[INFO] Saving plot with path step-by-step-tutorial/plots/mca.png

[65]:

((<Figure size 1400x800 with 11 Axes>,),
 (<AxesSubplot:title={'center':'Us Vs mode 1'}>,
  <AxesSubplot:title={'center':'Us Vs mode 2'}>,
  <AxesSubplot:title={'center':'Us Vs mode 3'}>,
  <GeoAxesSubplot:title={'center':'RUY mode 1. SCF=47.1%'}>,
  <GeoAxesSubplot:title={'center':'RUY mode 2. SCF=12.0%'}>,
  <GeoAxesSubplot:title={'center':'RUY mode 3. SCF=9.3%'}>,
  <GeoAxesSubplot:title={'center':'RUZ mode 1. SCF=47.1%'}>,
  <GeoAxesSubplot:title={'center':'RUZ mode 2. SCF=12.0%'}>,
  <GeoAxesSubplot:title={'center':'RUZ mode 3. SCF=9.3%'}>))

### 5. Crossvalidation

Crossvalidation is the same as MCA but leaving one year out of the analysis. The code for xarray-MCA can be modified to run crossvalidation. Here we show how it is run with spy4cast. It has the same .save, .load and .plot functions as MCA

[66]:

cross = Crossvalidation(predictor_preprocessed, predictand_preprocessed, nm, alpha)

[INFO] Applying Crossvalidation
    Shapes: Z(24750, 23)
            Y(1100, 23)
    Regions: Z MA (36.00ºN, 37.00ºN - 5.30ºW, 2.00ºW)
            Y AS (5.00ºN, 25.00ºN - 75.00ºW, 20.00ºW)
        year: 1 of 23
        year: 2 of 23
        year: 3 of 23
        year: 4 of 23
        year: 5 of 23
        year: 6 of 23
        year: 7 of 23
        year: 8 of 23
        year: 9 of 23
        year: 10 of 23
        year: 11 of 23
        year: 12 of 23
        year: 13 of 23
        year: 14 of 23
        year: 15 of 23
        year: 16 of 23
        year: 17 of 23
        year: 18 of 23
        year: 19 of 23
        year: 20 of 23
        year: 21 of 23
        year: 22 of 23
        year: 23 of 23

        Took: 150.828 seconds

[67]:

CROSS_PREFIX = 'cross_'

cross.save(CROSS_PREFIX, DATA_FOLDER)

[INFO] Saving Crossvalidation data in `step-by-step-tutorial/data/cross_*.npy`

[68]:

# cross = Crossvalidation.load(CROSS_PREFIX, DATA_FOLDER, dsy=predictor_preprocessed, dsz=predictand_preprocessed)

[69]:

CROSS_PLOT_NAME = 'cross.png'

cross.plot(save_fig=True, folder=PLOTS_FOLDER, name=CROSS_PLOT_NAME, map_levels=[-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6], rmse_levels=[0, .2, .4, .5, .6, .8, 1, 2])

[INFO] Saving plot with path step-by-step-tutorial/plots/cross.png

[69]:

((<Figure size 1600x800 with 6 Axes>,),
 (<GeoAxesSubplot:title={'center':'ACC map'}>,
  <AxesSubplot:title={'center':'ACC time series'}>,
  <GeoAxesSubplot:title={'center':'RMSE map'}>,
  <AxesSubplot:title={'center':'RMSE time series'}>))

Crossvalidation provides a different version of diagnosis plot with version=2 by which you also need to provide mca=mca

[70]:

cross.plot(save_fig=True, folder=PLOTS_FOLDER, name=CROSS_PLOT_NAME, version=2, mca=mca)

[INFO] Saving plot with path step-by-step-tutorial/plots/cross.png

[70]:

((<Figure size 2400x850 with 16 Axes>,),
 (<AxesSubplot:title={'center':'Mode 1, |ruv|=0.50 p=0.985'}, xlabel='Years', ylabel='nº std'>,
  <AxesSubplot:title={'center':'Mode 1, |ruv|=0.50 , scf=47.10%'}, xlabel='Years', ylabel='|ruv|'>,
  <GeoAxesSubplot:title={'center':'ACC 1 modes'}>,
  <AxesSubplot:title={'center':'Skill (bars), \n MSESS (lines)'}>,
  <AxesSubplot:title={'center':'Mode 2, |ruv|=0.70 p=1.000'}, xlabel='Years'>,
  <AxesSubplot:title={'center':'Mode 2, |ruv|=0.70 , scf=11.96%'}, xlabel='Years'>,
  <GeoAxesSubplot:title={'center':'ACC 2 modes'}>,
  <GeoAxesSubplot:title={'center':'ACC [1 2]) - [1]'}>,
  <AxesSubplot:title={'center':'Mode 3, |ruv|=0.49 p=0.982'}, xlabel='Years'>,
  <AxesSubplot:title={'center':'Mode 3, |ruv|=0.49 , scf=9.26%'}, xlabel='Years'>,
  <GeoAxesSubplot:title={'center':'ACC 3 modes'}>,
  <GeoAxesSubplot:title={'center':'ACC [1 2 3]) - [1 2]'}>))