Python | 超前滞后分析

需求幻想家
• 阅读 182

Nino SST Indices (Nino 1+2, 3, 3.4, 4; ONI and TNI)

有几个指标用于监测热带太平洋,所有这些指标都是基于海表温度(SST)异常在一个给定的区域的平均值。通常,异常是相对于30年的周期来计算的。厄尔尼诺3.4指数(Niño 3.4 index)和海洋厄尔尼诺指数(Oceanic Niño Index (ONI))是定义厄尔尼诺和拉尼娜事件最常用的指数。其他index用于帮助描述每个事件的独特性质。

Python | 超前滞后分析

El Niño 1、2、3和4 的范围如下所示:

Python | 超前滞后分析

Niño 1+2 (0-10S,90W-80W) :

  • 厄尔尼诺1 + 2区域是Niño海表温度区域中最小和最东部的区域,与南美洲沿海地区相对应,在那里El Niño,首先被当地居民识别发现。这一指数往往是Niño海表温度指数中变化最大的。

Niño 3(5N-5S,150W-90W) :

  • 这个区域曾经是监测和预测El Niño现象的主要焦点,但是研究人员后来了解到,ENSO 海洋-大气耦合相互作用的关键区域位于更西部(Trenberth,1997)。因此,Niño 3.4和 ONI 成为定义厄尔尼诺和拉尼娜事件的首选
    Niño 3.4(5N-5S,170W-120W) :
  • Niño3.4异常可以被认为代表从日期线到南美海岸的横跨太平洋的赤道 SST 的平均值。Niño 3.4 index通常使用5个月的连续平均值,El Niño or La Niña是在Niño 3.4 SST 超过 +/-0.4 C 6个月或更长时间内确定的。

ONI (5N-5S,170W-120W) :

  • ONI 使用与厄尔尼诺3.4指数相同的区域。ONI 采用的是连续3个月的平均值,要被归类为全面的厄尔尼诺或拉尼娜现象,异常必须至少连续5个月超过 + 0.5 C 或 -0.5 C。这是美国国家海洋和大气管理局使用的定义。

Niño 4(5N-5S,160E-150W) :

  • 厄尔尼诺4指数捕捉赤道太平洋中部的海表温度异常。这个区域的变化往往比其他厄尔尼诺区域少。

为了定义每个 El Niño or La Niña 事件的独特性,Trenberth 和 Stepaniak (2001)认为,应该将厄尔尼诺3.4指数与他们引入的一个指数结合使用,这个指数被称为跨尼诺指数(Trans-Niño Index,TNI)。TNI 被定义为厄尔尼诺1 + 2和厄尔尼诺4区域之间标准化海表温度异常的差异。TNI 因此测量了赤道太平洋中部和东部海表温度异常的梯度。当海表温度梯度特别大时(比如,由于尼诺4区域的正异常和尼诺1 + 2区域的负异常) ,一些研究人员将该事件归类为“太平洋中部的厄尔尼诺现象”或“El Niño Modoki”,尽管将这种类型的事件作为一个单独的现象进行区分是一个问题或争论。

Nino 3.4 index 的具体计算步骤

  • 1、选择厄尔尼诺3.4区域(5N-5S,170W-120W)的sst月平均数据
  • 2、计算Nino 3.4区域内每个月的异常
  • 3、计算每个月异常的区域平均值
  • 4、去线性趋势化处理(也可以不去)

代码实现过程:


# == netcdf file name and location"
fnc = 'oisst_monthly.nc'
dmask = xr.open_dataset('lsmask.nc')
print(dmask)

ds = xr.open_dataset(fnc)
print(ds)

# === Climatology and Anomalies
sst = ds.sst.where(dmask.mask.isel(time=0) == 1)
clm = sst.sel(time=slice('1982-01-01','2020-12-01')).groupby('time.month').mean(dim='time')
anm = (sst.groupby('time.month') - clm)
#print(clm)
# -- Detorending
def detrend_dim(da, dim, deg=1):
    # detrend along a single dimension
    p = da.polyfit(dim=dim, deg=deg)
    fit = xr.polyval(da[dim], p.polyfit_coefficients)
    return da - fit

