1
import numpy as np
2
import xarray as xr
3
import matplotlib.pyplot as plt
4
import cartopy.crs as ccrs
5
from sacpy import EOF
6

7
# 1. 数据准备
8
da = xr.open_dataset('sst.nc')['sst']
9
da_anom = da.groupby('time.month') - da.groupby('time.month').mean('time')
10

11
# 2. EOF 分析
12
coslat = np.cos(np.deg2rad(da_anom['lat'].values))
13
eof = EOF(da_anom, weights=np.sqrt(coslat))
14
eof.solve()
15
vf = eof.get_varperc()
16
print(f"EOF1 方差解释: {vf[0]:.2%}")
17

18
# 3. 获取结果
19
pc1 = eof.get_pc(npt=1, scaling="std")       # (1, time)
20
pt1 = eof.get_pt(npt=1, scaling="mstd")       # (1, lat, lon)
21
corr_map = eof.correlation_map(npt=1)          # (1, lat, lon)
22

23
# 4. 绘图（使用 sacpy 的绘图扩展）
24
import sacpy.Map  # 注册 GeoAxesSubplot 扩展方法
25

26
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=190))
27
c = ax.scontourf(da_anom['lon'], da_anom['lat'], corr_map[0],
28
                 cmap='RdBu_r')
29
ax.init_map(coastlines=True)
30
plt.colorbar(c, label='Correlation')
31
plt.title('EOF1 Correlation Map')
32
plt.show()

1
import numpy as np
2
import xarray as xr
3
import matplotlib.pyplot as plt
4
import cartopy.crs as ccrs
5
from sacpy import SVD, get_anom
6

7
# 1. 数据准备
8
sst = xr.open_dataset('sst.nc')['sst']
9
precip = xr.open_dataset('precip.nc')['precip']
10

11
sst_anom = get_anom(sst)         # (time, lat, lon)
12
precip_anom = get_anom(precip)   # (time, lat, lon)
13

14
# 2. SVD 分析
15
svd = SVD(sst_anom, precip_anom)
16
svd.solve()
17

18
# 3. 获取结果
19
pc_left, pc_right = svd.get_pc(npt=3, norm="std")
20
hetero_left, hetero_right = svd.get_heterogeneous_map(npt=2)
21

22
# 4. 绘图 - 第1 SVD 模态的非齐次相关图
23
import sacpy.Map
24

25
fig, axes = plt.subplots(1, 2, subplot_kw={
26
    'projection': ccrs.PlateCarree(central_longitude=190)
27
})
28

29
ax = axes[0]
30
c = ax.scontourf(sst_anom['lon'], sst_anom['lat'], hetero_left[0],
31
                 cmap='RdBu_r', extent='both')
32
ax.init_map(coastlines=True)
33
ax.set_title('Heterogeneous Map (SST ← Precip PC)')
34

35
ax = axes[1]
36
c = ax.scontourf(precip_anom['lon'], precip_anom['lat'], hetero_right[0],
37
                 cmap='RdBu_r', extent='both')
38
ax.init_map(coastlines=True)
39
ax.set_title('Heterogeneous Map (Precip ← SST PC)')
40

41
plt.show()

1
import numpy as np
2
import xarray as xr
3
import matplotlib.pyplot as plt
4
import cartopy.crs as ccrs
5
from sacpy import LinReg, MultLinReg
6
import sacpy.Map
7

8
# 一元回归 + 绘图
9
da = xr.open_dataset('sst.nc')['sst']
10
enso = np.random.randn(da.shape[0])
11

12
reg = LinReg(enso, da)
13
reg.mask(0.05)
14

15
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=190))
16
c = ax.scontourf(da['lon'], da['lat'], reg.slope1, cmap='RdBu_r')
17
ax.sig_plot(da['lon'], da['lat'], reg.p_value, marker='..')
18
ax.init_map(coastlines=True)
19
plt.colorbar(c, label='Regression slope')
20
plt.title('ENSO Regression on SST (p<0.05 dotted)')
21
plt.show()

1
import numpy as np
2
import xarray as xr
3
import matplotlib.pyplot as plt
4
import cartopy.crs as ccrs
5
from sacpy import STMV, get_anom
6
import sacpy.Map
7

8
# 数据
9
sst = xr.open_dataset('sst.nc')['sst']
10
early = sst.sel(time=slice('1950', '1980'))
11
late = sst.sel(time=slice('1985', '2015'))
12

13
# 双样本 T 检验
14
stmv = STMV(late, early)
15

16
# 绘图
17
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=190))
18
c = ax.scontourf(sst['lon'], sst['lat'], stmv.mean, cmap='RdBu_r')
19
ax.sig_plot(sst['lon'], sst['lat'], stmv.p_value, marker='..')
20
ax.init_map(coastlines=True)
21
plt.colorbar(c, label='SST Difference (late - early)')
22
plt.title('SST Change between 1950-1980 and 1985-2015')
23
plt.show()

1
import numpy as np
2
import matplotlib.pyplot as plt
3
from sacpy import Spectral
4

