1
from netCDF4 import Dataset, num2date
2

3
f = Dataset('sst.mnmean.nc')          # 打开文件
4
sst = f.variables['sst'][:]            # 读取变量数组
5
lon = f.variables['lon'][:]            # 读取经度
6
lat = f.variables['lat'][:]            # 读取纬度
7
time = f.variables['time']             # 获取时间变量对象
8
time = num2date(time[:], time.units, calendar='standard')
9
                                        # 将数值时间转为datetime对象

1
import xarray as xr
2

3
ds = xr.open_dataset('sst.mnmean.nc')  # 打开数据集
4
sst = ds['sst']                         # 获取 DataArray
5
lon = ds['lon']                         # 坐标自动关联
6
lat = ds['lat']

1
# 按时间范围选取（字符串切片）
2
sst = ds['sst'].sel(time=slice('1950', '2020'))
3

4
# 按具体日期选取
5
sst.loc['1997-12-01', :, :]            # 选取某一时刻
6

7
# 按时间范围选取（等间距切片）
8
nino34.loc['1950-12-01':'2019-12-01':12]  # 每隔12步取一次

1
# 按经纬度范围选取
2
nino34 = ds.sst.loc[:, 5:-5, 190:240]   # lat: 5S~5N, lon: 190~240
3
TNA    = ds.sst.loc[:, 25:5, 290:350]   # lat: 5N~25N, lon: 290~350
4
sst_trop = sst.loc[30:-30, 120:300]     # 热带区域

1
hgt500 = hgt.sel(level=500)   # 选择 500 hPa 层

1
# 等间距抽取（每12个月取一次，从第12个元素开始）
2
nino34_DJF = (np.array(nino34[11:-11:12]) +   # 12月 (索引11开始)
3
              np.array(nino34[12::12]) +       # 1月  (索引12开始)
4
              np.array(nino34[13::12])) / 3    # 2月  (索引13开始)

1
z55_20 = hgt500_jan.sel(lat=55, lon=20, method='nearest')

1
# 插值到新时间坐标（将 hgt 插值到 index 的时间轴上）
2
hgt500 = hgt500.interp(time=index_da['time'])
3

4
# 按精确时间轴对齐两个 DataArray
5
hgt_winter, wp_winter = xr.align(hgt_winter, wp_winter, join='exact')
6

7
# 取时间交集
8
start_time = max(wp_da['time'].values[0], hgt500['time'].values[0])
9
end_time   = min(wp_da['time'].values[-1], hgt500['time'].values[-1])
10
wp_da   = wp_da.sel(time=slice(start_time, end_time))
11
hgt500  = hgt500.sel(time=slice(start_time, end_time))

1
# Nino3.4 指数：5°S–5°N, 170°W–120°W (即 190°E–240°E)
2
nino34 = sst.loc[:, 5:-5, 190:240].mean(dim=['lat', 'lon'])
3

4
# Nino1+2 指数：0°–10°S, 90°W–80°W (即 270°E–280°E)
5
nino12 = sst.loc[:, 10:0, 270:280].mean(dim=['lat', 'lon'])
6

7
# TNA 指数（热带北大西洋）：5°N–25°N, 70°W–10°W (即 290°E–350°E)
8
TNA = sst.loc[:, 25:5, 290:350].mean(dim=['lat', 'lon'])

1
clim_mean = nino34.groupby('time.month').mean('time')  # 逐月气候态均值
2
clim_max  = nino34.groupby('time.month').max('time')   # 逐月最大值
3
clim_min  = nino34.groupby('time.month').min('time')   # 逐月最小值
4

5
yerr_upper = clim_max - clim_mean   # 误差条上限
6
yerr_lower = clim_mean - clim_min   # 误差条下限

1
# 方法一：手动索引（推荐，最清晰）
2
nino34_DJF = (np.array(nino34.loc['1950-12-01':'2019-12-01':12]) +
3
              np.array(nino34.loc['1951-01-01':'2020-01-01':12]) +
4
              np.array(nino34.loc['1951-02-01':'2020-02-01':12])) / 3
5

6
# 方法二：按位置索引（需假设数据从1月开始，无缺测）
7
nino34_DJF = (np.array(nino34[11:-11:12]) +   # 12月
8
              np.array(nino34[12::12]) +       # 1月
9
              np.array(nino34[13::12])) / 3    # 2月
10

11
# 方法三：混合 groupby（取不同年份的月份再平均）
12
nino34_DJF_2 = (np.array(nino34.loc['1950':'2019'].groupby('time.month')[12]) +
13
                np.array(nino34.loc['1951':'2020'].groupby('time.month')[1]) +
14
                np.array(nino34.loc['1951':'2020'].groupby('time.month')[2])) / 3

1
# xarray 方式
2
sst_mean = sst.mean(dim='time')
3

4
# numpy 方式
5
sst = np.mean(sst, axis=0)   # 沿第0轴（时间）求均值

1
coslat = np.cos(np.deg2rad(lat)).to_numpy()
2
wgts = np.sqrt(coslat)[..., np.newaxis]   # 权重矩阵 (lat, 1)
3

4
solver = Eof(ssta, weights=wgts)           # 传入权重

1
zonal_mean = hgt.mean(dim='lon')   # 结果维度 (time, lat)

1
index_std = (index - index.mean(dim='time')) / index.std(dim='time')

1
# 组合指数计算（如 EU 遥相关指数）
2
EU_index = -0.25 * z55_20 + 0.50 * z55_75 - 0.25 * z40_145
3

4
# 距平也可以直接用减法
5
hgt_anom = hgt - hgt.mean(dim='time')
6