def cal_nino34(ds,lat1,lat2,lon1,lon2):
    ds = ds.sel(lat=slice(lat1,lat2),lon=slice(lon1,lon2))
    clm = ds.groupby('time.month').mean(dim='time')
    an_ds = ds.groupby('time.month') - clm
    an_ds = an_ds.mean(('lon','lat'))
    an_ds = detrend_dim(an_ds,dim='time',deg=1)
    return an_ds


nino3_4 = cal_nino34(sst,5,-5,190,240)
nino3_4.plot()

结果展示:

Python | 超前滞后分析

平滑

还可以对于指数进行平滑,使结果看起来更好看一点。如果指数平滑了,后续计算相关相关时也需要对于输入的其他变量进行统一平滑。


ninoSD=nino3_4/nino3_4.std(dim='time')
rninoSD=ninoSD.rolling(time=3, center=True).mean('time')

Python | 超前滞后分析

计算超前滞后

相关性和回归分析部分:

Leading:

  • 通过计算NINO指数(rninoSD)和时间滞后12个时间步长的rdanm之间的相关系数和回归系数。类似地,也计算了时间滞后6个时间步长的相关系数和回归系数。
    Simultaneous:
  • 计算NINO指数(rninoSD)和rdanm之间同期的相关系数和回归系数。
    Lagging:
  • 计算NINO指数(rninoSD)和时间滞后6、12、18个时间步长的rdanm之间的相关系数和回归系数。
# -- Running mean

ranm = anm.rolling(time=3, center=True).mean('time')
rdanm = detrend_dim(ranm,'time',1)

# -- Correlation & Regression


# Leading
corM12 = xr.corr(rninoSD, rdanm.shift(time=-12), dim="time")
regM12 = xr.cov( rninoSD, rdanm.shift(time=-12), dim="time")/rninoSD.var(dim='time',skipna=True).values
corM6 = xr.corr(rninoSD, rdanm.shift(time=-6), dim="time")
regM6 = xr.cov( rninoSD, rdanm.shift(time=-6), dim="time")/rninoSD.var(dim='time',skipna=True).values

# simultaneous
cor0 = xr.corr(rninoSD, rdanm, dim="time")
reg0 = xr.cov(rninoSD, rdanm, dim="time")/rninoSD.var(dim='time',skipna=True).values

# Laging
corP6 = xr.corr(rninoSD, rdanm.shift(time=6), dim="time")
regP6 = xr.cov( rninoSD, rdanm.shift(time=6), dim="time")/rninoSD.var(dim='time',skipna=True).values
corP12 = xr.corr(rninoSD, rdanm.shift(time=12), dim="time")
regP12 = xr.cov( rninoSD, rdanm.shift(time=12), dim="time")/rninoSD.var(dim='time',skipna=True).values
corP18 = xr.corr(rninoSD, rdanm.shift(time=18), dim="time")
regP18 = xr.cov( rninoSD, rdanm.shift(time=18), dim="time")/rninoSD.var(dim='time',skipna=True).values

绘图

# -- figure plot

