Nino SST Indices

Niño 3.4 (5N-5S, 170W-120W): Niño 3.4异常可以被认为代表从日变线左右到南美海岸横跨太平洋的赤道海温平均水平。Niño 3.4指数通常使用5个月的运行平均值，当Niño 3.4 sst连续6个月或以上超过±0.4℃时，则定义El Niño或La Niña事件。

经验正交函数分析方法(Empirical Orthogonal Function,缩写为EOF)，也称特征向量分析(Eigenvector Analysis)，或者主成分分析(Principal Component Analysis,缩写PCA)，是一种分析矩阵数据中的结构特征，提取主要数据特征量的一种方法。Lorenz在1950年代首次将其引入气象和气候研究，在地学及水声学等其他学科中得到了非常广泛的应用。地学数据分析中通常特征向量对应的是空间样本，所以也称空间特征向量或者空间模态；主成分对应的是时间变化，也称时间系数。因此地学中也将EOF分析称为时空分解。

F=EOF×PC+MeanF=EOF\times PC+Mean

算法解释如下：

Input

time*(map)矩阵，行表示时间，列表示空间

Output

EOF空间模态，PC时间模态

Step

0️⃣ 预处理：经度处理

F=cos(lat)F=\sqrt{cos(lat)}

这一步的目的是计算协方差达到：

C=cos(lat)C=cos(lat)

的效果。

本质上是因为分带在不同维度上的面积会发生变化，需要通过角度进行修正。

1️⃣ 计算矩阵距平

这步是主成分分析的关键，不利用距平的话在数学意义上就无法通过协方差矩阵进行特征分解。

F=(F−Mean)F=(F-Mean)

2️⃣ 计算协方差矩阵

在相互独立的正交基上对原始矩阵进行空间变换，寻求方差最大的部分

C=cov(F)C=cov(F)

3️⃣ 计算特征值和特征向量

特征值和特征向量分别表示了正交基上的缩放倍数和正交基，我们选取方差上最大的缩放倍数，也就是最大的特征值来对原始空间进行变换，得到的新数据能够保留最大的信息量。

EOF,L=eig(C)EOF,L=eig(C)

这个特征向量就是我们的EOF空间模态。

4️⃣计算时间模态

PC=EOFT×(F)PC=EOF^T\times(F)

原始数据在新的向量空间上的变换，那就是主成分PC，在这里也就是我们的时间特征。

结果如下：

空间模态:

image-20221015205552316 image-20221015205611304

对应的时间特征：

image-20221015205651316 image-20221015205733048

我们选择Nino 3.4指数进行相关性判别：

PC1 PC2 NIno
PC1 1 4.17e-09 -0.142
PC2 4.17e-09 1 0.831
Nino -0.142 0.831 1

可以发现，PC2与NIno 3.4指数具有强相关性(≥0.8)(\ge0.8),可以用来探究厄尔尼诺现象。

我们也绘制了每个月的空间气候直方图，可以看到，有些月份的直方图具有明显的高峰和聚集程度。

GIF

年份箱线图变化趋势

例如：07年全球温度变化跨度最大，14年气温最高，98年气温比较集中且相对较低，数据量集中分布与Q1Q_1和中位数之间，以低温为主，高温为辅，但96年的时候，相对高温地区(高于Q2Q_2)较多。

下载

空间自相关指数

API代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
# 利用空间统计量Moran和Geary计算遥感数据的自相关程度
import numpy as np
import pandas as pd

