Tropical Waves Identification and Analysis

Tropical Waves Identification and Analysis#

本單元主要示範如何利用Python工具判識熱帶波動訊號，而非提供熱帶波動理論完整的回顧與介紹。完整的理論介紹請參照下方所列的參考文獻。

熱帶波動理論簡介#

熱帶波動是主導熱帶地區季內尺度變異的系統，它的理論最早在 Matsuno (1966) 文章中提出，這篇文章從熱帶 \(\beta \) 面上的淺水方程 (shallow water equations) 出發，導出了熱帶波動頻散關係的理論解：

\[ \omega ^{2} - k^2 - \dfrac{k}{\omega } = 2n+1, \qquad n = 1, 2, ... \]

若將頻散關係方程式繪圖，可以得到如下圖的結果。

https://agupubs.onlinelibrary.wiley.com/cms/asset/15874aea-bc1f-4385-b84c-26698dd2748d/rog1687-fig-0002.png — Fig. 1 Dispersion curves for equatorial waves (up to \(n = 4\)) as a function of the nondimensional frequency, \(\omega ^{*}\), and nondimensional zonal wave number, \(k^{*}\), where \(\omega ^{*} \equiv \left( \omega /\beta \sqrt{gh_e} \right) ^{1/2}\), and \(k^{*} \equiv k\left( \sqrt{gh_e} / \beta \right) ^{1/2}\). For all but the Kelvin wave, these dispersion curves are solutions of equation (1). Westward propagating waves (relative to the zero basic state) appear on the left, and eastward propagating waves appear on the right. For consistency with equation (1), the Kelvin wave solution is labeled as \(n = −1\). (Source: Kiladis et al. (2009))#

在這張圖中，\(k^{*} > 0\) 的部分代表往東傳遞的波動，\(k^{*} < 0\) 則代表往西傳遞的波動。隨後，許多觀測實驗證實了熱帶確實有這些往東和往西的波動存在，以及當衛星觀測和再分析資料越來越盛行後，從波譜分析的結果也可以發現波譜的峰值區域和理論解非常吻合，甚至可以根據波譜分析的結果進行適當的濾波來得到熱帶波動的資料。而從淺水方程式導出的波動的動力結構如下圖所示：

https://agupubs.onlinelibrary.wiley.com/cms/asset/e5386eb7-680f-46b0-9887-1f9389a70417/rog1687-fig-0003.png — Fig. 2 Horizontal structures of a subset of the zonally propagating wave solutions to the shallow water equations on an equatorial \(\beta \) plane (equations (1)–(3)). Each is shown for a nondimensional zonal wave number, \(k^{*} = \pm 1\). All scales and fields have been nondimensionalized by taking the units of time and length as given in the caption of Figure 1. The equator runs through the center of each diagram. Hatching is for divergence, and shading is for convergence, with a 0.6 unit interval between successive levels. Unshaded contours are geopotential, with a contour interval of 0.5 units. Negative contours are dashed, and the zero contour is omitted. The maximum wind vectors in each panel are specified in the bottom right corner. (Source: Kiladis et al. (2009))#

以下會介紹兩種常用的分析方法，其一是 Wheeler and Kiladis (1999) 提出的時空濾波法 (space-time filtering technique)，其二是 Yang et al. (2003) 提出的雙曲柱狀函數 (PCF) 投影法 (2D spatial projection on the parabolic cylinder functions, 2DS-PCF) 方法。其他還有許多熱帶波動辨識的方法，詳細的介紹與回顧請見 Knippertz et al. (2022)。

雙曲柱狀函數投影法 (2DS-PCF)#

在Yang et al. (2003) 中，他們指出由於熱帶波動會受到背景場調節，產生都卜勒效應，波動的頻散關係會被改變，因此像Wheeler and Kiladis (1999) 的方法去限制一個窄的時空頻段來進行濾波，可能會造成有部分的波動訊號被忽略。因此，他們只對動力場 (u, v, z) 做一個寬帶濾波，波數的範圍是 \(1\leq \left| k\right| \leq 15\)，週期的範圍是 \(2\leq \left| T\right| \leq 30\) 天，再將寬帶濾波後的動力場分別投影至各個波動模態的理論解上，此時得到的就是熱帶波動，其中理論解釋是根據 Gill (1980) ，這篇文章假定表面有加熱的作用，加入淺水方程中，導出的解就是雙曲柱狀函數的形式。