def makefig(cor, reg,title, grid_space):
  # 修复 0 度和 360 度经度附近未显示数据的伪影
  cor = gvutil.xr_add_cyclic_longitudes(cor, 'lon')
  reg = gvutil.xr_add_cyclic_longitudes(reg, 'lon')
  # 添加等距柱面投影,中心经度为210°
  ax = fig.add_subplot(grid_space,
          projection=ccrs.PlateCarree(central_longitude=210))
  # 添加海岸线
  ax.coastlines(linewidth=0.5, alpha=0.6)
  
  # 设置坐标范围
  gvutil.set_axes_limits_and_ticks(ax,
                                 xlim=(-180, 180),
                                 ylim=(-90, 90),
                                 xticks=np.arange(-180, 180, 60),
                                 yticks=np.arange(-90, 90, 30))
  
  # Use geocat.viz.util convenience function to add minor and major tick lines
  gvutil.add_major_minor_ticks(ax, labelsize=10)
  
  # Use geocat.viz.util convenience function to make latitude, longitude tick labels
  gvutil.add_lat_lon_ticklabels(ax)
  
  #设置colorbar
  newcmp = cmaps.NCV_blu_red
  index = [5, 20,  35, 50, 65, 85, 95, 110, 125,  0, 0, 135, 150,  165, 180, 200, 210, 220, 235, 250 ]
  color_list = [newcmp[i].colors for i in index]
  # 设置colorbar中间颜色为白色
  color_list[9]=[ 1., 1., 1.]
  color_list[10]=[ 1., 1., 1.]
  
  # 定义填色图的参数
  kwargs = dict(
    vmin = -1.0,
    vmax = 1.0,
    levels = 21,
    colors=color_list,
    add_colorbar=False,  # allow for colorbar specification later
    transform=ccrs.PlateCarree(),  # ds projection
  )
  
  # 相关系数的填色图
  fillplot = cor.plot.contourf(ax=ax,  **kwargs)

  # 添加陆地、地形
  ax.add_feature(cfeature.LAND, facecolor='lightgray', zorder=1)
  ax.add_feature(cfeature.COASTLINE, edgecolor='gray', linewidth=0.5, zorder=1)
  
  # 设置等值线参数
  # Specify contour levels excluding 0
  delc=0.2
  levels = np.arange(-3, 0, delc)
  levels = np.append(levels, np.arange(delc, 3, delc))
  # 回归系数为等值线
  rad = reg.plot.contour(ax=ax,
                  colors='black',
                  alpha=0.8,
                  linewidths=1.0,
                  add_labels=False,
                  levels=levels,
                  transform=ccrs.PlateCarree())

  pe = [PathEffects.withStroke(linewidth=2.0, foreground="w")]
  plt.setp(rad.collections, path_effects=pe)
  # 设置标题、及其位置、大小
  gvutil.set_titles_and_labels(ax,
                             lefttitle=title,
                             lefttitlefontsize=16,
                             righttitle='',
                             righttitlefontsize=16,
                             xlabel="",
                             ylabel="")

  return ax, fillplot

# Show the plot

fig = plt.figure(figsize=(10, 12),dpi=200)
grid = fig.add_gridspec(ncols=2, nrows=3)
#grid = fig.add_gridspec(ncols=2, nrows=3, hspace=-0.20)


ax1, fill1 = makefig(corP18,regP18,'18-month lag', grid[0,0])
ax2, fill2 = makefig(corP12,regP12,'12-month lag', grid[1,0])
ax3, fill3 = makefig(corP6,regP6,'6-month lag', grid[2,0])
ax4, fill4 = makefig(cor0,reg0,'Simultaneous', grid[0,1])
ax5, fill5 = makefig(corM6,regM6,'6-month lead', grid[1,1])
ax6, fill6 = makefig(corM12,regM12,'12-month lead', grid[2,1])

fig.colorbar(fill6,
                 ax=[ax1, ax2, ax3, ax4, ax5, ax6],
                 drawedges=True,
                 orientation='horizontal',
                 shrink=0.5,
                 pad=0.05,
                 extendfrac='auto',
                 extendrect=True)
  
fig.suptitle('SST correlation & regression with Nino3.4', fontsize=18, y=0.9)

plt.draw()

Python | 超前滞后分析

https://climatedataguide.ucar.edu/climate-data/nino-sst-indic...

Trenberth, Kevin & National Center for Atmospheric Research Staff (Eds). Last modified 2024-03-20 "The Climate Data Guide: Nino SST Indices (Nino 1+2, 3, 3.4, 4; ONI and TNI).”

Trenberth, K.E. and Stepaniak, D.P. (2001) Indices of El Nino Evolution. Journal of Climate, 14, 1697-1701.
https://doi.org/10.1175/1520-0442(2001)014<1697:LIOENO>2.0.CO;2

https://climate.usu.edu/people/yoshi/pyclm101/index.html

本文由mdnice多平台发布

