python计算nino3.4指数

本文定义厄尔尼诺与拉尼娜事件的方法为美国气象中心（CPC）的定义方法，即以30年为一次气候模态，三个月滑动平均的Niño 3.4指数连续五个月大于0.5℃被认为是一次厄尔尼诺事件，连续五个月小于-0.5℃则被认为是一次拉尼娜事件，Niño 3.4区的范围为5°N–5°S, 120°–170°W。

 我们分析的是海表面温度(SST)，我选取的是1951-1980年的数据，全球海洋温度月数据可以从美国国家气象海洋管理局（NOAA）上下载，格式是科学文件格式.nc。

首先我们引用可能会使用的第三方库，第三方库的下载可参考这篇文章（简短干脆）如何在windows系统中安装anaconda的教程_conda install cdsapi-CSDN博客

#第三方库
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import xarray as xr
from cartopy.mpl import ticker
import netCDF4 as nc
import pandas as pd

这些第三方库基本是使用python读取nc文件画图使用的必需，可以提前引用。

我们读取自己的sst.nc文件，

#读取文件
path='C:\\sst.nc'
SST=xr.open_dataset(path)

查看我们的SST信息

可以看出我的数据是从1854年1月-2023年2月的全球数据，但是我们需要分析的是1951-1980年的数据（其实1950年之前的数据准确度都不太够，所以我分析的时候一般选取1950作为起点分析）

其中我们也可以发现我们的经度范围为0-360，不需要使用cdo对经度进行处理，如果大家的数据需要处理可以参考如何解决 cdo转换经度-180~180 为0~360时出现报错：cdo sellonlatbox: Unsupported grid type: generic的问题_经度360和-180-CSDN博客我们此处的分析范围是太平洋地区的厄尔尼诺和拉尼娜现象，所以还需要我们对经纬度进行选取

path='C:\\sst.nc'
SST=xr.open_dataset(path).sel(lat=slice(6,-6),lon=slice(190,240),time=slice("1950","1981"))
sst1=SST.sst[:]
sst2=np.array(sst1)
lat=SST.lat[:]
lon=SST.lon[:]
#将数据转化为array格式
sst=np.array(sst1)
lat=np.array(lat)
lon=np.array(lon)

此时我们的温度和经纬度数据都转化为普通的array数组数据，但是sst中存在nan(缺测值)会影响后续的计算或画图，所以要进行处理

#将nan值转换为0
sst=np.nan_to_num(sst)

接下来就是选取30年气候态，也就是1951-1980年的数据，由于我的sst数据维度分辨率是2°，我需要的是5°S-5°N，所以要对-6°-6°和-4°-4°做一次平均

#气候态，温度异常
a = np.zeros((360,7,26))
for j in range (12):
    for i in range (0,360,12):
        a[j,:,:]=a[j,:,:]+sst[i+j-1,:,:]
    a[j,:,:]=a[j,:,:]/30


sst_ano = np.zeros((360,7,26)) 
for j in range (12):
    for i in range (0,360,12):
        sst_ano[i+j-1,:,:]= sst[i+j-1,:,:]-a[j,:,:]

#{6:-6}
b = np.zeros(360)
for n  in range (360):
    for i in range (26):
        for j in range(7):
            b[n] = b[n]+sst_ano[n,j,i]
    b[n] = b[n]/182   

#{4:-4}
c = np.zeros(360)
for n in range(360):
    for i in range (1,26):
        for j in range (7):
            c[n] = c[n]+sst_ano[n,j,i]
    c[n] = c[n]/175



#去除线性趋势后算数平均
n1 = np.nan_to_num(b - np.nanmean(b))
n2 = np.nan_to_num(c - np.nanmean(c))
nino = (n1 + n2) / 2

 此时得到的nino还是比较粗糙的数据，最后要进行滑动3个月平均，python中有对应的函数，直接引用就可以

#滑动3个月平均
a=pd.Series(nino)
nino34 = a.rolling(3,min_periods=1).mean()

这样，我们想要的nino3.4指数就算出啦

完整代码如下

#第三方库
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import xarray as xr
from cartopy.mpl import ticker
import netCDF4 as nc
import pandas as pd

#数据的路径
path='C:\\sst.nc'
SST=xr.open_dataset(path).sel(lat=slice(6,-6),lon=slice(190,240),time=slice("1951","1980"))
sst1=SST.sst[:]
sst2=np.array(sst1)
lat=SST.lat[:]
lon=SST.lon[:]
#将数据转化为array格式
sst=np.array(sst1)
lat=np.array(lat)
lon=np.array(lon)

#将nan值转换为0
sst=np.nan_to_num(sst)

#气候态，温度异常
a = np.zeros((360,7,26))
for j in range (12):
    for i in range (0,360,12):
        a[j,:,:]=a[j,:,:]+sst[i+j-1,:,:]
    a[j,:,:]=a[j,:,:]/30


sst_ano = np.zeros((360,7,26)) 
for j in range (12):
    for i in range (0,360,12):
        sst_ano[i+j-1,:,:]= sst[i+j-1,:,:]-a[j,:,:]

#{6:-6}
b = np.zeros(360)
for n  in range (360):
    for i in range (26):
        for j in range(7):
            b[n] = b[n]+sst_ano[n,j,i]
    b[n] = b[n]/182   

#{4:-4}
c = np.zeros(360)
for n in range(360):
    for i in range (1,26):
        for j in range (7):
            c[n] = c[n]+sst_ano[n,j,i]
    c[n] = c[n]/175

#去除线性趋势后算数平均
n1 = np.nan_to_num(b - np.nanmean(b))
n2 = np.nan_to_num(c - np.nanmean(c))
nino = (n1 + n2) / 2

#滑动3个月平均
a=pd.Series(nino)
nino34 = a.rolling(3,min_periods=1).mean()