https://rmets.onlinelibrary.wiley.com/cms/asset/0cf13e1c-5d09-4004-9f66-8dc3b634d610/qj4338-fig-0002-m.jpg — Fig. 4 Latitudinal profiles of different equatorial wave solutions. Spherical solutions for (a) Kelvin wave zonal wind \(u^{\prime }\) and (b) mixed Rossby–gravity wave (MRG) meridional wind \(v^{\prime }\), both for \(D\) = 8 m and 90 m and for \(k\) = 1 and \(k\) = 20, as well as (c) the \(n\) = 1 equatorial Rossby wave (ER) \(u^{\prime }, v^{\prime }\), and geopotential \(\phi ^{\prime }\) for \(k\) = 1 and \(D\) = 90 m. These curves were produced using the software developed by Swarztrauber and Kasahara (1985). (d) Parabolic cylinder functions (PCFs) 0, 1, and 2 for \(y_0 = 6\). The 15° S–15° N belt used for averaging throughout the paper is shaded in grey. Source: Knippertz et al. 2023.#

Note

這個方法必須同時使用u, v, z三個變數！
這個方法不能濾MJO的訊號！

2DS-PCF方法的程式，是由Professor Steven Woolnough 和 Dr. Gui-Ying Yang (University of Reading) 提供，並在 Dr. Sandro Lubis的GitHub網頁下載 (sandrolubis/CCEWs-PCF-Filter)。以下示範對2016-2020年的u, v, z資料進行熱帶波動辨識 (此處示範使用Dr. Lubis提供的範例資料檔)。

'''
; Authors: Prof. Steven Woolnough and Dr. Gui-Ying Yang
;          Modified by Dr. Sandro W. Lubis (Nov 2021)
;          CCEW Filter via PCFs following Yang et al., (2003)
; Contact: slubis.geomar@gmail.com
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
'''
import numpy as np
import pandas as pd
import xarray as xr
import pandas as pd
from scipy import signal
from metpy.units import units

### Read data and calculate anomaly

lev = 850

# Open files 
ua = xr.open_dataarray('./data/u.850.20N-20S.0-360.20160101-20201231.nc')
va = xr.open_dataarray('./data/v.850.20N-20S.0-360.20160101-20201231.nc')
za = xr.open_dataarray('./data/phi.850.20N-20S.0-360.20160101-20201231.nc')

### 2DS-PCF starts (以下都不需要更動)
u = signal.detrend(ua.values, axis=0)
v = signal.detrend(va.values, axis=0)
z = signal.detrend(za.values, axis=0)

wgt_taper = signal.tukey(u.shape[0], alpha=0.1)

uT = np.transpose(u, (1, 2, 0)) * np.reshape(wgt_taper, (-1, u.shape[0]))
vT = np.transpose(v, (1, 2, 0)) * np.reshape(wgt_taper, (-1, v.shape[0]))
zT = np.transpose(z, (1, 2, 0)) * np.reshape(wgt_taper, (-1, z.shape[0]))

u = np.transpose(uT, (2, 0, 1))
v = np.transpose(vT, (2, 0, 1))
z = np.transpose(zT, (2, 0, 1))

g     = 9.8
beta  = 2.3e-11
radea = 6.371e+06
spd   = 86400.0
ww    = 2.0 * np.pi / spd

latmax = 20.0

kmin = 2
kmax = 20

pmin = 3.0
pmax = 30.0

# convert trapping scale to meters
y0 = 6.0
y0real = 2.0 * np.pi * radea * y0 / 360.0

ce = 2.0 * y0real**2 * beta

g_on_c = g / ce
c_on_g = ce / g

waves = np.array(['Kelvin', 'WMRG', 'R1', 'R2'])

# transform u,z to q, r using q=z*(g/c) + u; r=z*(g/c) - u 

q = z * g_on_c + u
r = z * g_on_c - u

qf = np.fft.fft2(q, axes=(0, 2))
vf = np.fft.fft2(v, axes=(0, 2))
rf = np.fft.fft2(r, axes=(0, 2))

nf = qf.shape[0]
nlat = qf.shape[1]
nk = qf.shape[2]

# Find frequencies and wavenumbers corresponding to pmin,pmax and kmin,kmax in coeff matrices

f = np.fft.fftfreq(nf)
k = np.fft.fftfreq(nk) * nk