点赞
收藏
评论区
推荐文章
Wesley13 Wesley13
3年前
MySQL部分从库上面因为大量的临时表tmp_table造成慢查询
背景描述Time:20190124T00:08:14.70572408:00User@Host:@Id:Schema:sentrymetaLast_errno:0Killed:0Query_time:0.315758Lock_
美凌格栋栋酱 美凌格栋栋酱
6个月前
Oracle 分组与拼接字符串同时使用
SELECTT.,ROWNUMIDFROM(SELECTT.EMPLID,T.NAME,T.BU,T.REALDEPART,T.FORMATDATE,SUM(T.S0)S0,MAX(UPDATETIME)CREATETIME,LISTAGG(TOCHAR(
Prodan Labs Prodan Labs
4年前
Kubernetes自定义调度器 — 初窥门径
通过上一篇文章对schedulerframework调度框架已经有了大致了解,根据我们的实际生产的一些问题(如计算服务没有被调度到实际CPU最优的节点)和需求,来实现一个简单的基于CPU指标的自定义调度器。自定义调度器通过kubernetes资源指标服务metricsserver来获取各节点的当前的资源情况,并进行打分,然后把Pod调度到分数最高的节
Aidan075 Aidan075
4年前
数据分析面试-业务分析篇
一.指标建模1.1互联网产品常用指标1.2选好数据指标的通用方法论三部曲1.从业务的最终目的出发梳理业务模块拆解:多问几个how目的:我要卖货手段:通过图文来买货支撑手段的手段:通过社区创作的优美的图文来卖货2.判断业务模块所属类型业务类型二.数据分析方法论2.1对比分析绝对值(本身具备价值的数字)eg:销售金融阅读数比例值(在具体环境中看
Karen110 Karen110
4年前
人工智能数学基础4:离差、平均差、方差、标准差、协方差、皮尔森相关系数
一、离差(Deviation)离差即标志变动度,又称“偏差”,是观测值或估计量的平均值与真实值之间的差,是反映数据分布离散程度的量度之一,或说是反映统计总体中各单位标志值差别大小的程度或离差情况的指标,常写作:即参与计算平均数的变量值与平均数之差。离差的性质有二:(1)离差的代数和等于0;(2)参与计算平均数的各变量值与平均数之差的平均和,小于这些变量值
Wesley13 Wesley13
3年前
PIC中档单片机汇编指令详解(5)
位操作指令详述BCF数据寄存器指定位清0语法形式:BCFf,b操作数:f为数据寄存器的低7位地址(0x00~0x7F)B为数据位编号(0~7)执行时间:一个指令周期执行过程:使数据寄存器f的的b位清0状态标志影响:无说明:该指令可对任何数据寄存器的任意一个位置清0,常用于标志位的设定和清除,或者把某一管脚置成低电平。指
Stella981 Stella981
3年前
Prometheus监控学习笔记之PromQL简单示例
0x00简单的时间序列选择返回度量指标http_requests_total的所有时间序列样本数据:http_requests_total返回度量指标名称为http_requests_total,标签分别是job"apiserver",handler"/api/comments"
Stella981 Stella981
3年前
Python之time模块的时间戳、时间字符串格式化与转换
Python处理时间和时间戳的内置模块就有time,和datetime两个,本文先说time模块。关于时间戳的几个概念时间戳,根据1970年1月1日00:00:00开始按秒计算的偏移量。时间元组(struct_time),包含9个元素。 time.struct_time(tm_y
Wesley13 Wesley13
3年前
mysql设置时区
mysql设置时区mysql\_query("SETtime\_zone'8:00'")ordie('时区设置失败,请联系管理员!');中国在东8区所以加8方法二:selectcount(user\_id)asdevice,CONVERT\_TZ(FROM\_UNIXTIME(reg\_time),'08:00','0
Wesley13 Wesley13
3年前
2019年北航OO第1单元(表达式求导)总结
2019年北航OO第1单元(表达式求导)总结1基于度量的程序结构分析量化指标及分析以下是三次作业的量化指标统计:!(https://img2018.cnblogs.com/blog/1616496/201903/1616496201903251247313
Wesley13 Wesley13
3年前
Java 应用发布后,需要关注的7个性能指标
在某个重大发布之后,都需要记录相应的指标,本文介绍了最重要的几个Java性能指标,包括响应时间和平均负载等。为理解应用程序在生产环境中如何运行,就需要遵循一些Java性能指标。在以前,当软件被发布后,开发者是没有方法去了解它在生产环境中的运行情况;而现在,几乎任一个你可以想到的指标都可以被监测和报告。时下,开发者面临的问题并不是缺乏信息,