7
# 纬向偏差
8
zonal_anom = hgt - hgt.mean(dim='lon')

1
import numpy as np
2

3
def LinearRegression(var, ts):
4
    s = var.shape                                    # (time, lat, lon)
5
    var = np.array(var).reshape(s[0], -1)            # → (time, npts)
6
    mask = np.isnan(np.mean(var, axis=0))            # 掩膜 NaN 点位
7
    X = np.vstack((ts, np.ones(len(ts)))).T          # 设计矩阵 (time, 2)
8
    r, *res = np.linalg.lstsq(X, var[:, ~mask], rcond=None)  # 最小二乘
9
    reg_map = np.zeros((2, s[1] * s[2])) * np.nan
10
    reg_map[:, ~mask] = r                            # 结果回填
11
    reg_map = np.reshape(reg_map, (2, s[1], s[2]))   # → (2, lat, lon)
12
    return reg_map                                   # [0]=斜率, [1]=截距

1
import sacpy as scp
2

3
LinReg = scp.LinReg(nino34, sst)
4
slope = LinReg.slope          # 回归系数场
5
p_value = LinReg.p_value      # p 值场

1
X = np.vstack((nino34, TNA, np.ones(len(nino34)))).T   # 两个自变量
2
r, *res = np.linalg.lstsq(X, y[:, ~mask], rcond=None)
3
# r[0] = Nino34 的回归系数, r[1] = TNA 的回归系数, r[2] = 截距

1
from scipy import stats
2

3
def linregress_1d(y, x):
4
    a, b, r, p, std = stats.linregress(x, y)
5
    return np.array([a, b, r, p, std])

1
from scipy import stats
2

3
el_nino_years = [1983, 1987, 1992, 1995, 1998, 2010, 2016]
4
la_nina_years = [1985, 1989, 1999, 2000, 2008, 2011]
5

6
# 提取夏季数据
7
summer_months = [6, 7, 8]
8
summer_t2m = t2m.sel(time=t2m['time.month'].isin(summer_months)) \
9
                .groupby('time.year').mean(dim='time')
10

11
# 挑选年份并计算合成平均
12
el_nino_composite = summer_t2m.sel(year=el_nino_years).mean(dim='year')
13
la_nina_composite = summer_t2m.sel(year=la_nina_years).mean(dim='year')

1
from scipy import stats
2

3
_, p_values = stats.ttest_1samp(composite_data, popmean=0.0, axis=0)
4

5
alpha = 0.05
6
significant = p_values < alpha

1
_, p_values = stats.ttest_ind(group1, group2, axis=0)

1
lon2d, lat2d = np.meshgrid(lon, lat)
2
ax.scatter(lon2d[significant], lat2d[significant],
3
           color='grey', s=2, transform=ccrs.PlateCarree())

1
sst = sst.to_numpy()
2
y = np.reshape(sst, (sst.shape[0], -1))   # (time, lat*lon)
3
# 回归计算后...
4
map = np.reshape(map, (sst.shape[1], sst.shape[2]))   # 恢复为 (lat, lon)

1
mask = np.isnan(np.mean(y, axis=0))       # 找出含 NaN 的格点
2
r, *res = np.linalg.lstsq(X, y[:, ~mask], rcond=None)  # 只对有效格点回归
3

4
# 创建全 NaN 数组，回填计算结果
5
reg_map = np.zeros((2, y.shape[1])) * np.nan
6
reg_map[:, ~mask] = r

1
# 将不显著的格点设为 NaN（绘图时不显示）
2
sst_tmp = np.where(map_sst.p_value >= 0.05, np.nan, map_sst.slope)
3
uwnd_tmp = np.where(map_uwnd.p_value >= 0.10, np.nan, map_uwnd.slope)

1
# xarray → numpy
2
sst_np = sst.to_numpy()
3

4
# numpy → xarray（保留原维度坐标）
5
sst_xr = xr.DataArray(map, dims=sst.dims, coords=sst.coords)

1
from eofs.standard import Eof
2

3
solver = Eof(ssta, weights=wgts)                        # 初始化
4
pt = solver.eofsAsCovariance(neofs=4)                    # 空间模态（协方差）
5
pc = solver.pcs(npcs=4, pcscaling=1)                     # 时间序列 PC
6
pcvar = solver.varianceFraction(neigs=4) * 100           # 解释方差百分比(%)

1
import sacpy as scp
2

3
eof = scp.EOF(ssta)
4
eof.solve()
5
pc = -eof.get_pc(npt=3)                 # PC（取负号可翻转符号，因人而异）
6
pt = -eof.get_pt(npt=3)                 # EOF 空间模态
7
pcvar = eof.get_varperc(npt=3) * 100    # 解释方差百分比(%)

1
dat = dat - 273.15    # 开尔文 → 摄氏度

1
wind_speed = np.sqrt(uwnd**2 + vwnd**2)

1
import sacpy as scp
2

3
ssta = scp.get_anom(sst, method=1)

1
LinReg = scp.LinReg(nino34, sst)
2
slope = LinReg.slope          # 回归系数场 (lat, lon)
3
p_value = LinReg.p_value      # p 值场 (lat, lon)

1
eof = scp.EOF(ssta)                  # 初始化
2
eof.solve()                          # 求解
3
pc = eof.get_pc(npt=3)               # 获取 PC (npt 个模态)
4
pt = eof.get_pt(npt=3)               # 获取特征向量（空间模态）
5
pcvar = eof.get_varperc(npt=3) * 100 # 解释方差百分比(%)