<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 EOF analysis - global hgt500
In statistics and signal processing, the method of empirical orthogonal function (EOF) analysis is a decomposition of a signal or data set in terms of orthogonal basis functions which are determined from the data. It is similar to performing a principal components analysis on the data, except that the EOF method finds both temporal projections and spatial patterns. The term is also interchangeable with the geographically weighted PCAs in geophysics (https
Step1:  2. Load ght500 data
Step2:  3. Detrend
Here we used the detrend from scipy.signal. See help
Step3:  4. Carry out EOF analysis
4.1 Create an EOF solver to do the EOF analysis
Cosine of latitude weights are applied before the computation of EOFs.
Step4:  4.3 Retrieve the leading EOFs
Expressed as the correlation between the leading PC time series and the input SST anomalies at each grid point, and the
leading PC time series itself.
Step5:  5. Visualize leading EOFs
5.1 Plot EOFs and PCs
Step6:  5.2 Check variances explained by leading EOFs
<ASSISTANT_TASK:>
Python Code:
% matplotlib inline
import numpy as np
from scipy import signal
import numpy.polynomial.polynomial as poly
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from eofs.standard import Eof
Explanation: EOF analysis - global hgt500
In statistics and signal processing, the method of empirical orthogonal function (EOF) analysis is a decomposition of a signal or data set in terms of orthogonal basis functions which are determined from the data. It is similar to performing a principal components analysis on the data, except that the EOF method finds both temporal projections and spatial patterns. The term is also interchangeable with the geographically weighted PCAs in geophysics (https://en.wikipedia.org/wiki/Empirical_orthogonal_functions). The spatial patterns are the EOFs, and can be thought of as basis functions in terms of variance. The associated temporal projections are the pricipal components (PCs) and are the temporal coefficients of the EOF patterns.
Python EOF library
The EOF library is used here from Dawson (http://ajdawson.github.io/eofs/).
eofs is a Python package for EOF analysis of spatial-temporal data. Using EOFs (empirical orthogonal functions) is a common technique to decompose a signal varying in time and space into a form that is easier to interpret in terms of spatial and temporal variance. Some of the key features of eofs are:
Suitable for large data sets: computationally efficient for the large output data sets of modern climate models.
Transparent handling of missing values: missing values are removed automatically during computations and placed back into output fields.
Automatic metadata: metadata from input fields is used to construct metadata for output fields.
No Compiler required: a fast implementation written in pure Python using the power of numpy, no Fortran or C dependencies.
Data Source
The hgt data is downloaded from https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.pressure.html.
We select hgt at 500mb level from 1979 to 2003, and convert them into yearly climatology using CDO.
cdo yearmean  -sellevel,500/500, -selyear,1979/2003 hgt.mon.mean.nc hgt500.mon.mean.nc
1. Load basic libraries
End of explanation
infile = 'data/hgt500.mon.mean.nc'
ncin = Dataset(infile, 'r')
hgt  = ncin.variables['hgt'][:,0,:,:]
lat  = ncin.variables['lat'][:]
lon  = ncin.variables['lon'][:]
ncin.close()
Explanation: 2. Load ght500 data
End of explanation
nt,nlat,nlon = hgt.shape    
hgt = hgt.reshape((nt,nlat*nlon), order='F')
hgt_detrend = signal.detrend(hgt, axis=0, type='linear', bp=0)
hgt_detrend = hgt_detrend.reshape((nt,nlat,nlon), order='F')
print(hgt_detrend.shape)
Explanation: 3. Detrend
Here we used the detrend from scipy.signal. See help: signal.detrend?
End of explanation
wgts   = np.cos(np.deg2rad(lat))
wgts   = wgts.reshape(len(wgts), 1)
solver = Eof(hgt_detrend, weights=wgts)
Explanation: 4. Carry out EOF analysis
4.1 Create an EOF solver to do the EOF analysis
Cosine of latitude weights are applied before the computation of EOFs.
End of explanation
eof1 = solver.eofs(neofs=10)
pc1  = solver.pcs(npcs=10, pcscaling=0)
varfrac = solver.varianceFraction()
lambdas = solver.eigenvalues()
Explanation: 4.3 Retrieve the leading EOFs
Expressed as the correlation between the leading PC time series and the input SST anomalies at each grid point, and the
leading PC time series itself.
End of explanation
parallels = np.arange(-90,90,30.)
meridians = np.arange(-180,180,30)
for i in range(0,2):
    fig = plt.figure(figsize=(12,9))
    plt.subplot(211)
    
    m = Basemap(projection='cyl', llcrnrlon=min(lon), llcrnrlat=min(lat), urcrnrlon=max(lon), urcrnrlat=max(lat))    
    x, y = m(*np.meshgrid(lon, lat))
    clevs = np.linspace(np.min(eof1[i,:,:].squeeze()), np.max(eof1[i,:,:].squeeze()), 21)
    cs = m.contourf(x, y, eof1[i,:,:].squeeze(), clevs, cmap=plt.cm.RdBu_r)
    m.drawcoastlines()  
    m.drawparallels(parallels, labels=[1,0,0,0])
    m.drawmeridians(meridians, labels=[1,0,0,1])
    cb = m.colorbar(cs, 'right', size='5%', pad='2%')
    cb.set_label('EOF', fontsize=12)
    plt.title('EOF ' + str(i+1), fontsize=16)
    plt.subplot(212)
    days = np.linspace(1979,2003,nt)
    plt.plot(days, pc1[:,i], linewidth=2)
    plt.axhline(0, color='k')
    plt.xlabel('Year')
    plt.ylabel('PC Amplitude')   
    plt.ylim(np.min(pc1.squeeze()), np.max(pc1.squeeze()))    
Explanation: 5. Visualize leading EOFs
5.1 Plot EOFs and PCs
End of explanation
plt.figure(figsize=(11,6))
eof_num = range(1, 16)
plt.plot(eof_num, varfrac[0:15], linewidth=2)
plt.plot(eof_num, varfrac[0:15], linestyle='None', marker="o", color='r', markersize=8)
plt.axhline(0, color='k')
plt.xticks(range(1, 16))
plt.title('Fraction of the total variance represented by each EOF')
plt.xlabel('EOF #')
plt.ylabel('Variance Fraction')
plt.xlim(1, 15)
plt.ylim(np.min(varfrac), np.max(varfrac)+0.01)
Explanation: 5.2 Check variances explained by leading EOFs
End of explanation
<END_TASK>