python气象信息系统工程 第十一章

第11章 常用气象统计方法与检验

11.1 基本气候状态统计量

11.1.1 中心趋势统计量

1
2
3
4
5
6
7
import numpy as np
station1 = np.array([27,29,28,26,24,25,28])
station2 = np.array([32,29,30,33,np.nan,33,36])
print(np.mean(station1))
print(np.mean(station2))
print(np.nanmean(station1))
print(np.nanmean(station2))
26.714285714285715
nan
26.714285714285715
32.166666666666664

1
2
3
4
5
6
7
import numpy as np
data = np.array([[1,3,6], [7,5,3]])
print(np.percentile(data, 50))
print(np.percentile(data, 50, axis=0, interpolation='lower'))
print(np.percentile(data, 50, axis=1))
print(np.percentile(data, 50, axis=1, keepdims=True))

4.0
[1 3 3]
[3. 5.]
[[3.]
 [5.]]

11.1.2 变化幅度统计量

1
2
3
4
5
import numpy as np
data = np.array([1,3,6,7,5,3])
data_anomaly = data - data.mean()
print(data_anomaly)

[-3.16666667 -1.16666667  1.83333333  2.83333333  0.83333333 -1.16666667]

1
2
3
4
5
6
7
import numpy as np
data = np.array([1,3,6,7,5,3,np.nan])
std1 =  np.std(data)
std2 =  np.nanstd(data)
print(std1)
print(std2)

nan
2.034425935955617

11.1.3 相关统计量

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
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from scipy.stats import pearsonr
import cartopy.mpl.ticker as cticker
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook4259/data.nc')
a = f['a']
b = f['b']
lat = f['lat']
lon = f['lon']
#创建两个大小为[nlat,nlon]的数组存放相关系数及p值
r = np.zeros((b.shape[1],b.shape[2]))
p = np.zeros((b.shape[1],b.shape[2]))
#循环计算相关系数
for i in np.arange(b.shape[1]):
   for j in np.arange(b.shape[2]):
      r[i,j],p[i,j] = pearsonr(b[:,i,j],a)
