Python气象数据处理与绘图：经验正交函数(EOF)与旋转经验正交函数(REOF)

1、导入模块

1from pyEOF import *
2import xarray as xr
3import numpy as np
4import pandas as pd
5import matplotlib.pyplot as plt

2、定义绘图函数

 1# create a function for visualization convenience
 2def visualization(da, pcs, eofs_da, evf):
 3    fig = plt.figure(figsize = (12,12))
 4
 5    ax = fig.add_subplot(n+1,2,1)
 6    da.mean(dim=["lat","lon"]).plot(ax=ax)
 7    ax.set_title("average air temp")
 8
 9    ax = fig.add_subplot(n+1,2,2)
10    da.mean(dim="time").plot(ax=ax)
11    ax.set_title("average air temp")
12
13    for i in range(1,n+1):
14        pc_i = pcs["PC"+str(i)].to_xarray()
15        eof_i = eofs_da.sel(EOF=i)["air"]
16        frac = str(np.array(evf[i-1]*100).round(2))
17
18        ax = fig.add_subplot(n+1,2,i*2+1)
19        pc_i.plot(ax=ax)
20        ax.set_title("PC"+str(i)+" ("+frac+"%)")
21
22        ax = fig.add_subplot(n+1,2,i*2+2)
23        eof_i.plot(ax=ax,
24                   vmin=-0.75, vmax=0.75, cmap="RdBu_r",
25                   cbar_kwargs={'label': ""})
26        ax.set_title("EOF"+str(i)+" ("+frac+"%)")
27
28    plt.tight_layout()
29    plt.show()

3、读取气温数据

 1# load the DataArray
 2da = xr.tutorial.open_dataset('air_temperature')["air"]
 3print(da)
 4
 5# create a mask
 6mask = da.sel(time=da.time[0])
 7mask = mask.where(mask<250).isnull().drop("time")
 8
 9# get the DataArray with mask
10da = da.where(mask)
11da.sel(time=da.time[99]).plot()
12plt.show()

1# convert DataArray to DataFrame
2df = da.to_dataframe().reset_index() # get df from da
3display(df.head(5))
4print("DataFrame Shape:",df.shape)

1# reshape the dataframe to be [time, space]
2df_data = get_time_space(df, time_dim = "time", lumped_space_dims = ["lat","lon"])
3display(df_data.head(5))
4print("DataFrame Shape:",df_data.shape)

4、REOF分解

 1n = 4
 2pca = df_eof(df_data,pca_type="varimax",n_components=n)
 3
 4eofs = pca.eofs(s=2, n=n) # get eofs
 5eofs_da = eofs.stack(["lat","lon"]).to_xarray() # make it convenient for visualization
 6pcs = pca.pcs(s=2, n=n) # get pcs
 7evfs = pca.evf(n=n) # get variance fraction
 8
 9# plot
10visualization(da, pcs, eofs_da, evfs)

5、传统EOF分解

 1n = 4 # define the number of components
 2
 3pca = df_eof(df_data) # implement EOF
 4
 5eofs = pca.eofs(s=2, n=n) # get eofs
 6eofs_da = eofs.stack(["lat","lon"]).to_xarray() # make it convenient for visualization
 7pcs = pca.pcs(s=2, n=n) # get pcs
 8evfs = pca.evf(n=n) # get variance fraction
 9
10# plot
11visualization(da, pcs, eofs_da, evfs)

6、eofs模块分解结果

 1from eofs.standard import Eof
 2from sklearn.preprocessing import StandardScaler
 3solver = Eof(StandardScaler().fit_transform(df_data.values))
 4
 5s_pcs = pd.DataFrame(data=solver.pcs(npcs=4, pcscaling=2),
 6                     columns = pcs.columns,
 7                     index = pcs.index)
 8
 9
10s_eofs = pd.DataFrame(data = solver.eofs(neofs=4, eofscaling=2),
11                      columns = eofs.columns,
12                      index = eofs.index)
13s_eofs_da = s_eofs.stack(["lat","lon"]).to_xarray() # make it convenient for visualization
14
15s_evfs = solver.varianceFraction(neigs=4)
16
17# plot
18visualization(da, s_pcs, s_eofs_da, s_evfs)