<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Identify Coupled Patterns between SLP and SST through Maximum Covariance Analysis
Maximum Correlation Analysis (MCA; Bretherton et al., 1992) is similar to Empirical Orthogonal Function Analysis (EOF) in that they both deal with the decomposition of a covariance matrix. In EOF, this is a covariance matrix based on a single spatio-temporal field, while MCA is based on the decomposition of a "cross-covariance" matrix derived from two variables. The resulting expansion coefficients (ECs) associated with the left and right singular vectors can then be projected onto the original data to obtain homogeneous and heterogeneous regression maps. 
This notebook applies MCA method to real data—the tropical Pacific sea level pressure (SLP) and sea surface temperature (SST) fields—and identify the coupled patterns from the data, including the famous tropical Pacific climate variability known as the El Nin o–Southern Oscillation (ENSO). ENSO has warm El Nino states and cool La Nina states, with changes found not only in the SST but also in the SLP.
The MCA approach is chosen because it is able to capture patterns of maximum covariance between two variables; it has been found to reasonably capture atmospheric and oceanic processes (Wilks, 2015). It is a robust
method to investigate dominant modes of interaction, because it favors a better understanding of the relationship between groups of variables (Frankignoul et al., 2011). 
The SLP and SST data are downloaded from http
Step1:  2. Load data
Spatial domain
Step2:  2.2 SST
Step3:  2.3 Preprocess
Convert 3D SLP and SST to 2D arrays and only use none-NaN values of SST.
Step4:  3. Carry out Maximum Covariance Analysis
The key step is to apply SVD to decompose the covariance matrix of SLP and SST. It is worth noting that np.linalg.svd will produce V.T, instead of V. In practice, have to transpose it.
3.1 MCA
Step5:  3.2 Postprocess
Calculate cumulative fraction of squares covariance explained
Extract the leading SLP MCA pattern and EC
Extract the leading SST MCA pattern and EC
Normalize MCD mode 1 by standardizing the EC1, so patterns correspond to a 1-std deviation variation in EC1. The positive EC1 of SST should match the Nino SSTA.
3.2.1 Calculate cumulative fraction of squares covariance explained
Step6:  3.2.2 Extract the leading SLP MCA pattern and EC
Step7:  3.2.3 Extract the leading SST MCA pattern and EC
Step8:  4 Visualize MCA results
4.1 Plot cumulative fraction of squares covariance explained
Step9:  4.2 Plot the leading SLP/SST MCA spatial pattern and EC
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import datetime as dt
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.dates as mdates
mpl.rcParams['figure.figsize'] = 8.0, 4.0 
mpl.rcParams['font.size'] = 13  
Explanation: Identify Coupled Patterns between SLP and SST through Maximum Covariance Analysis
Maximum Correlation Analysis (MCA; Bretherton et al., 1992) is similar to Empirical Orthogonal Function Analysis (EOF) in that they both deal with the decomposition of a covariance matrix. In EOF, this is a covariance matrix based on a single spatio-temporal field, while MCA is based on the decomposition of a "cross-covariance" matrix derived from two variables. The resulting expansion coefficients (ECs) associated with the left and right singular vectors can then be projected onto the original data to obtain homogeneous and heterogeneous regression maps. 
This notebook applies MCA method to real data—the tropical Pacific sea level pressure (SLP) and sea surface temperature (SST) fields—and identify the coupled patterns from the data, including the famous tropical Pacific climate variability known as the El Nin o–Southern Oscillation (ENSO). ENSO has warm El Nino states and cool La Nina states, with changes found not only in the SST but also in the SLP.
The MCA approach is chosen because it is able to capture patterns of maximum covariance between two variables; it has been found to reasonably capture atmospheric and oceanic processes (Wilks, 2015). It is a robust
method to investigate dominant modes of interaction, because it favors a better understanding of the relationship between groups of variables (Frankignoul et al., 2011). 
The SLP and SST data are downloaded from http://www.esrl.noaa.gov/psd/gcos_wgsp/Gridded/data.hadslp2.html and http://www.esrl.noaa.gov/psd/thredds/catalog/Datasets/kaplan_sst/catalog.html, respectively.
1. Load all needed libraries
End of explanation
ds1 = xr.open_dataset('data/slp.mnmean.hadslp2.nc')
slp = ds1.slp.sel(lat=slice(30, -30), lon=slice(180, 290), time=slice('1950-01-01','2005-12-31'))
lon_slp = ds1.lon.sel(lon=slice(180, 290))
lat_slp = ds1.lat.sel(lat=slice(30, -30))
dates = ds1.time.sel(time=slice('1950-01-01','2005-12-31')).values
# climatology
slp_clm = slp.groupby('time.month').mean(dim='time')
# anomaly
slp_anom = slp.groupby('time.month') - slp_clm
Explanation: 2. Load data
Spatial domain: -180W ~ -70W, and 30N ~ -30S
Period: 1950~2005
Raw SLP and SST are converted to monthly anomalies
2.1 SLP
End of explanation
ds2 = xr.open_dataset('data/sst.mon.anom.kaplan.nc')
sst_anom = ds2.sst.sel(lat=slice(-30, 30), lon=slice(180, 290), time=slice('1950-01-01','2005-12-31'))
lat_sst = ds2.lat.sel(lat=slice(-30, 30))
lon_sst = ds2.lon.sel(lon=slice(180, 290))
Explanation: 2.2 SST
End of explanation
slp2d = slp_anom.values
ntime, nrow_slp, ncol_slp = slp2d.shape
slp2d = np.reshape(slp2d, (ntime, nrow_slp*ncol_slp), order='F')
sst2d = sst_anom.values
ntime, nrow_sst, ncol_sst = sst2d.shape
sst2d = np.reshape(sst2d, (ntime, nrow_sst*ncol_sst), order='F')
nonMissingIndex = np.where(np.isnan(sst2d[0]) == False)[0]
sst2dNoMissing = sst2d[:, nonMissingIndex]
Explanation: 2.3 Preprocess
Convert 3D SLP and SST to 2D arrays and only use none-NaN values of SST.
End of explanation
Cxy = np.dot(slp2d.T, sst2dNoMissing)/(ntime-1.0)
U, s, V = np.linalg.svd(Cxy, full_matrices=False)
V = V.T
Explanation: 3. Carry out Maximum Covariance Analysis
The key step is to apply SVD to decompose the covariance matrix of SLP and SST. It is worth noting that np.linalg.svd will produce V.T, instead of V. In practice, have to transpose it.
3.1 MCA
End of explanation
scf = s**2./np.sum(s**2.0)
Explanation: 3.2 Postprocess
Calculate cumulative fraction of squares covariance explained
Extract the leading SLP MCA pattern and EC
Extract the leading SST MCA pattern and EC
Normalize MCD mode 1 by standardizing the EC1, so patterns correspond to a 1-std deviation variation in EC1. The positive EC1 of SST should match the Nino SSTA.
3.2.1 Calculate cumulative fraction of squares covariance explained
End of explanation
# SLP MCA pattern
U1 = np.reshape(U[:,0, None], (nrow_slp, ncol_slp), order='F')
# EC1 of SLP
a1 = np.dot(slp2d, U[:,0, np.newaxis])
# normalize
U1_norm = U1*np.std(a1)
a1_norm = a1/np.std(a1)
Explanation: 3.2.2 Extract the leading SLP MCA pattern and EC
End of explanation
# SST MCA pattern
V1 = np.ones([nrow_sst*ncol_sst,1]) * np.NaN
V1 = V1.astype(V.dtype)
V1[nonMissingIndex,0] = V[:,0]
V1 = V1.reshape([nrow_sst,ncol_sst], order='F')
# EC1 of SST
b1 = np.dot(sst2dNoMissing, V[:,0, np.newaxis])
# normalize
V1_norm = V1*np.std(b1)
b1_norm = b1/np.std(b1)
Explanation: 3.2.3 Extract the leading SST MCA pattern and EC
End of explanation
plt.plot(np.cumsum(scf),'x') 
plt.xlabel('SVD mode')
plt.ylabel('Cumulative squares covariance fraction')
plt.ylim([0.7,1.1])
plt.xlim([-0.5, 40])
Explanation: 4 Visualize MCA results
4.1 Plot cumulative fraction of squares covariance explained
End of explanation
gs = gridspec.GridSpec(2, 2)
gs.update(wspace=0.1, hspace=0.15)
fig = plt.figure(figsize = (14,10))
levels = np.arange(-1.0, 1.01, 0.05)
# SLP Pattern
ax0 = fig.add_subplot(gs[0,0], projection=ccrs.PlateCarree())
x1, y1 = np.meshgrid(lon_slp, lat_slp)
cs = ax0.contourf(x1, y1, U1_norm,
                        levels=levels,
                        transform=ccrs.PlateCarree(),
                        cmap='RdBu_r')
    