5
# 准备数据
6
ts = sst_anom.mean(dim=('lat', 'lon')).values  # 区域平均 SSTA
7

8
# 谱分析
9
spec = Spectral(ts, M_length=64, remove_annual=True)
10

11
# 计算谱
12
f, Pxx = spec.get_spectra()
13
red = spec.get_red_noise()
14
sig = spec.get_rs_sig(p_crit=0.95)
15

16
# 绘图
17
plt.figure(figsize=(10, 4))
18
plt.semilogx(1/f, f*Pxx, 'b-', linewidth=1.5, label='Power Spectrum')
19
plt.semilogx(1/f, f*red, 'r--', label='Red Noise')
20
plt.semilogx(1/f, f*sig, 'r:', linewidth=0.5, label='95% Significance')
21
plt.xlabel('Period')
22
plt.ylabel('Frequency x Power')
23
plt.legend()
24
plt.title('Power Spectrum')
25
plt.show()

1
import numpy as np
2
import matplotlib.pyplot as plt
3
from sacpy import CrossSpectral
4

5
# 两个时间序列
6
ts1 = enso_index
7
ts2 = precip_index
8

9
cross = CrossSpectral(ts1, ts2, M_length=64)
10
f, Pxx, Pyy = cross.get_spectra()
11
cov, phase, coh = cross.get_cross_spectra()
12
coh_sig = cross.get_coh_sig(p_crit=0.95)
13

14
fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
15

16
axes[0].semilogx(1/f, f*Pxx, label='x')
17
axes[0].semilogx(1/f, f*Pyy, label='y')
18
axes[0].set_ylabel('Power')
19
axes[0].legend()
20

21
axes[1].semilogx(1/f, coh)
22
axes[1].axhline(coh_sig, color='r', linestyle='--', label='95% Sig')
23
axes[1].set_ylabel('Coherence')
24
axes[1].set_ylim(0, 1)
25

26
axes[2].semilogx(1/f, phase)
27
axes[2].set_ylabel('Phase (°)')
28
axes[2].set_xlabel('Period')
29

30
plt.suptitle('Cross-Spectral Analysis')
31
plt.show()

1
import numpy as np
2
import xarray as xr
3
import matplotlib.pyplot as plt
4
import cartopy.crs as ccrs
5
import cartopy.feature as cfeature
6
import sacpy.Map
7
from sacpy.Map import get_levels
8

9
# 准备数据
10
da = xr.open_dataset('sst.nc')['sst'].isel(time=0)  # 2D 场
11
lon, lat = da['lon'].values, da['lat'].values
12

13
# 创建地图
14
proj = ccrs.PlateCarree(central_longitude=180)
15
fig, axes = plt.subplots(2, 2, figsize=(14, 8),
16
                         subplot_kw={'projection': proj})
17
axes = axes.flatten()
18

19
# 子图1：智能填色 + 陆地 + 显著性打点
20
ax = axes[0]
21
c = ax.scontourf(lon, lat, da.values, cmap='RdBu_r',
22
                 percentile=98, num_level=13)
23
ax.add_feature(cfeature.LAND, facecolor='lightgray')
24
ax.init_map(coastlines=False)
25
plt.colorbar(c, ax=ax, label='SST')
26

27
# 子图2：pcolormesh
28
ax = axes[1]
29
m = ax.spcolormesh(lon, lat, da.values, cmap='coolwarm')
30
ax.init_map()
31
plt.colorbar(m, ax=ax)
32

33
# 子图3：使用 DataArray.splot
34
ax = axes[2]
35
m, _ = da.splot(ax=ax, kw1={'cmap': 'RdBu_r'})
36
ax.set_title('splot method')
37

38
# 子图4：不同投影
39
ax = axes[3]
40
proj2 = ccrs.Orthographic(central_longitude=180, central_latitude=0)
41
ax.remove()
42
ax = fig.add_subplot(2, 2, 4, projection=proj2)
43
c = ax.scontourf(lon, lat, da.values, cmap='RdBu_r')
44
ax.coastlines()
45
ax.set_global()
46
plt.colorbar(c, ax=ax, shrink=0.6)
47

48
plt.tight_layout()
49
plt.show()

1
import xarray as xr
2
import numpy as np
3
from sacpy import get_anom, spec_moth_dat, spec_moth_yrmean, convert_lon
4

5
# 加载数据
6
ds = xr.open_dataset('era5_monthly.nc')
7

8
# 经度转换（0-360 → -180~180）
9
ds = ds.assign_coords(lon=convert_lon(ds['lon'].values)).sortby('lon')
10

11
# 计算距平
12
t_anom = get_anom(ds['t2m'], method=0, freq='month')
13

14
# 提取冬季并计算年平均值
15
djf_t = spec_moth_yrmean(t_anom, "DJF")  # 每年冬季平均
16

17
# 提取夏季
18
jja_t = get_anom(spec_moth_dat(ds['t2m'], "JJA"))
19

20
print(djf_t)