fmin = np.where(f >= 1.0 / pmax)[0][0]
fmax = (np.where(f > 1.0 / pmin)[0][0]) - 1

f1p = fmin
f2p = fmax + 1
f1n = nf - fmax
f2n = nf - fmin + 1

k1p = kmin
k2p = kmax + 1
k1n = nk - kmax
k2n = nk - kmin + 1

# Define the parobolic cylinder functions

spi2 = np.sqrt(2.0 * np.pi)
dsq = np.array([spi2, spi2, 2.0 * spi2, 6.0 * spi2])
d = np.zeros((dsq.size, nlat))
y = ua.Y.values / y0
ysq = y**2

d[0, :] = np.exp(-ysq / 4.0)
d[1, :] = y * d[0, :]
d[2, :] = (ysq - 1.0) * d[0, :]
d[3, :] = y * (ysq - 3.0) * d[0, :]

dlat = np.abs(ua.Y.values[1] - ua.Y.values[0]) * np.pi / 180.0

qf_Kel = np.zeros((nf, nk), dtype='complex')

qf_mode = np.zeros((dsq.size, nf, nk), dtype='complex')
vf_mode = np.zeros((dsq.size, nf, nk), dtype='complex')
rf_mode = np.zeros((dsq.size, nf, nk), dtype='complex')

# reorder the spectral coefficents to make the latitudes the last dimension

qf = np.transpose(qf, (0, 2, 1))
vf = np.transpose(vf, (0, 2, 1))
rf = np.transpose(rf, (0, 2, 1))