def getMoranV(path,val=1,t=0,method="Moran"):
    '''
    :param path: 图像路径
    :param t: 卷积核类型
    :param method: 空间自相关方法
    :param val: 空间权重
    :return:
    '''
    # 读取数据
    wb1 = path
    data = wb1.values
    n, m = wb1.shape

    # 定义卷积核

    # queen 卷积
    dir2 = [[1, 0], [-1, 0], [0, 1], [0, -1], [1, 1], [1, -1], [-1, 1], [-1, -1]]
    # Rook 卷积
    dir1 = [[1, 0], [0, 1], [-1, 0], [0, -1]]
    # Bishop 卷积
    dir3=[[1,1],[1,-1],[-1,1],[-1,-1]]

    if t==0:
        dir=dir1
    elif t==1:
        dir=dir2
    else:
        dir=dir3

    # 开辟O(n^4)权重矩阵
    # 这玩意很耗内存，所以可以考虑用局部计算代替全局计算
    a = [[] for _ in range(n*m)]
    a_weight=0
    for i in range(n):
        for j in range(m):
            # 找到他在二维上的临界点
            loc=i*m+j
            for detX, detY in dir:
                x = i + detX
                y = j + detY
                # 转化为一维
                # 边界处理
                if 0 <= x < n and 0 <= y < m:
                    a[loc].append([x,y,val])
                    a_weight+=val

    def moranIndex(data,weight):
        # 莫兰指数计算
        # s_0=np.sum(np.sum(weight))
        s_0=a_weight
        n,m=len(data),len(data[0])
        x_hat=np.mean(data)
        up_sum=0
        down_sum=0
        for i in range(n):
            for j in range(m):
                # 下一轮
                loc=i*m+j
                for v in weight[loc]:
                    up_sum+=v[2]*(data[i][j]-x_hat)*(data[v[0]][v[1]]-x_hat)
                down_sum+=(data[i][j]-x_hat)**2
        return (n*m/s_0)*(up_sum/down_sum)

    def Geary_C(data,weight):
        # Geary's C计算
        s_0=a_weight*2
        n, m = len(data), len(data[0])
        x_hat = np.mean(data)
        up_sum = 0
        down_sum = 0
        for i in range(n):
            for j in range(m):
                # 下一轮
                loc = i * m + j
                for v in weight[loc]:

                    up_sum += v[2] * (data[i][j] - data[v[0]][v[1]])**2
                down_sum += (data[i][j] - x_hat) ** 2
        return ((n*m-1)/(s_0))*(up_sum/down_sum)

    if (M:=method.upper())=="MORAN":
        return moranIndex(data,a)
    elif M=="GEARY":
        return Geary_C(data,a)

def getHelp():
    print("Method: {'Moran','Geary'}")
    print("t: {0:Rook,1:Queen,2:bishop}")

if __name__ == '__main__':
    getHelp()
    path = r"\data.xlsx"
    print("墨兰指数为: %s"%getMoranV(path,val=1,t=0))
    print("Geary's C为: %s"%getMoranV(path,val=1,t=0,method="Geary"))

处理

逐月遍历数据即可

1
2
3
List=[]
for i in range(data.shape[0]): 	
    List.append(getMoranV(pd.DataFrame(data[i]),val=1,t=1))

结果可视化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
plt.figure(figsize=(20,6))
plt.plot(range((v:=data.shape[0])),List,'r-')
Mean=sum(List)/(len(List))
plt.plot(range(v),[Mean]*v,'b',linestyle='--',alpha=0.2)
plt.title("基于Queen邻接的SST空间自相关月变化")
plt.ylabel("时间")
plt.xlabel("日期")
T=[]
M,Y=1,1990
for i in range(360):
    T.append("%s年%s月"%(Y,M))
    if M%12==0:
        M=1
        Y+=1
    else:
        M+=1
plt.xticks(range(360)[::20],T[::20])
plt.show()

下载

绘制逐月动图

1
2
3
4
5
6
7
8
List=[[]]
t=0
while t<data.shape[0]:

    if len(List[-1])==12:
        List.append([])
    List[-1].append(getMoranV(pd.DataFrame(data[t]),val=1,t=1))
    t+=1

我们通过二维列表的方式对年份加以区分，这也是常见的手段。

1
2
3
4
5
6
7
8
9
10
11
12
13
# 逐年绘制
path=r"\yourpath"
timeList=[]
for i in range(len(List)):
    lis=List[i]
    plt.figure(figsize=(20,6))
    plt.plot(range(12),lis,'r-')
    Mean=sum(lis)/12
    plt.plot(range(12),[Mean]*12,'b',linestyle='--',alpha=0.2)
    plt.title("%s年"%(i+1990))
    plt.xticks(range(12),["%s月"%(m+1) for m in range(12)])
    plt.savefig((val:=path+"\\"+"%s.png"%i),dpi=200)
    timeList.append(val)

接着就是正常的绘制，我们需要将图像存储下来，并记录存储位置。

通过imageio模块制作动图。

1
2
3
4
5
6
IMG=[]
import imageio
path=r"C:\Users\lenovo\Desktop\Ocean\Spatial"
for i in timeList:
    IMG.append(imageio.imread(i))
imageio.mimsave(path+"\\"+"GIF1.gif",IMG,"GIF",duration=1)

GIF1