#绘图
fig = plt.figure(figsize=(15, 12))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
ax.set_extent([-90,30,0,90])
ax.set_xticks(np.arange(-90,60,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(0,120,30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制相关系数分布图
c1 = ax.contourf(lon, lat, r,levels=np.arange(-1,1.1,0.1),extend='both',cmap='bwr',zorder=1,transform=ccrs.PlateCarree())
#对通过90%显著性检验的区域打点，即将p值为0～0.1的区域覆盖上“…”标记
c2 = ax.contourf(lon,lat, p, [0,0.1,1] , zorder=1,hatches=['...', None],colors="none", transform=ccrs.PlateCarree())
#添加色标
fig.colorbar(c1,ax=ax,fraction=0.032)
Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1


/opt/conda/lib/python3.7/site-packages/scipy/stats/stats.py:3399: PearsonRConstantInputWarning: An input array is constant; the correlation coefficent is not defined.
  warnings.warn(PearsonRConstantInputWarning())





<matplotlib.colorbar.Colorbar at 0x7f0d56cee1d0>

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
import xarray as xr
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from scipy.stats import linregress
import cartopy.mpl.ticker as cticker
#读取数据
f = xr.open_dataset('/home/mw/input/pythonbook4259/data.nc')
a = f['a']
b = f['b']
lat = f['lat']
lon = f['lon']
#创建两个大小为[nlat,nlon]的数组存放相关系数及p值
slope = np.zeros((b.shape[1],b.shape[2]))
intercept = np.zeros((b.shape[1],b.shape[2]))
r = np.zeros((b.shape[1],b.shape[2]))
p = np.zeros((b.shape[1],b.shape[2]))
std_err = np.zeros((b.shape[1],b.shape[2]))#循环计算相关系数
for i in np.arange(b.shape[1]):
   for j in np.arange(b.shape[2]):
      slope[i,j],intercept[i,j],r[i,j],p[i,j],std_err[i,j] = linregress(b[:,i,j],a)
#绘图
fig = plt.figure(figsize=(15, 12))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
ax.set_extent([-90,30,0,90])
ax.set_xticks(np.arange(-90,60,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(0,120,30), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
#绘制相关系数分布图
c1 = ax.contourf(lon, lat, slope,levels=np.arange(-1,1.1,0.1),extend='both',cmap='bwr',zorder=1,transform=ccrs.PlateCarree())
#对通过90%显著性检验的区域打点，即将p值为0～0.1的区域覆盖上“…”标记
c2 = ax.contourf(lon,lat, p, [0,0.1,1] , zorder=1,hatches=['...', None],colors="none", transform=ccrs.PlateCarree())
#添加色标
fig.colorbar(c1,ax=ax,fraction=0.032)

/opt/conda/lib/python3.7/site-packages/scipy/stats/_stats_mstats_common.py:130: RuntimeWarning: invalid value encountered in double_scalars
  slope = r_num / ssxm
/opt/conda/lib/python3.7/site-packages/scipy/stats/_stats_mstats_common.py:142: RuntimeWarning: divide by zero encountered in double_scalars
  sterrest = np.sqrt((1 - r**2) * ssym / ssxm / df)





<matplotlib.colorbar.Colorbar at 0x7f0d5679f2d0>

11.1.4 数据标准化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np
#原始值
x = np.array([15.4,14.6,15.8,14.8,15.0,15.1,15.1,15.0,15.2,15.4,
            14.8,15.0,15.1,14.7,16.0,15.7,15.4,14.5,15.1,15.3,
            15.5,15.1,15.6,15.1,15.1,14.9,15.5,15.3,15.3,15.4,
            15.7,15.2,15.5,15.5,15.6,15.1,15.1,16.0,16.0,16.8,
            16.2,16.2,16.0,15.6,15.9,16.2,16.7,15.8,16.2,15.9,
            15.8,15.5,15.9,16.8,15.5,15.8,15.0,14.9,15.3,16.0,
            16.1,16.5,15.5,15.6,16.1,15.6,16.0,15.4,15.5,15.2,
            15.4,15.6,15.1,15.8,15.5,16.0,15.2,15.8,16.2,16.2,
            15.2,15.7,16.0,16.0,15.7,15.9,15.7,16.7,15.3,16.1])
#标准化
x_std = (x - x.mean(axis=0))/x.std(axis=0)
fig = plt.figure(figsize=(15,15))
#原始序列
ax1 = fig.add_axes([0.1, 0.1, 0.4, 0.3])
ax1.plot(np.arange(1900,1990,1),x,'k')
ax1.set_title('Original value')
#标准化后
ax2 = fig.add_axes([0.55, 0.1, 0.4, 0.3])
ax2.plot(np.arange(1900,1990,1),x_std,'k')
ax2.set_title('Normalized values')

Text(0.5, 1.0, 'Normalized values')

11.2 气候变化趋势分析

11.2.1 拟合

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import numpy as np
import matplotlib.pyplot as plt
#创建数据
x = np.arange(2000,2011,1)
y = np.array([3.9, 4.4, 10.8, 10.3, 11.2, 13.1, 14.1,  9.9, 13.9, 15.1, 12.5]).reshape(-1,1)
#创建回归模型
a1, b1 = np.polyfit(x, y, deg=1)  #对x、y进行拟合
y_fit1 = a1*x + b1  #拟合出的直线
a2, b2,c2 = np.polyfit(x, y, deg=2)  #对x、y进行拟合
y_fit2 = a2*x*x +b2*x+c2  #拟合出的直线
#绘图
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(1,2,1)
ax1.plot(x, y_fit1, '-')
ax1.plot(x, y, 'o', color='tab:brown')
ax1.text(2001,6,'y = {}x+{}'.format(a1.round(2)[0],b1.round(2)[0]),fontsize=12)
ax2 = fig.add_subplot(1,2,2)
ax2.plot(x, y_fit2, '-')
ax2.plot(x, y, 'o', color='tab:brown')
ax2.text(2001,6,'y = {}x+{}x+{}'.format(a2.round(2)[0],b2.round(2)[0],c2.round(2)[0]),fontsize=12)
plt.show()

11.2.2 滑动平均

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
import numpy as np
import matplotlib.pyplot as plt
#创建气候序列y
y = np.array([6.08, 4.56, 5.63, 5.31, 5.15, 5.44, 4.65, 4.24, 7.3,  5.86, 4.51, 6.28, 5.55, 5.35,
5.12, 4.76, 4.35, 3.76, 4.74, 5.55, 4.54, 5.74, 5.54, 3.67, 4.77, 4.9,  3.06, 3.9,
4.18, 5.44, 5.21, 3.86, 3.96, 4.47, 4.37, 4.86, 4.43, 3.63, 3.98, 3.94, 5.09, 4.48,
4.05, 4.81, 4.07, 4.48, 4.46, 3.95, 5.24, 3.54, 3.11, 5.07, 6.09, 4.59, 4.55, 4.7,
3.43, 4.37, 4.79, 3.64, 4.3,  3.5 ])
#创建一个9年滑动平均的卷积核
v= np.repeat(1/9, 9)
#滑动平均
y_smooth9_full = np.convolve(y,v,'full')
y_smooth9_same = np.convolve(y,v,'same')
y_smooth9_valid = np.convolve(y,v,'valid')
#绘图
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(2,2,1)
ax1.plot(np.arange(0,62,1),y,label='y')
ax1.set_xlim(0,60)
ax1.set_ylim(2,8)
ax1.set_title('Original value')
ax2 = fig.add_subplot(2,2,2)
ax2.plot(np.arange(-4,62+4,1),y_smooth9_full,label='full')
ax2.set_xlim(0,60)
ax2.set_ylim(2,8)
ax2.set_title('mode = full')
ax3 = fig.add_subplot(2,2,3)
ax3.plot(np.arange(0,62,1),y_smooth9_same,label='same')
ax3.set_xlim(0,60)
ax3.set_ylim(2,8)
ax3.set_title('mode = same')
ax4 = fig.add_subplot(2,2,4)
ax4.plot(np.arange(0+4,62-4,1),y_smooth9_valid,label='valid')
ax4.set_xlim(0,60)
ax4.set_ylim(2,8)
ax4.set_title('mode = valid')
#防止图题被遮挡
plt.tight_layout()
plt.show()

11.2.3 去趋势

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from scipy import signal
import numpy as np
y = np.array([6.08, 4.56, 5.63, 5.31, 5.15, 5.44, 4.65, 4.24, 7.3,  5.86, 4.51, 6.28, 5.55, 5.35,
5.12, 4.76, 4.35, 3.76, 4.74, 5.55, 4.54, 5.74, 5.54, 3.67, 4.77, 4.9,  3.06, 3.9,
4.18, 5.44, 5.21, 3.86, 3.96, 4.47, 4.37, 4.86, 4.43, 3.63, 3.98, 3.94, 5.09, 4.48,
4.05, 4.81, 4.07, 4.48, 4.46, 3.95, 5.24, 3.54, 3.11, 5.07, 6.09, 4.59, 4.55, 4.7,
3.43, 4.37, 4.79, 3.64, 4.3,  3.5 ])
y_detrend_linear = signal.detrend(y,axis = 0,type='linear')
y_detrend_constant = signal.detrend(y,axis = 0,type='constant')
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(1,2,1)
ax1.plot(np.arange(0,62,1),y,'k',label='y')
ax1.plot(np.arange(0,62,1),y_detrend_linear+y.mean(),'r',label='linear')
ax1.legend()
ax2 = fig.add_subplot(1,2,2)
ax2.plot(np.arange(0,62,1),y,'k',label='y')
ax2.plot(np.arange(0,62,1),y_detrend_constant,'r',label='constant')
ax2.legend()
plt.show()


11.2.4 滤波

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
from scipy import signal
import numpy as np
y = np.array([ 1.74082189, -0.12322827,  1.18896494,  0.79653332,  0.60031752,  0.95595867,
 -0.01285688, -0.51565989,  3.23696743,  1.47102516, -0.18454571,  1.98609165,
  1.09085703,  0.84558727,  0.56352705,  0.12204148, -0.38076152, -1.10430731,
  0.09751451,  1.09085703, -0.14775525,  1.3238633,  1.07859355,  -1.2146787,
  0.13430497,  0.29373032, -1.96275147, -0.93261848, -0.58924082,  0.95595867,
  0.67389844, -0.98167243, -0.85903755, -0.23359967, -0.35623455,  0.24467636,
 -0.28265362, -1.26373265, -0.83451058, -0.88356453,  0.52673659, -0.22133618,
 -0.74866616,  0.18335892, -0.72413918, -0.22133618, -0.24586315, -0.87130104,
  0.71068891, -1.37410405, -1.90143403,  0.50220961,  1.75308538, -0.08643781,
 -0.13549176,  0.04846056, -1.50900241, -0.35623455,  0.15883195, -1.25146917,
 -0.44207896, -1.423158  ])
####高通滤波####
b1, a1 = signal.butter(3, 2/3, 'highpass')
high_series = signal.filtfilt(b1, a1, y)
####带通滤波####获得3～10年的信号
b2, a2 = signal.butter(3, [2/10,2/3], 'bandpass')
band_series = signal.filtfilt(b2, a2, y)
####低通滤波####获得10年以上的信号
b3, a3 = signal.butter(3, 2/10, 'lowpass')
low_series = signal.filtfilt(b3, a3, y)
#绘图
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(1,3,1)
ax1.plot(np.arange(0,62,1),y,'k',label='y')
ax1.plot(np.arange(0,62,1),high_series,'r',label='highpass')
ax1.legend()
ax2 = fig.add_subplot(1,3,2)
ax2.plot(np.arange(0,62,1),y,'k',label='y')
ax2.plot(np.arange(0,62,1),band_series,'r',label='bandpass')
ax2.legend()
ax3 = fig.add_subplot(1,3,3)
ax3.plot(np.arange(0,62,1),y,'k',label='y')
ax3.plot(np.arange(0,62,1),low_series,'r',label='lowpass')
ax3.legend()
plt.show()

11.3 气候序列突变检验

11.3.1 滑动t检验

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
import numpy as np
import matplotlib.pyplot as plt
def slidet(inputdata,step): #inputdata为输入的一维气候序列；step为滑动步长，整型
   inputdata = np.array(inputdata)
   n = inputdata.shape[0]
   t = np.zeros(n)
   t1 = np.empty(n)
   n1 = step   
   n2 = step
   n11 = 1 / n1
   n22 = 1 / n2
   m = np.sqrt(n11 + n22)
   for i in range (step, n-step-1):
      x1_mean = np.mean(inputdata[i-step : i])   
      x2_mean = np.mean(inputdata[i : i+step])
      s1 = np.var(inputdata[i-step : i])        
      s2 = np.var(inputdata[i : i+step])
      s = np.sqrt((n1 * s1 + n2 * s2) / (n1 + n2 - 2))
      t[i-step] = (x2_mean - x1_mean) / (s * m)
      t1 = np.roll(t , step-1)
      t1[:step]=np.nan  
      t1[n-step+1:]=np.nan          
   return  t1  
#示例
a = [15.4,14.6,15.8,14.8,15.0,15.1,15.1,15.0,15.2,15.4,
     14.8,15.0,15.1,14.7,16.0,15.7,15.4,14.5,15.1,15.3,
     15.5,15.1,15.6,15.1,15.1,14.9,15.5,15.3,15.3,15.4,
     15.7,15.2,15.5,15.5,15.6,15.1,15.1,16.0,16.0,16.8,
     16.2,16.2,16.0,15.6,15.9,16.2,16.7,15.8,16.2,15.9,
     15.8,15.5,15.9,16.8,15.5,15.8,15.0,14.9,15.3,16.0,
     16.1,16.5,15.5,15.6,16.1,15.6,16.0,15.4,15.5,15.2,
     15.4,15.6,15.1,15.8,15.5,16.0,15.2,15.8,16.2,16.2,
     15.2,15.7,16.0,16.0,15.7,15.9,15.7,16.7,15.3,16.1,16.2]
t = slidet(a,11)
#绘制图形
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(np.arange(1900,1991,1),t,'r')
ax1.axhline(2.85) #显著性水平可通过查阅对应自由度的t分布表得到
ax1.axhline(-2.85)
plt.show()

11.3.2 曼-肯德尔法

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
def mktest(x):
    n = len(x)
    #定义累计量序列Sk，长度n，初始值为0
    Sk = np.zeros(n)
    s = 0    
    #计算正序列UFK
    UF = np.zeros(n)
    for i in range(2,n):
        for j in range(1,i):
            if x[i]>x[j]:
                s += 1
        Sk[i] = s
        E = i * (i - 1)/4
        Var = i * (i - 1) * (2 * i + 5)/72
        UF[i] = (Sk[i] - E)/np.sqrt(Var)
    #定义逆序列累计量序列Sk2，长度n，初始值为0
    Sk2 = np.zeros(n)
    s = 0
    #计算逆序列，长度为n，初值为0
    x2 = x[::-1]    
    UB = np.zeros(n)
    for i in range(2, n):
        for j in range(1, i):
            if x2[i] > x2[j]:
                s += 1
        Sk2[i] = s
        E = i * (i - 1) / 4
        Var = i * (i - 1) * (2 * i + 5) / 72
        UB[i] = -(Sk2[i] - E) / np.sqrt(Var)
    UB = UB[::-1]
    return UF, UB

import numpy as np
import matplotlib.pyplot as plt

x = [15.4,14.6,15.8,14.8,15.0,15.1,15.1,15.0,15.2,15.4,
     14.8,15.0,15.1,14.7,16.0,15.7,15.4,14.5,15.1,15.3,
     15.5,15.1,15.6,15.1,15.1,14.9,15.5,15.3,15.3,15.4,
     15.7,15.2,15.5,15.5,15.6,15.1,15.1,16.0,16.0,16.8,
     16.2,16.2,16.0,15.6,15.9,16.2,16.7,15.8,16.2,15.9,
     15.8,15.5,15.9,16.8,15.5,15.8,15.0,14.9,15.3,16.0,
     16.1,16.5,15.5,15.6,16.1,15.6,16.0,15.4,15.5,15.2,
     15.4,15.6,15.1,15.8,15.5,16.0,15.2,15.8,16.2,16.2,
     15.2,15.7,16.0,16.0,15.7,15.9,15.7,16.7,15.3,16.1,16.2]


uf,ub = mktest(x)
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(np.arange(1900,1991,1),uf,'r',label='UFk')
ax1.plot(np.arange(1900,1991,1),ub,'b',label='UBk')
ax1.legend()
ax1.axhline(1.96)
ax1.axhline(-1.96)
plt.show()

11.4 气候变量场时空结构的分离（经验正交函数分解）

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
from eofs.standard import Eof
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
#读入数据
f_u = xr.open_dataset('/home/mw/input/pythonbook4259/uwnd.mon.mean.nc')
u = f_u['uwnd'].loc[f_u.time.dt.month.isin([12,1,2])].loc['1961-12-01':'2017-02-28'].loc[:,1000:100,70:20,60:120]
lon_u = u['lon']
lat_u = u['lat']
level = u['level']
u_array = np.array(u).reshape(56,3,level.shape[0],lat_u.shape[0],lon_u.shape[0]).mean((1,4)).transpose((0,2,1))
#EOF分析
coslat = np.cos(np.deg2rad(lat_u))
coslat = np.array(coslat)
wgts = np.sqrt(coslat)[..., np.newaxis]
solver = Eof(u_array, weights=wgts)
u_eof = solver.eofsAsCorrelation(neofs=2)
u_pc = solver.pcs(npcs=2, pcscaling=1)
u_var = solver.varianceFraction()
u_eof = u_eof.transpose((0,2,1))
#设置PC序列的颜色
color1 = []
color2 = []
for i in range(1961,2017):
   if u_pc[i-1961,0] >=0:
      color1.append('red')
   elif u_pc[i-1961,0] <0:
      color1.append('blue')
   if u_pc[i-1961,1] >=0:
      color2.append('red')
   elif u_pc[i-1961,1] <0:
      color2.append('blue')   
#绘图    
fig1 = plt.figure(figsize=(15,15))
years = range(1961, 2017)
#EOF1
f1_ax1 = fig1.add_axes([0.1, 0.42, 0.3, 0.2])
f1_ax1.set_title('(a) EOf1',loc='left')
f1_ax1.set_title( '%.2f%%' % (u_var[0]*100),loc='right')
f1_ax1.set_yscale('symlog')
f1_ax1.set_xticklabels([r'20$^\degree$',r'30$^\degree$N', r'40$^\degree$N',
                    r'50$^\degree$N',r'60$^\degree$N', r'70$^\degree$N'])
f1_ax1.set_yticks([1000, 500,300, 200, 100])
f1_ax1.set_yticklabels(['1000','500','300','200','100'])
f1_ax1.invert_yaxis()
f1_ax1.set_ylabel('Height (hPa)',fontsize=18)
f1_ax1.set_xlabel('Latitude',fontsize=18)
c1 = f1_ax1.contourf(lat_u,level ,u_eof[0,:,:], levels=np.arange(-0.8,0.9,0.1), extend = 'both',
             zorder=0, cmap=plt.cm.RdBu_r)
#PC1
f1_ax2 = fig1.add_axes([0.45, 0.42, 0.3, 0.2])
f1_ax2.set_title('(b) PC1',loc='left')
f1_ax2.set_ylim(-3,3)
f1_ax2.axhline(0,linestyle="--")
f1_ax2.bar(years,u_pc[:,0],color=color1)
#EOF2
f1_ax3 = fig1.add_axes([0.1, 0.17, 0.3, 0.2])
f1_ax3.set_title('(c) EOf1',loc='left')
f1_ax3.set_title( '%.2f%%' % (u_var[1]*100),loc='right')
f1_ax3.set_yscale('symlog')
f1_ax3.set_xticklabels([r'20$^\degree$',r'30$^\degree$N', r'40$^\degree$N',
                    r'50$^\degree$N',r'60$^\degree$N', r'70$^\degree$N'])
f1_ax3.set_yticks([1000, 500,300, 200, 100])
f1_ax3.set_yticklabels(['1000','500','300','200','100'])
f1_ax3.invert_yaxis()
f1_ax3.set_ylabel('Height (hPa)',fontsize=18)
f1_ax3.set_xlabel('Latitude',fontsize=18)
c2 = f1_ax3.contourf(lat_u,level ,u_eof[1,:,:], levels=np.arange(-0.8,0.9,0.1), extend = 'both',
             zorder=0, cmap=plt.cm.RdBu_r)
#PC2
f1_ax4 = fig1.add_axes([0.45, 0.17, 0.3, 0.2])
f1_ax4.set_title('(d) PC2',loc='left')
f1_ax4.set_ylim(-3,3)
f1_ax4.axhline(0,linestyle="--")
f1_ax4.bar(years,u_pc[:,1],color=color2)
#色标
position=fig1.add_axes([0.1, 0.1, 0.3, 0.017])
fig1.colorbar(c1,cax=position,orientation='horizontal',format='%.1f',)
plt.show()

/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:42: UserWarning: FixedFormatter should only be used together with FixedLocator
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:62: UserWarning: FixedFormatter should only be used together with FixedLocator
/小书匠/python