for m in np.arange(dsq.size):
    if m == 0:
        qf_Kel[f1n:f2n, k1p:k2p] = np.sum(qf[f1n:f2n, k1p:k2p, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)
        qf_Kel[f1p:f2p, k1n:k2n] = np.sum(qf[f1p:f2p, k1n:k2n, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)
    
    qf_mode[m, f1n:f2n, k1n:k2n] = np.sum(qf[f1n:f2n, k1n:k2n, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)
    qf_mode[m, f1p:f2p, k1p:k2p] = np.sum(qf[f1p:f2p, k1p:k2p, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)
    
    vf_mode[m, f1n:f2n, k1n:k2n] = np.sum(vf[f1n:f2n, k1n:k2n, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)
    vf_mode[m, f1p:f2p, k1p:k2p] = np.sum(vf[f1p:f2p, k1p:k2p, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)
    
    rf_mode[m, f1n:f2n, k1n:k2n] = np.sum(rf[f1n:f2n, k1n:k2n, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)
    rf_mode[m, f1p:f2p, k1p:k2p] = np.sum(rf[f1p:f2p, k1p:k2p, :] * d[m, :] * dlat, axis=-1) / (dsq[m] / y0)

uf_wave = np.zeros((waves.size, nf, nlat, nk), dtype='complex')
vf_wave = np.zeros((waves.size, nf, nlat, nk), dtype='complex')
zf_wave = np.zeros((waves.size, nf, nlat, nk), dtype='complex')

for w in np.arange(waves.size):
    if waves[w] == 'Kelvin':
        for j in np.arange(nlat):
            uf_wave[w, :, j, :] = 0.5 * qf_Kel * d[0, j]
            zf_wave[w, :, j, :] = 0.5 * qf_Kel * d[0, j] * c_on_g
    
    if waves[w] == 'WMRG':
        for j in np.arange(nlat):
            uf_wave[w, :, j, :] = 0.5 * qf_mode[1, :, :] * d[1, j]
            vf_wave[w, :, j, :] = 0.5 * vf_mode[0, :, :] * d[0, j]
            zf_wave[w, :, j, :] = 0.5 * qf_mode[1, :, :] * d[1, j] * c_on_g
    
    if waves[w] == 'R1':
        for j in np.arange(nlat):
            uf_wave[w, :, j, :] = 0.5 * (qf_mode[2, :, :] * d[2, j] - rf_mode[0, :, :] * d[0, j])
            vf_wave[w, :, j, :] = 0.5 * vf_mode[1, :, :] * d[1, j]
            zf_wave[w, :, j, :] = 0.5 * (qf_mode[2, :, :] * d[2, j] + rf_mode[0, :, :] * d[0, j]) * c_on_g
    
    if waves[w] == 'R2':
        for j in np.arange(nlat):
            uf_wave[w, :, j, :] = 0.5 * (qf_mode[3, :, :] * d[3, j] - rf_mode[1, :, :] * d[1, j])
            vf_wave[w, :, j, :] = 0.5 * vf_mode[2, :, :] * d[2, j]
            zf_wave[w, :, j, :] = 0.5 * (qf_mode[3, :, :] * d[3, j] + rf_mode[1, :, :] * d[1, j]) * c_on_g

u_Kelvin = xr.DataArray(np.real(np.fft.ifft2(uf_wave[0, :, :, :], axes=(0, 2))), coords=[ua['T'], ua['Y'], ua['X']], dims=['T', 'Y', 'X'], name='u')
v_Kelvin = xr.DataArray(np.real(np.fft.ifft2(vf_wave[0, :, :, :], axes=(0, 2))), coords=[va['T'], va['Y'], va['X']], dims=['T', 'Y', 'X'], name='v')
z_Kelvin = xr.DataArray(np.real(np.fft.ifft2(zf_wave[0, :, :, :], axes=(0, 2))), coords=[za['T'], za['Y'], za['X']], dims=['T', 'Y', 'X'], name='z')

u_Kelvin.attrs['long_name'] = 'Kelvin Waves in '+ str(lev) + ' hPa Zonal Wind'
v_Kelvin.attrs['long_name'] = 'Kelvin Waves in '+ str(lev) + ' hPa Meridional Wind'
z_Kelvin.attrs['long_name'] = 'Kelvin Waves in '+ str(lev) + ' hPa Geopotential Height'

u_Kelvin.attrs['units'] = 'm/s'
v_Kelvin.attrs['units'] = 'm/s'
z_Kelvin.attrs['units'] = 'm'

u_WMRG = xr.DataArray(np.real(np.fft.ifft2(uf_wave[1, :, :, :], axes=(0, 2))), coords=[ua['T'], ua['Y'], ua['X']], dims=['T', 'Y', 'X'], name='u')
v_WMRG = xr.DataArray(np.real(np.fft.ifft2(vf_wave[1, :, :, :], axes=(0, 2))), coords=[va['T'], va['Y'], va['X']], dims=['T', 'Y', 'X'], name='v')
z_WMRG = xr.DataArray(np.real(np.fft.ifft2(zf_wave[1, :, :, :], axes=(0, 2))), coords=[za['T'], za['Y'], za['X']], dims=['T', 'Y', 'X'], name='z')

u_WMRG.attrs['long_name'] = 'Westward Mixed Rossby-Gravity Waves in '+ str(lev) + ' hPa Zonal Wind'
v_WMRG.attrs['long_name'] = 'Westward Mixed Rossby-Gravity Waves in '+ str(lev) + ' hPa Meridional Wind'
z_WMRG.attrs['long_name'] = 'Westward Mixed Rossby-Gravity Waves in '+ str(lev) + ' hPa Geopotential Height'

u_WMRG.attrs['units'] = 'm/s'
v_WMRG.attrs['units'] = 'm/s'
z_WMRG.attrs['units'] = 'm'

u_R1 = xr.DataArray(np.real(np.fft.ifft2(uf_wave[2, :, :, :], axes=(0, 2))), coords=[ua['T'], ua['Y'], ua['X']], dims=['T', 'Y', 'X'], name='u')
v_R1 = xr.DataArray(np.real(np.fft.ifft2(vf_wave[2, :, :, :], axes=(0, 2))), coords=[va['T'], va['Y'], va['X']], dims=['T', 'Y', 'X'], name='v')
z_R1 = xr.DataArray(np.real(np.fft.ifft2(zf_wave[2, :, :, :], axes=(0, 2))), coords=[za['T'], za['Y'], za['X']], dims=['T', 'Y', 'X'], name='z')

u_R1.attrs['long_name'] = 'n = 1 Equatorial Rossby Waves in '+ str(lev) + ' hPa Zonal Wind'
v_R1.attrs['long_name'] = 'n = 1 Equatorial Rossby Waves in '+ str(lev) + ' hPa Meridional Wind'
z_R1.attrs['long_name'] = 'n = 1 Equatorial Rossby Waves in '+ str(lev) + ' hPa Geopotential Height'

u_R1.attrs['units'] = 'm/s'
v_R1.attrs['units'] = 'm/s'
z_R1.attrs['units'] = 'm'

u_R2 = xr.DataArray(np.real(np.fft.ifft2(uf_wave[3, :, :, :], axes=(0, 2))), coords=[ua['T'], ua['Y'], ua['X']], dims=['T', 'Y', 'X'], name='u')
v_R2 = xr.DataArray(np.real(np.fft.ifft2(vf_wave[3, :, :, :], axes=(0, 2))), coords=[va['T'], va['Y'], va['X']], dims=['T', 'Y', 'X'], name='v')
z_R2 = xr.DataArray(np.real(np.fft.ifft2(zf_wave[3, :, :, :], axes=(0, 2))), coords=[za['T'], za['Y'], za['X']], dims=['T', 'Y', 'X'], name='z')

u_R2.attrs['long_name'] = 'n = 2 Equatorial Rossby Waves in '+ str(lev) + ' hPa Zonal Wind'
v_R2.attrs['long_name'] = 'n = 2 Equatorial Rossby Waves in '+ str(lev) + ' hPa Meridional Wind'
z_R2.attrs['long_name'] = 'n = 2 Equatorial Rossby Waves in '+ str(lev) + ' hPa Geopotential Height'

u_R2.attrs['units'] = 'm/s'
v_R2.attrs['units'] = 'm/s'
z_R2.attrs['units'] = 'm' 

到此熱帶波動投影的程序已經結束。以下也可以輸出資料：

u_Kelvin.to_netcdf('./data/yhs03_ccews/u.850.kelvin.2016-2020.nc',unlimited_dims='T')
v_Kelvin.to_netcdf('./data/yhs03_ccews/v.850.kelvin.2016-2020.nc',unlimited_dims='T')
z_Kelvin.to_netcdf('./data/yhs03_ccews/z.850.kelvin.2016-2020.nc',unlimited_dims='T')

u_WMRG.to_netcdf('./data/yhs03_ccews/u.850.wmrg.2016-2020.nc',unlimited_dims='T')
v_WMRG.to_netcdf('./data/yhs03_ccews/v.850.wmrg.2016-2020.nc',unlimited_dims='T')
z_WMRG.to_netcdf('./data/yhs03_ccews/z.850.wmrg.2016-2020.nc',unlimited_dims='T')

u_R1.to_netcdf('./data/yhs03_ccews/u.850.r1.2016-2020.nc',unlimited_dims='T')
v_R1.to_netcdf('./data/yhs03_ccews/v.850.r1.2016-2020.nc',unlimited_dims='T')
z_R1.to_netcdf('./data/yhs03_ccews/z.850.r1.2016-2020.nc',unlimited_dims='T')

u_R2.to_netcdf('./data/yhs03_ccews/u.850.r2.2016-2020.nc',unlimited_dims='T')
v_R2.to_netcdf('./data/yhs03_ccews/v.850.r2.2016-2020.nc',unlimited_dims='T')
z_R2.to_netcdf('./data/yhs03_ccews/z.850.r2.2016-2020.nc',unlimited_dims='T')

我們繪製2019年12月16日這4種波動的平面圖。

import xarray as xr 
import numpy as np 
import cmaps
from matplotlib import pyplot as plt
from cartopy import crs as ccrs   
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

lon_formatter = LONGITUDE_FORMATTER
lat_formatter = LATITUDE_FORMATTER

time = '2019-12-15 12:00:00'
lon1, lon2 =  30, 210
lats, latn = -20,  20
lb = ['a','b','d','c','e']

# Define the figure and each axis for the 3 rows and 3 columns
proj = ccrs.PlateCarree(central_longitude=180) 
fig, axes = plt.subplots(nrows=3,ncols=2,
                         subplot_kw={'projection': proj}, 
                         figsize=(12,8))

# axs is a 2 dimensional array of `GeoAxes`.  We will flatten it into a 1-D array
ax = axes.flatten()

#### Plot Total fields 
cf1 = (za.sel(T=time,X=slice(lon1,lon2))
         .plot.contourf(x="X", y="Y", 
                        levels=[-25,-20,-15,-10,-5,0,5,10,15,20,25],
                        cmap='RdYlBu',
                        add_colorbar=False, extend='both',
                        transform=ccrs.PlateCarree(),
                        ax=ax[0]))
tot_wnd = (xr.merge([ua,va]).sel(T=time,X=slice(lon1,lon2))
                            .interp(X=np.arange(31,211,5),Y=np.arange(-19,24,5)))
qv = (tot_wnd.plot.quiver(ax=ax[0],
                          transform=ccrs.PlateCarree(),
                          x='X', y='Y',
                          u='u', v='v',
                          width=0.003 ,headaxislength=3,headlength=6,headwidth=7,
                          scale=200, colors="black",add_guide=False
                          ))
ax[0].set_title('')
ax[0].set_title("(a) Total Fields",loc='left')
qk = 5
qv_key = ax[0].quiverkey(qv,0.95,1.05,qk,(str(qk)+' m s$^{-1}$'),labelpos='N', labelsep =0.05, color='black')
    
# KW
cf = (z_Kelvin.sel(T=time,X=slice(lon1,lon2))
              .plot.contourf(x="X", y="Y", 
                             levels=[-5,-4,-3,-2,-1,-0.5,0.5,1,2,3,4,5],
                             cmap='RdYlBu',
                             add_colorbar=False, extend='both',
                             transform=ccrs.PlateCarree(),
                             ax=ax[2]))
kw_wnd = (xr.merge([u_Kelvin,v_Kelvin]).sel(T=time,X=slice(lon1,lon2))
                                       .interp(X=np.arange(31,211,7),Y=np.arange(-19,24,5)))
qv = (kw_wnd.plot.quiver(ax=ax[2],
                         transform=ccrs.PlateCarree(),
                         x='X', y='Y',
                         u='u', v='v',
                         width=0.005 ,headaxislength=3,headlength=6,headwidth=7,
                         scale=30, colors="black",add_guide=False
                         ))
ax[2].set_title('')
ax[2].set_title('(b) Kelvin Waves',loc='left')

qk = 2
qv_key = ax[2].quiverkey(qv,0.95,1.05,qk,(str(qk)+' m s$^{-1}$'),labelpos='N', labelsep =0.05, color='black')
    
# R1
cf = (z_R1.sel(T=time,X=slice(lon1,lon2))
           .plot.contourf(x="X", y="Y", 
                          levels=[-5,-4,-3,-2,-1,-0.5,0.5,1,2,3,4,5],
                          cmap='RdYlBu',
                          add_colorbar=False, extend='both',
                          transform=ccrs.PlateCarree(),
                          ax=ax[3]))
r1_wnd = (xr.merge([u_R1,v_R1]).sel(T=time,X=slice(lon1,lon2))
                               .interp(X=np.arange(31,211,7),Y=np.arange(-19,24,5)))
qv = (r1_wnd.plot.quiver(ax=ax[3],
                         transform=ccrs.PlateCarree(),
                         x='X', y='Y',
                         u='u', v='v',
                         width=0.005 ,headaxislength=3,headlength=6,headwidth=7,
                         scale=30, colors="black",add_guide=False
                         ))
ax[3].set_title('')
ax[3].set_title('(d) R1 waves',loc='left')

qk = 2
qv_key = ax[3].quiverkey(qv,0.95,1.05,qk,(str(qk)+' m s$^{-1}$'),labelpos='N', labelsep =0.05, color='black')

# R2
cf = (z_R2.sel(T=time,X=slice(lon1,lon2))
           .plot.contourf(x="X", y="Y", 
                          levels=[-5,-4,-3,-2,-1,-0.5,0.5,1,2,3,4,5],
                          cmap='RdYlBu',
                          add_colorbar=False, extend='both',
                          transform=ccrs.PlateCarree(),
                          ax=ax[4]))
r2_wnd = (xr.merge([u_R2,v_R2]).sel(T=time,X=slice(lon1,lon2))
                               .interp(X=np.arange(31,211,7),Y=np.arange(-19,24,5)))
qv = (r2_wnd.plot.quiver(ax=ax[4],
                         transform=ccrs.PlateCarree(),
                         x='X', y='Y',
                         u='u', v='v',
                         width=0.005 ,headaxislength=3,headlength=6,headwidth=7,
                         scale=30, colors="black",add_guide=False
                         ))
ax[4].set_title('')
ax[4].set_title('(d) R2 waves',loc='left')

# WMRG
cf = (z_WMRG.sel(T=time,X=slice(lon1,lon2))
           .plot.contourf(x="X", y="Y", 
                          levels=[-5,-4,-3,-2,-1,-0.5,0.5,1,2,3,4,5],
                          cmap='RdYlBu',
                          add_colorbar=False, extend='both',
                          transform=ccrs.PlateCarree(),
                          ax=ax[5]))
r1_wnd = (xr.merge([u_WMRG,v_WMRG]).sel(T=time,X=slice(lon1,lon2))
                               .interp(X=np.arange(31,211,7),Y=np.arange(-19,24,5)))
qv = (r1_wnd.plot.quiver(ax=ax[5],
                         transform=ccrs.PlateCarree(),
                         x='X', y='Y',
                         u='u', v='v',
                         width=0.005 ,headaxislength=3,headlength=6,headwidth=7,
                         scale=30, colors="black",add_guide=False
                         ))
ax[5].set_title('')
ax[5].set_title('(e) Westward MRG waves',loc='left')
    
qk = 2
qv_key = ax[4].quiverkey(qv,0.95,1.05,qk,(str(qk)+' m s$^{-1}$'),labelpos='N', labelsep =0.05, color='black')

for i in range(0,6):
    # Draw the coastines for each subplot
    ax[i].coastlines(color="dimgray")
    ax[i].set_extent([lon1,lon2,lats,latn], ccrs.PlateCarree())   
    gl = ax[i].gridlines(crs=ccrs.PlateCarree(), draw_labels=True, 
                         xlocs=np.arange(-180,240,30), ylocs=np.arange(-20, 30,10), 
                         x_inline=False, y_inline=False, linewidth=0)
    gl.right_labels = False
    gl.top_labels = False 
    ax[i].set_ylabel(' ')             # 設定坐標軸名稱。
    ax[i].set_xlabel(' ')
    ax[i].set_title(' ')

ax[1].remove()

plt.suptitle('GPH (Shade.), 850-hPa Wind Anml. (vectors)',y=0.95,size='large',weight='bold')

#plt.show()

Text(0.5, 0.95, 'GPH (Shade.), 850-hPa Wind Anml. (vectors)')

Error in callback <function _draw_all_if_interactive at 0x189a96680> (for post_execute):

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/pyplot.py:120, in _draw_all_if_interactive()
    118 def _draw_all_if_interactive():
    119     if matplotlib.is_interactive():
--> 120         draw_all()

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/_pylab_helpers.py:132, in Gcf.draw_all(cls, force)
    130 for manager in cls.get_all_fig_managers():
    131     if force or manager.canvas.figure.stale:
--> 132         manager.canvas.draw_idle()

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/backend_bases.py:2082, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
   2080 if not self._is_idle_drawing:
   2081     with self._idle_draw_cntx():
-> 2082         self.draw(*args, **kwargs)

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/backends/backend_agg.py:400, in FigureCanvasAgg.draw(self)
    396 # Acquire a lock on the shared font cache.
    397 with RendererAgg.lock, \
    398      (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
    399       else nullcontext()):
--> 400     self.figure.draw(self.renderer)
    401     # A GUI class may be need to update a window using this draw, so
    402     # don't forget to call the superclass.
    403     super().draw()

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     93 @wraps(draw)
     94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95     result = draw(artist, renderer, *args, **kwargs)
     96     if renderer._rasterizing:
     97         renderer.stop_rasterizing()

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     69     if artist.get_agg_filter() is not None:
     70         renderer.start_filter()
---> 72     return draw(artist, renderer)
     73 finally:
     74     if artist.get_agg_filter() is not None:

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/figure.py:3150, in Figure.draw(self, renderer)
   3147 finally:
   3148     self.stale = False
-> 3150 DrawEvent("draw_event", self.canvas, renderer)._process()

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/backend_bases.py:1263, in Event._process(self)
   1261 def _process(self):
   1262     """Generate an event with name ``self.name`` on ``self.canvas``."""
-> 1263     self.canvas.callbacks.process(self.name, self)

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/cbook/__init__.py:309, in CallbackRegistry.process(self, s, *args, **kwargs)
    307 except Exception as exc:
    308     if self.exception_handler is not None:
--> 309         self.exception_handler(exc)
    310     else:
    311         raise

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/cbook/__init__.py:96, in _exception_printer(exc)
     94 def _exception_printer(exc):
     95     if _get_running_interactive_framework() in ["headless", None]:
---> 96         raise exc
     97     else:
     98         traceback.print_exc()

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/cbook/__init__.py:304, in CallbackRegistry.process(self, s, *args, **kwargs)
    302 if func is not None:
    303     try:
--> 304         func(*args, **kwargs)
    305     # this does not capture KeyboardInterrupt, SystemExit,
    306     # and GeneratorExit
    307     except Exception as exc:

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/cartopy/mpl/gridliner.py:492, in Gridliner._draw_event(self, event)
    491 def _draw_event(self, event):
--> 492     self._draw_gridliner(renderer=event.renderer)

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/cartopy/mpl/gridliner.py:917, in Gridliner._draw_gridliner(self, nx, ny, renderer)
    915         visible = False
    916     kw = self._get_text_specs(segment_angle, loc, xylabel)
--> 917     kw['transform'] = self._get_padding_transform(
    918         segment_angle, loc, xylabel)
    919 kw.update(label_style)
    921 # Get x and y in data coords

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/cartopy/mpl/gridliner.py:1160, in Gridliner._get_padding_transform(self, padding_angle, loc, xylabel, padding_factor)
   1157 dy = padding_factor * padding * np.sin(padding_angle * np.pi / 180)
   1159 # Final transform
-> 1160 return mtrans.offset_copy(
   1161     self.axes.transData, fig=self.axes.figure,
   1162     x=dx, y=dy, units='points')

File ~/opt/anaconda3/envs/p3/lib/python3.10/site-packages/matplotlib/transforms.py:2974, in offset_copy(trans, fig, x, y, units)
   2972     return trans + Affine2D().translate(x, y)
   2973 if fig is None:
-> 2974     raise ValueError('For units of inches or points a fig kwarg is needed')
   2975 if units == 'points':
   2976     x /= 72.0

ValueError: For units of inches or points a fig kwarg is needed

_images/fac80019118904e1296da20b67a0b8db29f6cca3c30a2dcec1e46dff1c2d0d02.png

Caution

ERA5下載的資料是geopotential，但是計算熱帶波動時要使用geopotential height，請使用metpy套件進行轉換，如：

from metpy import units 
import metpy.calc

zha = (metpy.calc.geopotential_to_height(za * ((units.m)**2 / (units.s)**2))
            .metpy.dequantify())

參考文獻#

Gill, A. E., 1980: Some simple solutions for heat‐induced tropical circulation. Q. J. R. Meteorol. Soc., 106, 447–462, https://doi.org/10.1002/qj.49710644905.
Kiladis, G. N., M. C. Wheeler, P. T. Haertel, K. H. Straub, and P. E. Roundy, 2009: Convectively coupled equatorial waves. Rev. Geophys., 47, 1–42, https://doi.org/10.1029/2008RG000266.
Knippertz, P., and Coauthors, 2022: The intricacies of identifying equatorial waves. Q. J. R. Meteorol. Soc., 148, 2814–2852, https://doi.org/10.1002/qj.4338.
Matsuno, T., 1966: Quasi-Geostrophic Motions in the Equatorial Area. J. Meteorol. Soc. Japan. Ser. II, 44, 25–43, https://doi.org/10.2151/jmsj1965.44.1_25.
Tsai, W. Y., M.-M. Lu, C.-H. Sui, and P.-H. Lin, 2020: MJO and CCEW modulation on South China Sea and Maritime Continent boreal winter subseasonal peak precipitation. Terr. Atmos. Ocean. Sci., 31, 177–195, https://doi.org/10.3319/TAO.2019.10.28.01.
Wheeler, M. C., and G. N. Kiladis, 1999: Convectively Coupled Equatorial Waves: Analysis of Clouds and Temperature in the Wavenumber–Frequency Domain. J. Atmos. Sci., 56, 374–399, https://doi.org/10.1175/1520-0469(1999)056<0374:CCEWAO>2.0.CO;2.
Yang, G.-Y., B. Hoskins, and J. Slingo, 2003: Convectively Coupled Equatorial Waves: A New Methodology for Identifying Wave Structures in Observational Data. J. Atmos. Sci., 60, 1637–1654, https://doi.org/10.1175/1520-0469(2003)060<1637:CCEWAN>2.0.CO;2.

Tropical Waves Identification and Analysis

Contents

Tropical Waves Identification and Analysis#

熱帶波動理論簡介#

時空濾波 (Space-time filtering)#

雙曲柱狀函數投影法 (2DS-PCF)#

參考文獻#