ノーマルなEOFは定在的な変動しか表現できません。空間的に伝播する現象が卓越している場合には複素EOF(complex EOF)解析が適しています。

n行(時間数)m列(格子点数)のデータ行列dataを作ります。平均値を引いておくとかの前処理は済んでいるとします。

ヒルベルト変換を行うためノーマルなEOFと違ってデータの有限性に起因する歪みを軽減するためにまずテーパリングを行います。中央銀行による金融資産買い入れの段階的縮小のことではありません。
以下では例としてコサインテーパーをかけてみます。
from scipy.signal import hilbert, tukey
import copy

wdw = tukey( n, alpha=0.2 ) # 両側合わせて10*(alpha) % のコサインテーパー
data_w = copy.deepcopy( data )
for i in range( m ):
        data_w[:, i] = data_w[:, i] * wdw
tukeyはコサインテーパーウィンドウを返します。alpha=0.2とすると両側合わせて20%のcos20テーパーになります。もちろん別の窓関数でも構いません。Scipyには様々な種類の窓関数があります。

ヒルベルト変換を行い複素数の共分散を求めて固有値問題を解きます。
data_w = hilbert( data_w, axis=0 )
bunsan = np.dot( np.conj( data_w.T ), data_w ) / n
values, vectors = np.linalg.eigh( bunsan ) 
vectors = vectors[:, ::-1]
values = values[::-1]
固有値valuesは実数、固有ベクトルvectorsは複素数です。
固有値が小さいほうから順に並んでいるので、大きいほうから並ぶように反転させています。
各モードについて、寄与率を計算して、半周期ぶんの固有ベクトル(空間関数)の変化を描画します。
for mode in range(0, 3):
    kiyo = 100 * values[mode] / np.sum(values) # 寄与率
    print('第', (mode+1), 'モード 寄与率= ', kiyo, ' %')

    score = np.dot( data, vectors[:, mode] ) # 主成分=時係数(複素数)
    amp_score = abs( score * np.conj( score ) ) # 時係数の振幅
    phs_score = np.angle( score ) * 180 / np.pi # 時係数の位相

    # 位相が0度から180度まで45度毎に固有ベクトルを求める
    agl = np.arange(0, 181, 45) * np.pi / 180.
    ee = np.cos(agl) + np.sin(agl) * 1j # 大きさが1で偏角がaglの複素数

    for k in np.argsort(agl):
        vec = ee[k] * np.conj( vectors[:, mode] ) # 位相がagl[k]の時の固有ベクトル
( vec.realをマッピングする )
時係数の振幅と位相を一つの図の中に描いてみます。
    fig = plt.figure()
    ax1 = fig.subplots(1, 1)
    ax2 = ax1.twinx() # x軸を共有して右側にも軸を作る

    x = np.arange(0, n)
    ax1.plot(x, amp_score/np.max(amp_score), color='k', ls='-', lw = 3.)
    ax1.set_ylabel('Normalized amplitude of temporal function', fontsize=13)
    ax1.set_yticks(np.arange(0, 1.1, 0.25))
    ax1.set_ylim([0, 1])
    ax1.tick_params(axis='both', labelsize=13)
    ax1.grid(which='major', color='gray', linestyle='dotted')
    ax1.set_xlabel('Time', fontsize=14)

    ax2.plot(x, phs_score, color='b', ls='--', lw = 1.5)
    ax2.set_ylabel('Phase [$\u00B0$]', fontsize=14, color='b')
    ax2.set_ylim([-180, 180])
    ax2.set_yticks(np.arange(-180, 181, 45))
    ax2.tick_params(axis='y', colors='b', labelsize=13)

    plt.tight_layout()
    plt.show()