cb=fig.colorbar(cs, ax=ax0, shrink=0.8, aspect=20)    
ax0.coastlines()
ax0.set_global()
ax0.set_extent([-180, -70, -19, 19])
ax0.set_title('Normalized SLP MCA Mode 1')
  
    
# SST Pattern
ax1 = fig.add_subplot(gs[0,1], projection=ccrs.PlateCarree())
x2, y2 = np.meshgrid(lon_sst, lat_sst)
cs2 = ax1.contourf(x2, y2, V1_norm,
                        levels=levels, 
                        transform=ccrs.PlateCarree(),
                        cmap='RdBu_r')
    
cb=fig.colorbar(cs, ax=ax1, shrink=0.8, aspect=20)    
ax1.coastlines()
ax1.set_global()
ax1.set_extent([-180, -70, -19, 19])
ax1.set_title('Normalized SST MCA Mode 1')
# EC1
ax2 = fig.add_subplot(gs[1,:]) 
ax2.plot(dates, a1_norm, label='SLP')
ax2.plot(dates, b1_norm, label='SST')
r = np.corrcoef(a1[:,0], b1[:,0])[0, 1]
ax2.set_title('Expansion Coefficients: SFC = '+ str(round(scf[0],2)) + ', R = ' + str(round(r,2)))
ax2.legend()
ax2.set_ylim([-4,4])
ax2.format_xdata = mdates.DateFormatter('%Y')
fig.autofmt_xdate()
Explanation: 4.2 Plot the leading SLP/SST MCA spatial pattern and EC
End of explanation
<END_TASK>