docs: add MCA example

docs/auto_examples/2cross/plot_mca.ipynb

@@ -0,0 +1,144 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n# Maximum Covariance Analysis\n\nMaximum Covariance Analysis (MCA) between two data sets.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Load packages and data:\nimport numpy as np\nimport xarray as xr\nimport matplotlib.pyplot as plt\nfrom matplotlib.gridspec import GridSpec\nfrom cartopy.crs import Orthographic, PlateCarree\nfrom cartopy.feature import LAND\n\nfrom xeofs.xarray import MCA"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Create 2 different DataArrays\n\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "t2m = xr.tutorial.load_dataset('air_temperature')['air']\nda1 = t2m.isel(lon=slice(0, 26))\nda2 = t2m.isel(lon=slice(27, None))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Perform MCA\n\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "mca = MCA(\n X=da1, Y=da2,\n n_modes=20,\n dim='time',\n norm=False,\n weights_X='coslat',\n weights_Y='coslat'\n)\nmca.solve()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Get singular vectors, projections (PCs), homogeneous and heterogeneous\npatterns:\n\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "singular_vectors = mca.singular_vectors()\npcs = mca.pcs()\nhom_pats, pvals_hom = mca.homogeneous_patterns()\nhet_pats, pvals_het = mca.heterogeneous_patterns()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Create a mask to identifiy where p-values are below 0.05\n\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "hom_mask = [values < 0.05 for values in pvals_hom]\nhet_mask = [values < 0.05 for values in pvals_het]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Plot some relevant quantities of mode 2.\n\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "lonlats = [\n np.meshgrid(pvals_hom[0].lon.values, pvals_hom[0].lat.values),\n np.meshgrid(pvals_hom[1].lon.values, pvals_hom[1].lat.values)\n]\nproj = [\n Orthographic(central_latitude=30, central_longitude=-120),\n Orthographic(central_latitude=30, central_longitude=-60)\n]\nkwargs1 = {\n 'cmap' : 'BrBG', 'vmin' : -.05, 'vmax': .05, 'transform': PlateCarree()\n}\nkwargs2 = {\n 'cmap' : 'RdBu', 'vmin' : -1, 'vmax': 1, 'transform': PlateCarree()\n}\n\nmode = 2\n\nfig = plt.figure(figsize=(7, 14))\ngs = GridSpec(5, 2)\nax1 = [fig.add_subplot(gs[0, i], projection=proj[i]) for i in range(2)]\nax2 = [fig.add_subplot(gs[1, i], projection=proj[i]) for i in range(2)]\nax3 = [fig.add_subplot(gs[2, i], projection=proj[i]) for i in range(2)]\nax4 = [fig.add_subplot(gs[3, i]) for i in range(2)]\n\nfor i, a in enumerate(ax1):\n singular_vectors[i].sel(mode=mode).plot(ax=a, **kwargs1)\n\nfor i, a in enumerate(ax2):\n hom_pats[i].sel(mode=mode).plot(ax=a, **kwargs2)\n a.scatter(\n lonlats[i][0], lonlats[i][1], hom_mask[i].sel(mode=mode).values * .5,\n color='k', alpha=.5, transform=PlateCarree()\n )\nfor i, a in enumerate(ax3):\n het_pats[i].sel(mode=mode).plot(ax=a, **kwargs2)\n a.scatter(\n lonlats[i][0], lonlats[i][1], het_mask[i].sel(mode=mode).values * .5,\n color='k', alpha=.5, transform=PlateCarree()\n )\n\nfor i, a in enumerate(ax4):\n pcs[i].sel(mode=mode).plot(ax=a)\n a.set_xlabel('')\n\n\nfor a in np.ravel([ax1, ax2, ax3]):\n a.coastlines(color='.5')\n a.add_feature(LAND)\n\nplt.tight_layout()\nplt.savefig('mca.jpg')"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.8.13"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}

docs/auto_examples/2cross/plot_mca.py

@@ -0,0 +1,108 @@
+"""
+Maximum Covariance Analysis
+===========================
+
+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.xarray import MCA
+
+#%%
+# 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(
+ X=da1, Y=da2,
+ n_modes=20,
+ dim='time',
+ norm=False,
+ weights_X='coslat',
+ weights_Y='coslat'
+)
+mca.solve()
+
+#%%
+# Get singular vectors, projections (PCs), homogeneous and heterogeneous
+# patterns:
+
+singular_vectors = mca.singular_vectors()
+pcs = mca.pcs()
+hom_pats, pvals_hom = mca.homogeneous_patterns()
+het_pats, pvals_het = mca.heterogeneous_patterns()
+
+#%%
+# 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' : -.05, 'vmax': .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 * .5,
+ color='k', alpha=.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 * .5,
+ color='k', alpha=.5, transform=PlateCarree()
+ )
+
+for i, a in enumerate(ax4):
+ pcs[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('mca.jpg')

docs/auto_examples/2cross/plot_mca.py.md5

@@ -0,0 +1 @@
+91ab56bff33c723e4b800b87cc9251c7