climSIRS¶

Seceptible-Infectiou-Recovered-Sceptible (SIRS) model coupled to climate.

Wenchang Yang (wenchang@prnceton.edu)
Department of Geosciences, Princeton University

In [1]:

# init
%matplotlib inline
import sys
wython = '/tiger/scratch/gpfs/GEOCLIM/wenchang/wython'
if wython not in sys.path:
    sys.path.append(wython)
%run -im wystart
# https://www.dataquest.io/blog/jupyter-notebook-tips-tricks-shortcuts/
# from IPython.core.interactiveshell import InteractiveShell
# InteractiveShell.ast_node_interactivity = 'last' # options: 'all', 'none'
# https://stackoverflow.com/questions/41125690/matplotlib-notebook-showing-a-blank-histogram
# in case to switch back to the notebook backend, run the 2 lines below, 
# reload(plt)
# %matplotlib notebook
# import xaddon, xfilter, xtc

[imported]: os.path, sys, os, datetime, glob
[imported]: xarray as xr, numpy as np, matplotlib.pyplot as plt, pandas as pd
[imported]: from matplotlib.pyplot import figure, close
[executed]: plt.ion()
[config]: plt.rcParams['figure.dpi'] = 128
[config]: plt.rcParams['figure.figsize'] = [6.4, 6.4*9/16]:
[imported]: from misc.colormap_parula import parula
[config]: xr.set_options(cmap_sequential=parula)
[iPython config]: InlineBackend.figure_format ='retina'
[added]: PROJ_LIB = /home/wenchang/anaconda3/share/proj

In [2]:

from scipy.integrate import odeint
from scipy.integrate import solve_ivp

regions = dict(
    London=(51.5074, 0.1278),
    NYC=(40.7128, -74.0060+360),
    LA=(34.0522, -118.2437+360),
    Wuhan=(30.5928, 114.3055),
    Delhi=(28.7041, 77.1025),
    Maracaibo=(10.6427, -71.6125+360),
    Kinshasa=(-4.4419, 15.2663),
    Jakarta=(-6.2088, 106.8456),
    Johannesburg=(-26.2041, 28.0473),
    Buenos_Aires=(-34.6037, -58.3816+360),
    Melbourne=(-37.8136, 144.9631)
    
)

In [3]:

from numpy import pi as Pi
def get_q(case='cos', qmin=0.002, qmax=0.014, qrange=None, phi0=None, qyear=None):
    '''obtain specific humidity annual cycle'''
#     case = 'merra2'
#     qmin, qmax = 0.002, 0.018# units: kg/kg   
    if case == 'cos':
        if phi0 is None:
            phi0 = -Pi
        if qrange is None:
            qrange = qmax - qmin
        else: #overide qmax according to qmin and qrange
            qmax = qmin + qrange
        t = np.arange(0, 365)
        t = xr.DataArray(t, dims='time', coords=[t])
        t.attrs['units'] = 'day'
        qmean = (qmin + qmax)/2
        cycle = np.cos(2*Pi*t/365 + phi0) # normalized to have range of [0, 1]
        q = qmean + (qrange/2)*cycle
        q.name = 'q'
        q.attrs['units'] = 'kg/kg'
        s = f'cos-curve q: '
        if not hasattr(qmin, 'ndim'):
            s += f'qmin = {qmin}; '
        if not hasattr(qmax, 'ndim'):
            s += f'qmax = {qmax}; '
        if not hasattr(phi0, 'ndim'):
            s += f'phi0 = {phi0/Pi:.1f}xPi; '
        q.attrs['note'] = s
    elif case == 'era5':
        ifile = '/tigress/wenchang/data/era5/analysis_wy/2m_specific_humidity/daily/clim/q2m.run7clim.2010-2019.nc'
        q = ( 
            xr.open_dataset(ifile)
            ['q2m']
            .rename(dayofyear='time')
            .load()
        )
        q.attrs['note'] = case
    elif case == 'era5_single_year':
        ifiles = f'/tigress/wenchang/data/era5/analysis_wy/2m_specific_humidity/daily/era5.2m_specific_humidity.daily.{qyear}-??.nc'
        q = ( 
            xr.open_mfdataset(ifiles, combine='by_coords')
            ['q2m']
#             .load()
        )
        q.attrs['note'] = f'{case} {qyear}'
    elif case == 'merra2':
        ifile = '/tigress/gvecchi/ANALYSIS/DISEASE/COVID-19/merra_weekly_clim.nc'
        q2m = (
    xr.open_dataset('/tigress/gvecchi/ANALYSIS/DISEASE/COVID-19/merra_weekly_clim.nc')
    .rename(WEEK_QV2M='q2m', TT='time', LAT='latitude', LON='longitude')['q2m']
)
        q = ( 
            xr.open_dataset(ifile)
            .rename(WEEK_QV2M='q2m', TT='time', LAT='latitude', LON='longitude')
            ['q2m']
            .pipe(lambda x: x.assign_coords(time=x.time.dt.dayofyear))
            .reindex(time=np.arange(1, 366), method='nearest')
            .load()
        )
        # shift from 0-center to 180E-center
        if q.longitude.min() < 0:
            q = q.roll(longitude=q.longitude.size//2)
            q = q.assign_coords(
                longitude=q.longitude.where(q.longitude>=0, other=q.longitude+360).values
            )
        q.attrs['note'] = case
        
    return q
q = get_q()
q.plot()
plt.title(q.note)

# q.isel(time=0).plot(cmap='plasma')
# plt.title(case)
# plt.figure()
# for city,(latx, lonx) in regions.items():
#     q.sel(latitude=latx, longitude=lonx, method='nearest').plot(label=city)
# print()
# plt.legend()
# plt.title(case)

Out[3]:

Text(0.5, 1.0, 'cos-curve q: qmin = 0.002; qmax = 0.014; phi0 = -1.0xPi; ')

In [4]:

climdepens ={
    'HKU1': dict(a=-227.5, R0_min=1.4, R0_max=2.5),
    'OC43': dict(a=-32.5, R0_min=1.4, R0_max=2.5)
}
def q2R0(q, a, R0_min, R0_max):
    '''R0 dependence on climate (specific humidity)'''
    R0 = (R0_max - R0_min) * np.exp(a*q) + R0_min
    if isinstance(R0, xr.DataArray):
        R0.name = 'R0'
    
    return R0

In [5]:

def R0_at_t(R0, t, cycle=True):
    '''R0 as function of t (time) and specific humidity (q).
    q has dims of (n_time, n_grid).'''
    #R0_min, R0_max = 1.4, 2.5
    #a = -227.5 for HKU1 and -32.5 for OC43
    if cycle:# cycle through the first 365 daily specific humidity
        R0t = R0.isel(time=int( np.mod(t, 365) ), drop=True)
    else:
        R0t = R0.isel(time=int(t), drop=True)
    return R0t
# func_R0(0, R0t)

In [7]:

Ldis = {
    'HKU1': 66.25*7,
    'OC43': 62.5*7
}

How to do ODE in Python?¶

An example from scpy document: https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html.

In [8]:

# pendulum model
def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, 
            -b*omega - c*np.sin(theta)]
    return dydt

b, c = 0.25, 5.0
y0 = [np.pi - 0.1, 0.0]
t = np.linspace(0, 30, 30*10+1)
sol = odeint(pend, y0, t, args=(b, c))

# visualize
fig, axes = plt.subplots(1, 2, figsize=(6,3))

ax = axes[0]
ax.plot(t, sol[:, 0], 'C0', label='theta(t)')
ax.plot(t, sol[:, 1], 'C1', label='omega(t)')
ax.legend(loc='best')
ax.set_xlabel('t')
ax.grid('on')

ax = axes[1]
ax.scatter(sol[0, 0], sol[0, 1], color='k', marker='v')
ax.plot(sol[:, 0], sol[:, 1], color='C2')
ax.set_xlabel('theta(t)')
ax.set_ylabel('omega(t)')

plt.tight_layout()

In [28]:

# Lorenz system: https://en.wikipedia.org/wiki/Lorenz_system
# images: http://paulbourke.net/fractals/lorenz/
from mpl_toolkits.mplot3d import Axes3D

sigma = 10.0
rho = 28.0
beta = 8.0 / 3.0

def f(state, t):
    x, y, z = state  # Unpack the state vector
    return sigma * (y - x), x * (rho - z) - y, x * y - beta * z  # Derivatives

state0 = [1.0, 1.0, 1.0]
t = np.arange(0.0, 100.0, 0.01)

states = odeint(f, state0, t)

fig = plt.figure(figsize=(6,4))
ax = fig.gca(projection='3d')
ax.plot(states[:, 0], states[:, 1], states[:, 2], lw=0.5, color='C1')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('Lorenz System')
plt.tight_layout()

Apply to the SIRS model¶

Non-dimensionalized y and t¶

In [61]:

def model_nondimSIRS(yy, tt, tao, R0):
    '''non-dimensionized version of SIRS model from Baker et al., 2020.
    tao = L/D
    R0 = beta*D
    '''
    ss, ii = yy
    dydt = [
        (1-ss-ii)/tao - R0*ii*ss, 
        R0*ii*ss - ii
    ]
    return dydt

In [119]:

tao, R0 = 66.25*7/5, 2
y0 = [1.0-1e-5, 1e-5]
t = np.linspace(0, 5*365//5, 5*365//5*20+1)
sol = odeint(model_nondimSIRS, y0, t, args=(tao, R0))

In [120]:

fig, axes = plt.subplots(1, 2, figsize=(6,3))
ax = axes[0]
ax.plot(t, sol[:, 0], 'C1', label='S(t)/N')
ax.plot(t, sol[:, 1], 'C3', label='I(t)/N')
ax.legend(loc='best')
ax.set_xlabel('t[D]')
ax.grid('on')

ax = axes[1]
ax.scatter(sol[0, 0], sol[0, 1], color='k', marker='^')
ax.plot(sol[:, 0], sol[:, 1], color='C0')
ax.set_xlabel('S(t)/N')
ax.set_ylabel('I(t)/N')

plt.tight_layout()

In [121]:

def sirsplot(R0=2, tao=66.25*7/5, I0=1e-5, T=None, nt=None):
    #tao, R0 = 65*7/5, 1.5
    if T is None:
        T = 5*365//5
    if nt is None:
        nt = T*5+1
    y0 = [1.0-I0, I0]
    t = np.linspace(0, T, nt)
    sol = odeint(model_nondimSIRS, y0, t, args=(tao, R0))

    fig, axes = plt.subplots(1, 2, figsize=(6,3))
    ax = axes[0]
    ax.plot(t, sol[:, 0], 'C1', label='S(t)/N')
    ax.plot(t, sol[:, 1], 'C3', label='I(t)/N')
    ax.legend(loc='best')
    ax.set_xlabel('t[D]')
    ax.grid('on')

    ax = axes[1]
    ax.scatter(sol[0, 0], sol[0, 1], color='k', marker='^')
    ax.plot(sol[:, 0], sol[:, 1], color='C0')
    ax.set_xlabel('S(t)/N')
    ax.set_ylabel('I(t)/N')
    ax.text(0.98, 0.98, f'R0 = {R0}\ntao = {tao}\nI0={I0}',
            ha='right', va='top', transform=ax.transAxes)

    plt.tight_layout()
sirsplot()

In [75]:

sirsplot(R0=1.8)

In [17]:

sirsplot(R0=1.5)

In [18]:

sirsplot(I0=1e-4)

Non-dimensionalized y only¶

In [92]:

def model_SIRS(y, t, R0, L, D):
    '''Normalized version SIRS model from Baker et al., 2020.
    L: immunity length, 60 weeks
    D: infection period, 5 days
    R0 = beta*D
    '''
    ss, ii = y
    dydt = [
        (1-ss-ii)/L - R0*ii*ss/D, 
        R0*ii*ss/D - ii/D
    ]
    return dydt

In [94]:

def run_SIRS(R0=2, L=66.25*7, D=5, I0=1e-5, T=10*365, dt=1):
#     R0 = 2
#     L = 60*7
#     D = 5
#     I0 = 1e-5
#     T = 10*365
#     dt = 1
    
    y0 = [1.0-I0, I0]
    t = np.linspace(0, T, int(T/dt)+1)
    sol = odeint(model_SIRS, y0, t, args=(R0, L, D))
    
    da = xr.DataArray(sol, dims=['t', 'si'], coords=[t, ['S/N', 'I/N']])
    da.t.attrs['units'] = 'day'
    
    return da
da = run_SIRS()
da.sel(si='I/N').plot()

Out[94]:

[<matplotlib.lines.Line2D at 0x7ff515a43490>]

In [95]:

def siplot(R0=2, L=66.25*7, D=5, I0=1e-5, T=10*365, dt=1):
    #tao, R0 = 65*7/5, 1.5

#     y0 = [1.0-I0, I0]
#     t = np.linspace(0, T, int(T/dt)+1)
#     sol = odeint(SIRS_RH, y0, t, args=(R0, L, D))
#     da = xr.DataArray(sol, dims=['t', 'y'], coords=[t, ['S/N', 'I/N']])
#     da.t.attrs['units'] = 'day'
    
    da = run_SIRS(R0=R0, L=L, D=D, I0=I0, T=T, dt=dt)
    df = da.to_pandas()

    fig, axes = plt.subplots(1, 2, figsize=(7, 3))
    ax = axes[0]
    df.plot(ax=ax, secondary_y='I/N')
    xticks = np.arange(0, t[-1]+1, 365)
    ax.set_xticks(xticks)
    ax.set_xticklabels((xticks/365).astype('int'))
    ax.set_xlabel('time [year]')
#     ax.grid('on')

    ax = axes[1]
    df.plot(ax=ax, x='S/N', y='I/N', color='C2', legend=False)
    ax.scatter(sol[0, 0], sol[0, 1], color='k', marker='^')
    ax.set_xlabel('S/N')
    ax.set_ylabel('I/N')
#     ax.grid('on')
    ax.text(0.98, 0.98, f'R0 = {R0}\nL = {L}\nD= {D}\nI0={I0}',
            ha='right', va='top', transform=ax.transAxes)
#     ax.set_ylim(0, None)

    plt.tight_layout()
siplot()

In [96]:

siplot(R0=1.5)

In [24]:

siplot(R0=5)

In [25]:

siplot(R0=5, dt=0.5)

In [26]:

siplot(R0=1.1, T=10*365)

In [27]:

siplot(R0=1.2, T=10*365)

sensitivity to R0¶

In [97]:

R0s = np.arange(2, 1.2, -0.1)
das = []
for R0 in R0s:
    das.append(run_SIRS(R0=R0, T=5*365, I0=1e-5, D=5))
da = xr.concat(das, dim=pd.Index(R0s, name='R0'))

fig, ax = plt.subplots(1,1,figsize=(6,3))
colors = plt.cm.Spectral(np.linspace(0,1,da.R0.size))
for R0,c in zip(da.R0, colors):
    da.sel(si='I/N', R0=R0).plot(color=c, label=f'R0 = {R0.item():.1f}')
# da.sel(vname='I/N').plot(hue='R0')
ax.set_title('')
ax.set_ylabel('I/N')
ax.set_xlim(0, 365*1)

plt.legend()
plt.tight_layout()

sensitivity to D¶

In [98]:

Ds = np.arange(5, 21, 5)
das = []
for D in Ds:
    das.append(run_SIRS(R0=2, T=5*365, I0=1e-5, D=D))
da = xr.concat(das, dim=pd.Index(Ds, name='D'))

fig, ax = plt.subplots(1,1,figsize=(6,3))
colors = plt.cm.Spectral(np.linspace(0,1,da.D.size))
for D,c in zip(da.D, colors):
    da.sel(si='I/N', D=D).plot(color=c, label=f'D = {D.item()}')
ax.set_title('')
ax.set_ylabel('I/N')
plt.xlim(0, 365*3)
plt.legend()
plt.tight_layout()

sensitivity to L¶

In [99]:

Ls = 7*np.arange(60, 30, -5)
das = []
for L in Ls:
    das.append(run_SIRS(R0=2, T=5*365, I0=1e-5, L=L))
da = xr.concat(das, dim=pd.Index(Ls, name='L'))

fig, ax = plt.subplots(1,1,figsize=(6,3))
colors = plt.cm.Spectral_r(np.linspace(0,1,da.L.size))
for L,c in zip(da.L, colors):
    da.sel(si='I/N', L=L).plot(color=c, label=f'L = {L.item()}')
# da.sel(vname='I/N').plot(hue='R0')
ax.set_title('')
ax.set_ylabel('I/N')
plt.xlim(0, 365*3)
plt.legend()
plt.tight_layout()

sensitivity to dt¶

In [219]:

fig, ax = plt.subplots(1,1,figsize=(6,3))
dts = [1/8, 1/4, 1/2, 1, 1*2, 1*4]
colors = plt.cm.Spectral_r(np.linspace(0,1,len(dts)))
for dt,c in zip(dts, colors):
    da = run_SIRS(dt=dt)
    da.sel(si='I/N').sel(t=slice(52, 62)).plot(color=c, label=f'dt = {dt}')
# da.sel(vname='I/N').plot(hue='R0')
ax.set_title('')
ax.set_ylabel('I/N')
# plt.xlim(50, 60)
plt.legend()
plt.tight_layout()

SIRS model: vectorized version¶

Do multiple Ros simutaneously.

In [170]:

def model_vecSIRS(y, t, R0, L, D):
    '''Normalized version SIRS model from Baker et al., 2020.
    L: immunity length, 60 weeks
    D: infection period, 5 days
    R0 = beta*D
    '''
    N = y.size
    ss = y[:N//2]
    ii = y[N//2:]
    dydt = np.array([
        (1-ss-ii)/L - R0*ii*ss/D, 
        R0*ii*ss/D - ii/D
    ]).ravel()
    return dydt

In [176]:

def run_vecSIRS(R0=2, L=60*7, D=5, I0=1e-5, T=10*365, dt=1):
#     R0 = 2
#     L = 60*7
#     D = 5
#     I0 = 1e-5
#     T = 10*365
#     dt = 1
    
    zeros = I0*L*D*R0*0
    y0 = np.array([1.0-I0+zeros, I0+zeros]).ravel()
    t = np.linspace(0, T, int(T/dt)+1)
    sol = odeint(model_vecSIRS, y0, t, args=(R0, L, D))
    sol = sol.reshape(t.size, 2, y0.size//2)
    
    da = xr.DataArray(sol, dims=['t', 'si', 'case'],coords=[t, ['S/N', 'I/N'], np.arange(1, y0.size//2+1)])
    da.t.attrs['units'] = 'day'
    
    return da
R0=np.arange(2,1,-0.1)
I0 = R0*0 + 1e-5
da = run_SIRSvec(I0=I0, R0=R0)
da = da.rename(case='R0').assign_coords(R0=R0)
print(da)

<xarray.DataArray (t: 3651, si: 2, R0: 10)>
array([[[9.99990000e-01, 9.99990000e-01, 9.99990000e-01, ...,
         9.99990000e-01, 9.99990000e-01, 9.99990000e-01],
        [1.00000000e-05, 1.00000000e-05, 1.00000000e-05, ...,
         1.00000000e-05, 1.00000000e-05, 1.00000000e-05]],

       [[9.99985566e-01, 9.99985833e-01, 9.99986094e-01, ...,
         9.99987322e-01, 9.99987554e-01, 9.99987780e-01],
        [1.22181222e-05, 1.19750938e-05, 1.17371103e-05, ...,
         1.06184326e-05, 1.04081083e-05, 1.02019920e-05]],

       [[9.99980153e-01, 9.99980847e-01, 9.99981514e-01, ...,
         9.99984484e-01, 9.99985012e-01, 9.99985520e-01],
        [1.49287123e-05, 1.43406275e-05, 1.37762151e-05, ...,
         1.12751164e-05, 1.08328692e-05, 1.04080593e-05]],

       ...,

       [[5.00072619e-01, 5.26291764e-01, 5.55417489e-01, ...,
         7.68322689e-01, 8.32731461e-01, 9.09011175e-01],
        [5.87643030e-03, 5.58596293e-03, 5.21756192e-03, ...,
         2.71041919e-03, 2.00627133e-03, 1.03960376e-03]],

       [[5.00073455e-01, 5.26289210e-01, 5.55420363e-01, ...,
         7.68326464e-01, 8.32724019e-01, 9.09017431e-01],
        [5.87660201e-03, 5.58590918e-03, 5.21730531e-03, ...,
         2.70978066e-03, 2.00597975e-03, 1.03958624e-03]],

       [[5.00074251e-01, 5.26286679e-01, 5.55423276e-01, ...,
         7.68330355e-01, 8.32716656e-01, 9.09023674e-01],
        [5.87677564e-03, 5.58585003e-03, 5.21705415e-03, ...,
         2.70914498e-03, 2.00568465e-03, 1.03957015e-03]]])
Coordinates:
  * t        (t) float64 0.0 1.0 2.0 3.0 ... 3.648e+03 3.649e+03 3.65e+03
  * si       (si) <U3 'S/N' 'I/N'
  * R0       (R0) float64 2.0 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1

In [177]:

R0 = np.arange(2, 1.2, -0.1)
I0 = 1e-5
da = run_vecSIRS(I0=I0, R0=R0)
da = da.rename(case='R0').assign_coords(R0=R0)

fig, ax = plt.subplots(1,1,figsize=(6,3))
colors = plt.cm.Spectral(np.linspace(0,1,da.R0.size))
for R0,c in zip(da.R0, colors):
    da.sel(si='I/N', R0=R0).plot(color=c, label=f'R0 = {R0.item():.1f}')
# da.sel(vname='I/N').plot(hue='R0')
ax.set_title('')
ax.set_ylabel('I/N')
ax.set_xlim(0, 365*1)

plt.legend()
plt.tight_layout()

climSIRS: SIRS coupled to climate¶

In [48]:

fig = plt.figure(figsize=(5, 4))
q = np.arange(0, 20+1)
dis = 'OC43'
R0 = q2R0(q/1000, **climdepens[dis])
plt.plot(q, R0, label=f'{dis}')

dis = 'HKU1'
R0 = q2R0(q/1000, **climdepens[dis])
plt.plot(q, R0, label=f'{dis}')

plt.legend()
plt.xlabel('q[g/kg]')
plt.ylabel('R0')
plt.tight_layout()

q and R0 annual cycle: idealized case¶

In [76]:

# idealized q annual cycle
dis = 'HKU1'
q = get_q()
R0 = q2R0(q, **climdepens[dis])
df = pd.DataFrame({'q': q, 'R0': R0}, index=pd.Index(R0.time, name='t[day]'))

df.plot(secondary_y='R0')
plt.title(q.note)

Out[76]:

Text(0.5, 1.0, 'cos-curve q: qmin = 0.002; qmax = 0.014; phi0 = 1.0xPi')

In [77]:

# idealized q annual cycle
dis = 'HKU1'
q = get_q(phi0=0)
R0 = q2R0(q, **climdepens[dis])
df = pd.DataFrame({'q': q, 'R0': R0}, index=pd.Index(R0.time, name='t[day]'))

df.plot(secondary_y='R0')
plt.title(q.note)

Out[77]:

Text(0.5, 1.0, 'cos-curve q: qmin = 0.002; qmax = 0.014; phi0 = 0.0xPi')

In [78]:

# idealized q annual cycle
dis = 'HKU1'
q = get_q(phi0=-np.pi/2)
R0 = q2R0(q, **climdepens[dis])
df = pd.DataFrame({'q': q, 'R0': R0}, index=pd.Index(R0.time, name='t[day]'))

df.plot(secondary_y='R0')
plt.title(q.note)

Out[78]:

Text(0.5, 1.0, 'cos-curve q: qmin = 0.002; qmax = 0.014; phi0 = -0.5xPi')

In [79]:

# idealized q annual cycle
dis = 'HKU1'
q = get_q(phi0=-np.pi*3/2)
R0 = q2R0(q, **climdepens[dis])
df = pd.DataFrame({'q': q, 'R0': R0}, index=pd.Index(R0.time, name='t[day]'))

df.plot(secondary_y='R0')
plt.title(q.note)

Out[79]:

Text(0.5, 1.0, 'cos-curve q: qmin = 0.002; qmax = 0.014; phi0 = -1.5xPi')

realistic q¶

In [142]:

q_era5 = get_q('era5')
q = q_era5
q.isel(time=0).plot(cmap='plasma')
plt.title(q.note)

Out[142]:

Text(0.5, 1.0, 'era5')

In [141]:

q_merra2 = get_q('merra2')
q = q_merra2
q.isel(time=0).plot(cmap='plasma')
plt.title(q.note)

/home/wenchang/anaconda3/lib/python3.7/site-packages/xarray/core/dataarray.py:2823: FutureWarning: roll_coords will be set to False in the future. Explicitly set roll_coords to silence warning.
  shifts=shifts, roll_coords=roll_coords, **shifts_kwargs

Out[141]:

Text(0.5, 1.0, 'merra2')

In [145]:

q = q_era5
for city,(latx, lonx) in regions.items():
    q.sel(latitude=latx, longitude=lonx, method='nearest').plot(label=city)
print()
plt.legend(bbox_to_anchor=(1,1), loc='upper left')
plt.title(q.note)
plt.tight_layout()
print(plt.ylim())

(0.0014669001218862831, 0.02305866334354505)

In [146]:

q = q_merra2
for city,(latx, lonx) in regions.items():
    q.sel(latitude=latx, longitude=lonx, method='nearest').plot(label=city)
print()
plt.legend(bbox_to_anchor=(1,1), loc='upper left')
plt.title(q.note)
plt.ylim(0.0014669001218862831, 0.02305866334354505)
plt.tight_layout()

In [80]:

# realistic q annual cycle
dis = 'HKU1'
city = 'NYC'
latx, lonx = regions[city]
q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
R0 = q2R0(q, **climdepens[dis])
df = pd.DataFrame({f'q_{city}': q, f'R0_{dis}': R0})

df.plot(secondary_y=f'R0_{dis}')
plt.title(q.note)

Out[80]:

Text(0.5, 1.0, 'era5')

In [82]:

# realistic q annual cycle
dis = 'OC43'
city = 'NYC'
latx, lonx = regions[city]
q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
R0 = q2R0(q, **climdepens[dis])
df = pd.DataFrame({f'q_{city}': q, f'R0_{dis}': R0})

df.plot(secondary_y=f'R0_{dis}')
plt.title(q.note)

Out[82]:

Text(0.5, 1.0, 'era5')

climSIRS using `odeint` backend¶

In [8]:

def model_climSIRS(y, t, L, D, R0):
    '''SIRS model right hand side from Baker et al., 2020.
    Input:
    -------
        y: y[0:y.size//2] is S/N; y[y.size//2:] is I/N
        t: time (units is 'day')
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        dydy: half d(S/N)/dt, half d(I/N)/dt
    '''
    N = y.size
    ss = y[:N//2]
    ii = y[N//2:]
    if hasattr(R0, 'time'): # R0 is time-dependent
        _R0 = R0_at_t(R0, t)
        if _R0.ndim < 1:
            _R0 = _R0.item()
    else:
        _R0 = R0
    dydt = np.array([
        (1-ss-ii)/L - _R0*ii*ss/D, 
        _R0*ii*ss/D - ii/D
    ]).ravel()
    
    return dydt

In [9]:

from scipy.integrate import odeint
def run_climSIRS(ii0=1e-5, ss0=None, T=5*365, dt=1, L=None, D=5, R0=2, dis='HKU1'):
    '''Integrate SIRS model from Baker et al., 2020.

    Input:
    -------
        I0: initial I/N
        T: time stop (start is 1; units is 'day')
        dt: time step
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        ds: solution for the SIRS model.
    '''
    if L is None and dis is not None:
        L = Ls[dis]
    if hasattr(R0, 'time'): # R0 time-dependent
        if R0.ndim > 1:# 
            grid = R0.stack(grid=R0.isel(time=0).dims).grid # grid info is saved
            zeros = np.zeros(grid.size) # broadcast array by R0.dims[1:]
        else: # only time dimension
            zeros = ii0*L*D*0 # broadcast array by ii0*L*D*a
    else: # R0 time-independent
        zeros = ii0*L*D*R0*0
    # default initial S/N value determined by initial I/N value
    if ss0 is None:
        ss0 = 1 - ii0
    # broadcast initial values by adding zeros
    ii0 = ii0 + zeros
    ss0 = ss0 + zeros
        
    y0 = np.hstack([ss0, ii0])
    t = np.arange(0, T, dt)
    sol = odeint(model_climSIRS, y0, t, args=(L, D, R0))
    ss = sol[:, 0:y0.size//2]
    ii = sol[:, y0.size//2:]
    da_ss = xr.DataArray(ss, dims=['t', 'grid'],coords=[t, np.arange(1, y0.size//2+1)])
    da_ss.attrs['long_name'] = 'S/N'
    da_ii = xr.DataArray(ii, dims=['t', 'grid'],coords=[t, np.arange(1, y0.size//2+1)])
    da_ii.attrs['long_name'] = 'I/N'
    
    if hasattr(R0, 'time') and R0.ndim > 1: # R0 has more than "time" dimensions
        da_ss = da_ss.assign_coords(grid=grid).unstack()
        da_ii = da_ii.assign_coords(grid=grid).unstack()
    
    ds = xr.Dataset(dict(ss=da_ss, ii=da_ii, year=da_ss.t/365))
    ds.t.attrs['units'] = 'day'
    
    return ds

latx, lonx = regions['NYC']
# q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
# q = get_q()
ds = run_climSIRS(R0=2, T=2*365, ii0=1/8e6, L=66.25*7)
ds.ii.plot()
ds

In [10]:

#R0 is a ndarray
dis='HKU1'
R0 = np.arange(2,1.3,-0.2)
ds = run_climSIRS(dis=dis, ii0=1/8e6, R0=R0, T=3*365)
ds = ds.assign_coords(grid=R0).rename(grid='R0')
ds.ii.assign_coords(t=ds.year).rename(t='year').plot(hue='R0')

Out[10]:

[<matplotlib.lines.Line2D at 0x7f81fab48110>,
 <matplotlib.lines.Line2D at 0x7f81fab48390>,
 <matplotlib.lines.Line2D at 0x7f81fab48550>,
 <matplotlib.lines.Line2D at 0x7f81fab48710>]

idealized q¶

In [37]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 1*365

# q is from ERA5 NYC location
cities = ['NYC', 'Johannesburg']
colors = ['k', 'k']
# city = 'NYC'
for city,c in zip(cities, colors):
    latx, lonx = regions[city]
    color = c
    if city == 'NYC':
        label = f'{city}({dis})'
        ls = '-'
    else:
        label = f'{city}'
        ls = '--'
    q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
    q.plot(ax=axes[0], color=color, ls=ls)

    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], color=color, ls=ls)

    ds = run_climSIRS(R0=R0, T=T, ii0=1/8e6, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                         color=color,
                                                         label=label,
                                                         ls=ls
                                                        )
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['qpeakJul', 'qpeakOct', 'qpeakOct', 'qpeakApr']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRS(R0=R0, T=T, ii0=1/8e6, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

The cell below is the same as the above except we run the four cases together. For reasons not quite understood, it is very slow and takes four times time to run.

In [53]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 3.5*365

# q is from ERA5 NYC location
city = 'NYC'
latx, lonx = regions[city]
color = 'k'
q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
q.plot(ax=axes[0], color=color)

R0 = q2R0(q, **climdepens[dis])
R0.plot(ax=axes[1], color=color)

ds = run_climSIRS(R0=R0, T=T, ii0=1/8e6, dis=dis)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                     color=color,
                                                     label=f'{dis}_{city}')
#q is idealized
from numpy import pi as Pi
qpeakmons = ['Jul', 'Oct', 'Jan', 'Apr']
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
phi0s = xr.DataArray(phi0s, dims='qpeakmon', coords=[qpeakmons])
qs = get_q(phi0=phi0s)
R0s = q2R0(qs, **climdepens[dis])
ds = run_climSIRS(R0=R0s, T=T, ii0=1/8e6, dis=dis)
for mon in qpeakmons:
    q = qs.sel(qpeakmon=mon)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = R0s.sel(qpeakmon=mon)
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds.sel(qpeakmon=mon).ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

In [135]:

# NYC q vs. idealized q for 'OC43'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'OC43'
T = 2*365

# q is from ERA5 NYC location
city = 'NYC'
latx, lonx = regions[city]
color = 'k'
q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
q.plot(ax=axes[0], color=color)

R0 = q2R0(q, **climdepens[dis])
R0.plot(ax=axes[1], color=color)

ds = run_climSIRS(q=q, T=T, ii0=1/8e6, dis=dis)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                     color=color,
                                                     label=f'{dis}_{city}')
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['Jul_qpeak', 'Oct_qpeak', 'Jan_qpeak', 'Apr_qpeak']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRS(q=q, T=T, ii0=1/8e6, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

idealized q in a warming climate¶

In [136]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 3.5*365
qfactor = 1.1

# q is from ERA5 NYC location
city = 'NYC'
latx, lonx = regions[city]
color = 'k'
q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
q.plot(ax=axes[0], color=color)
q.pipe(lambda q: q*qfactor).plot(ax=axes[0], color=color, ls='--')

R0 = q2R0(q, **climdepens[dis])
R0.plot(ax=axes[1], color=color)
q2R0(q*qfactor, **climdepens[dis]).plot(ax=axes[1], color=color, ls='--')

ds = run_climSIRS(q=q, T=T, ii0=1/8e6, dis=dis)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                     color=color,
                                                     label=f'{dis}_{city}')
run_climSIRS(q=q*qfactor, T=T, ii0=1/8e6, dis=dis) \
    .ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                     color=color,
                                                     ls='--')
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['Jul_qpeak', 'Oct_qpeak', 'Jan_qpeak', 'Apr_qpeak']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRS(q=q, T=T, ii0=1/8e6, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

In [145]:

names = ['NYC', 'Maracaibo', 'Delhi']
dss = []
for name in names:
    print(name)
    latx, lonx = regions[name]
    q = q2m.sel(latitude=[latx,], longitude=[lonx,], method='nearest')
    ds = run_climSIRS(q=q, T=5*365, I0=1/8e6, L=62.5*7, a=-32.5)
    ds.ii.plot()
    ds.ii.squeeze().isel(t=ds.ii.squeeze()==ds.ii.max())
    dss.append(ds.drop(['latitude', 'longitude']))
# ds_ = xr.concat(dss, dim=pd.Index(names, name='city'))
# ofile = '/tigress/wenchang/analysis/climate_and_disease/climSIRS/ClimSIRS.OC43.run7climq.cities.nc'
# print(ofile)
# ds_.to_netcdf(ofile)
# print('[saved]:', ofile)

NYC
Maracaibo
Delhi

In [146]:

# MERRA2 q
q2m = (
    xr.open_dataset('/tigress/gvecchi/ANALYSIS/DISEASE/COVID-19/merra_weekly_clim.nc')
    .rename(WEEK_QV2M='q2m', TT='time', LAT='latitude', LON='longitude')['q2m']
)
q2m = q2m.assign_coords(time=q2m.time.dt.dayofyear).reindex(time=np.arange(1, 366), method='nearest')
names = ['NYC', 'Maracaibo', 'Delhi']
dss = []
for name in names:
    print(name, end=' ')
    latx, lonx = regions[name]
    if lonx > 180:
        lonx -= 360
    print(latx, lonx)
    q = q2m.sel(latitude=[latx,], longitude=[lonx,], method='nearest')
    ds = run_climSIRS(q=q, T=5*365, I0=1/8e6, L=62.5*7, a=-32.5)
    ds.ii.plot()
    print( ds.ii.squeeze().isel(t=ds.ii.squeeze()==ds.ii.max()) )
    dss.append(ds.drop(['latitude', 'longitude']))
# ds_ = xr.concat(dss, dim=pd.Index(names, name='city'))
# ofile = '/tigress/wenchang/analysis/climate_and_disease/climSIRS/ClimSIRS.OC43.merra2wclimq.cities.nc'
# print(ofile)
# ds_.to_netcdf(ofile)
# print('[saved]:', ofile)

NYC 40.7128 -74.00599999999997
<xarray.DataArray 'ii' (t: 1)>
array([0.21764435])
Coordinates:
  * t          (t) float64 60.0
    latitude   float64 40.5
    longitude  float64 -73.75
Attributes:
    long_name:  I/N
Maracaibo 10.6427 -71.61250000000001
<xarray.DataArray 'ii' (t: 1)>
array([0.16936215])
Coordinates:
  * t          (t) float64 75.0
    latitude   float64 10.5
    longitude  float64 -71.88
Attributes:
    long_name:  I/N
Delhi 28.7041 77.1025
<xarray.DataArray 'ii' (t: 1)>
array([0.21486067])
Coordinates:
  * t          (t) float64 60.0
    latitude   float64 28.5
    longitude  float64 76.88
Attributes:
    long_name:  I/N

climSIRS using `solve_ivp` backend¶

In [10]:

# t instead of y is the first argument required by scipy.integrate.solve_ivp
def model_climSIRSpy(t, y, L, D, R0):
    return model_climSIRS(y=y, t=t, L=L, D=D, R0=R0)

In [26]:

# scipy.integrate.solve_ivp backend
from scipy.integrate import solve_ivp
def run_climSIRSpy(ii0=1e-5, ss0=None, T=5*365, dt=1, L=None, D=5, R0=2, dis='HKU1', ode_kws=None):
    '''Integrate SIRS model from Baker et al., 2020.

    Input:
    -------
        I0: initial I/N
        T: time stop (start is 1; units is 'day')
        dt: time step
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        ds: solution for the SIRS model.
    '''
    if ode_kws is None:
        ode_kws = {'method': 'RK45'}
    if L is None and dis is not None:
        L = Ls[dis]
    if hasattr(R0, 'time'): # R0 time-dependent
        if R0.ndim > 1:# 
            R0 = R0.stack(grid=R0.isel(time=0).dims)
            grid = R0.grid # grid info is saved
            zeros = np.zeros(grid.size) # broadcast array by R0.dims[1:]
        else: # only time dimension
            zeros = ii0*L*D*0 # broadcast array by ii0*L*D*a
    else: # R0 time-independent
        zeros = ii0*L*D*R0*0
    # default initial S/N value determined by initial I/N value
    if ss0 is None:
        ss0 = 1 - ii0
    # broadcast initial values by adding zeros
    ii0 = ii0 + zeros
    ss0 = ss0 + zeros
        
    y0 = np.hstack([ss0, ii0])
    t = np.arange(0, T, dt)
#     # scipy.integrate.odeint backend; call FORTRAN code ultimately
#     sol = odeint(model_climSIRS, y0, t, args=(L, D, R0))
#     ss = sol[:, 0:y0.size//2]
#     ii = sol[:, y0.size//2:]
#     da_ss = xr.DataArray(ss, dims=['t', 'grid'],coords=[t, np.arange(1, y0.size//2+1)])
#     da_ss.attrs['long_name'] = 'S/N'
#     da_ii = xr.DataArray(ii, dims=['t', 'grid'],coords=[t, np.arange(1, y0.size//2+1)])
#     da_ii.attrs['long_name'] = 'I/N'
    
    # scipy.integrate.solve_ivp backend, pure Python code
#     print(t)
    sol = solve_ivp(model_climSIRSpy,[0, T], y0, t_eval=t,  args=(L, D, R0), **ode_kws)
    if sol.status == -1:
        method = ode_kws['method']
        print(f'[failed]: solve_ivp(method="{method}")')
        sys.exit()
#     print(sol)
    ss = sol.y[0:y0.size//2, :]
    ii = sol.y[y0.size//2:, :]
#     print(ss.shape)
    da_ss = xr.DataArray(ss, dims=['grid', 't'],coords=[np.arange(1, y0.size//2+1), t])
    da_ss.attrs['long_name'] = 'S/N'
    da_ii = xr.DataArray(ii, dims=['grid', 't'],coords=[np.arange(1, y0.size//2+1), t])
    da_ii.attrs['long_name'] = 'I/N'
    
    if hasattr(R0, 'time') and R0.ndim > 1: # R0 has more than "time" dimensions
        da_ss = da_ss.assign_coords(grid=grid).unstack()
        da_ii = da_ii.assign_coords(grid=grid).unstack()
    
    ds = xr.Dataset(dict(ss=da_ss, ii=da_ii, year=da_ss.t/365))
    ds.t.attrs['units'] = 'day'
    
    return ds

latx, lonx = regions['NYC']
# q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
# q = get_q()
ds = run_climSIRSpy(R0=2, T=2*365, ii0=1/8e6, L=66.25*7)
ds.ii.plot()
ds

In [19]:

#R0 is a ndarray
dis='HKU1'
R0 = np.arange(2.4,1.1,-0.4)
ds = run_climSIRSpy(dis=dis, ii0=1/8e6, R0=R0, T=3*365)
ds = ds.assign_coords(grid=[f'{r:.1f}' for r in R0]).rename(grid='R0')
ds.ii.assign_coords(t=ds.year).rename(t='year').plot(hue='R0')

Out[19]:

[<matplotlib.lines.Line2D at 0x7f81fa6ec490>,
 <matplotlib.lines.Line2D at 0x7f81fa6ec6d0>,
 <matplotlib.lines.Line2D at 0x7f81fa6ec890>,
 <matplotlib.lines.Line2D at 0x7f81fa6eca50>]

idealized q¶

In [20]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 2*365
ii0 = 1e-5

# q is from ERA5 NYC location
cities = ['NYC', 'Johannesburg']
colors = ['k', 'k']
# city = 'NYC'
for city,c in zip(cities, colors):
    latx, lonx = regions[city]
    color = c
    if city == 'NYC':
        label = f'{city}({dis})'
        ls = '-'
    else:
        label = f'{city}'
        ls = '--'
    q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
    q.plot(ax=axes[0], color=color, ls=ls)

    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], color=color, ls=ls)

    ds = run_climSIRS(R0=R0, T=T, ii0=ii0, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                         color=color,
                                                         label=label,
                                                         ls=ls
                                                        )
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['qpeakJul', 'qpeakOct', 'qpeakOct', 'qpeakApr']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRS(R0=R0, T=T, ii0=ii0, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

In [27]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 2*365
ii0 = 1e-5

# q is from ERA5 NYC location
cities = ['NYC', 'Johannesburg']
colors = ['k', 'k']
# city = 'NYC'
for city,c in zip(cities, colors):
    latx, lonx = regions[city]
    color = c
    if city == 'NYC':
        label = f'{city}({dis})'
        ls = '-'
    else:
        label = f'{city}'
        ls = '--'
    q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
    q.plot(ax=axes[0], color=color, ls=ls)

    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], color=color, ls=ls)

    ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis, ode_kws=None)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                         color=color,
                                                         label=label,
                                                         ls=ls
                                                        )
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['qpeakJul', 'qpeakOct', 'qpeakOct', 'qpeakApr']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis, ode_kws=None)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

[failed]: solve_ivp(method="RK45")

An exception has occurred, use %tb to see the full traceback.

SystemExit

In [19]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 2*365
ii0 = 1e-5

# q is from ERA5 NYC location
cities = ['NYC', 'Johannesburg']
colors = ['k', 'k']
# city = 'NYC'
for city,c in zip(cities, colors):
    latx, lonx = regions[city]
    color = c
    if city == 'NYC':
        label = f'{city}({dis})'
        ls = '-'
    else:
        label = f'{city}'
        ls = '--'
    q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
    q.plot(ax=axes[0], color=color, ls=ls)

    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], color=color, ls=ls)

    ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis, ode_kws={'method': 'LSODA'})
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                         color=color,
                                                         label=label,
                                                         ls=ls
                                                        )
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['qpeakJul', 'qpeakOct', 'qpeakOct', 'qpeakApr']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis, ode_kws={'method': 'LSODA'})
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

In [21]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 2*365
ii0 = 1e-5

# q is from ERA5 NYC location
cities = ['NYC', 'Johannesburg']
colors = ['k', 'k']
# city = 'NYC'
for city,c in zip(cities, colors):
    latx, lonx = regions[city]
    color = c
    if city == 'NYC':
        label = f'{city}({dis})'
        ls = '-'
    else:
        label = f'{city}'
        ls = '--'
    q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
    q.plot(ax=axes[0], color=color, ls=ls)

    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], color=color, ls=ls)

    ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis, ode_kws={'method': 'DOP853'})
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                         color=color,
                                                         label=label,
                                                         ls=ls
                                                        )
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['qpeakJul', 'qpeakOct', 'qpeakOct', 'qpeakApr']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis, ode_kws={'method': 'DOP853'})
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

In [13]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 2*365

# q is from ERA5 NYC location
cities = ['NYC', 'Johannesburg']
colors = ['k', 'k']
# city = 'NYC'
for city,c in zip(cities, colors):
    latx, lonx = regions[city]
    color = c
    if city == 'NYC':
        label = f'{city}({dis})'
        ls = '-'
    else:
        label = f'{city}'
        ls = '--'
    q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
    q.plot(ax=axes[0], color=color, ls=ls)

    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], color=color, ls=ls)

    ds = run_climSIRSpy(R0=R0, T=T, ii0=1/8e6, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                         color=color,
                                                         label=label,
                                                         ls=ls
                                                        )
#q is idealized
from numpy import pi as Pi
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
labels = ['qpeakJul', 'qpeakOct', 'qpeakOct', 'qpeakApr']
for phi0,label in zip(phi0s, labels):
    q = get_q(phi0=phi0)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds = run_climSIRSpy(R0=R0, T=T, ii0=1/8e6, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=label)
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

idealized q: vectorized form¶

In [14]:

# NYC q vs. idealized q for 'HKU1'
fig, axes = plt.subplots(3, 1, figsize=(6,6))
dis = 'HKU1'
T = 1*365

# q is from ERA5 NYC location
cities = ['NYC', 'Johannesburg']
colors = ['k', 'k']
# city = 'NYC'
for city,c in zip(cities, colors):
    latx, lonx = regions[city]
    color = c
    if city == 'NYC':
        label = f'{city}({dis})'
        ls = '-'
    else:
        label = f'{city}'
        ls = '--'
    q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
    q.plot(ax=axes[0], color=color, ls=ls)

    R0 = q2R0(q, **climdepens[dis])
    R0.plot(ax=axes[1], color=color, ls=ls)

    ds = run_climSIRSpy(R0=R0, T=T, ii0=1/8e6, dis=dis)
    ds.ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2], 
                                                         color=color,
                                                         label=label,
                                                         ls=ls
                                                        )
#q is idealized
from numpy import pi as Pi
qpeakmons = ['Jul', 'Oct', 'Jan', 'Apr']
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
phi0s = xr.DataArray(phi0s, dims='qpeakmon', coords=[qpeakmons])
qs = get_q(phi0=phi0s)
R0s = q2R0(qs, **climdepens[dis])
ds = run_climSIRSpy(R0=R0s, T=T, ii0=1/8e6, dis=dis)
for mon in qpeakmons:
    q = qs.sel(qpeakmon=mon)
    q.plot(ax=axes[0], lw=1)
    axes[0].set_ylabel('q[kg/kg] annual cycle')
    
    R0 = R0s.sel(qpeakmon=mon)
    R0.plot(ax=axes[1], lw=1)
    axes[1].set_ylabel('R0 annual cycle')
    
    ds.sel(qpeakmon=mon).ii.assign_coords(t=ds.year).rename(t='year').plot(ax=axes[2],
                                                         lw=1, 
                                                         label=f'qpeak{mon}')
    axes[2].set_title('')


for i in range(2):
    axes[i].set_title('')
axes[2].legend(loc='upper right')
plt.tight_layout()

NYC q: impact of different years¶

In [38]:

city = 'NYC'
latx, lonx = regions[city]
qyears = range(2010, 2020)
q = [get_q('era5_single_year', qyear=qyear)
     .sel(latitude=latx, longitude=lonx, method='nearest')
     .isel(time=slice(0, 365))
     .pipe(lambda x: x.assign_coords(time=x.time.dt.dayofyear))
     .load()
     for qyear in qyears]
q = xr.concat(q, dim=pd.Index(qyears, name='qyear'))

qclim = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest').isel(time=slice(0, 365))
qclim

# q.rename(time='day_of_year').plot(hue='qyear')
# qclim.rename(time='day_of_year').plot(color='k', lw=2)

Out[38]:

xarray.DataArray

'q2m'

time: 365

0.0029838255 0.0028502066 0.0027963307 ... 0.0032594036 0.0031316439

array([0.00298383, 0.00285021, 0.00279633, 0.0026232 , 0.00249139,
       0.00244834, 0.00247756, 0.00253146, 0.00278112, 0.00289252,
       0.00296091, 0.00301259, 0.0030684 , 0.0031    , 0.00310943,
       0.00293562, 0.00287054, 0.00285219, 0.00282139, 0.00278855,
       0.00281523, 0.00276455, 0.00271965, 0.00272953, 0.00274707,
       0.00275541, 0.00274354, 0.00261775, 0.00258906, 0.00261268,
       0.00255315, 0.00256932, 0.00262775, 0.00262192, 0.00269351,
       0.00276731, 0.00269903, 0.00266943, 0.00266811, 0.00262026,
       0.00257257, 0.00253425, 0.00253323, 0.00267819, 0.00269291,
       0.0026987 , 0.00272116, 0.0028233 , 0.00296684, 0.0029939 ,
       0.00302629, 0.0032451 , 0.00343834, 0.00343027, 0.00338666,
       0.00333835, 0.00337992, 0.00336824, 0.00318758, 0.00299252,
       0.00298523, 0.00297473, 0.00293393, 0.00287977, 0.002873  ,
       0.00301059, 0.00319288, 0.00336579, 0.00346135, 0.00361224,
       0.00371678, 0.00374277, 0.00369648, 0.0035979 , 0.00354492,
       0.0035714 , 0.00359213, 0.00363137, 0.00363051, 0.00362291,
       0.00366346, 0.0037292 , 0.00370551, 0.00363743, 0.00367925,
       0.00382889, 0.00393062, 0.00402769, 0.00411019, 0.00430024,
       0.0044209 , 0.00437195, 0.00430771, 0.00438187, 0.00450932,
       0.00452559, 0.00448783, 0.00456437, 0.00472051, 0.00492918,
       0.00510272, 0.00510014, 0.00519037, 0.00538268, 0.00539052,
       0.00544057, 0.00549308, 0.00542517, 0.00550712, 0.00559975,
       0.00564346, 0.00573824, 0.00586848, 0.00597063, 0.00609938,
       0.00607843, 0.00598045, 0.00598353, 0.0061977 , 0.0063144 ,
       0.00644392, 0.0065402 , 0.00676817, 0.00701586, 0.00721704,
       0.00719994, 0.00726655, 0.00728235, 0.00731255, 0.00719592,
       0.00721992, 0.00736933, 0.00760355, 0.00773339, 0.00781016,
       0.00800058, 0.00832558, 0.00844975, 0.00849005, 0.00859783,
       0.00874655, 0.00892145, 0.00918967, 0.00948166, 0.00989131,
       0.01025824, 0.0105838 , 0.01090292, 0.01125845, 0.01130643,
       0.01119524, 0.0109133 , 0.01077728, 0.01060782, 0.01023766,
       0.0099617 , 0.00990987, 0.01001441, 0.01031724, 0.01036136,
       0.01056542, 0.01078143, 0.01072756, 0.01078174, 0.01098623,
       0.01118382, 0.01136664, 0.01131971, 0.01148064, 0.01180747,
       0.0121467 , 0.01225646, 0.0122347 , 0.01225388, 0.01229249,
       0.01234406, 0.01246966, 0.01237108, 0.01235381, 0.01254093,
       0.01275897, 0.01295836, 0.01305358, 0.01329406, 0.01356693,
       0.01382581, 0.01387897, 0.01385755, 0.01391336, 0.01402459,
       0.0139499 , 0.01394221, 0.01398067, 0.0140339 , 0.01416924,
       0.01428529, 0.0144612 , 0.01463841, 0.01466296, 0.01458733,
       0.01464051, 0.014631  , 0.01444313, 0.01413336, 0.01390837,
       0.01391606, 0.01392933, 0.01370429, 0.01351061, 0.01365661,
       0.01382541, 0.01386655, 0.01384586, 0.01374391, 0.01378906,
       0.01389483, 0.01377601, 0.01371613, 0.01363686, 0.01353133,
       0.01360953, 0.01373505, 0.01366582, 0.01365621, 0.01366529,
       0.01378994, 0.01389873, 0.01385606, 0.0137314 , 0.01379114,
       0.01394298, 0.01374797, 0.01329587, 0.01297569, 0.01270335,
       0.01260903, 0.01251966, 0.01223417, 0.01218839, 0.01241615,
       0.01254056, 0.01277539, 0.01291791, 0.0129063 , 0.01292148,
       0.01296604, 0.01283056, 0.01270987, 0.0124273 , 0.01219967,
       0.01212194, 0.01212507, 0.01196113, 0.01182092, 0.01160805,
       0.01140808, 0.01107915, 0.01073387, 0.01048678, 0.01042214,
       0.01041515, 0.01047954, 0.01048504, 0.01058835, 0.01063359,
       0.01055985, 0.01053594, 0.0106807 , 0.01071155, 0.01072162,
       0.01059966, 0.0104521 , 0.01024613, 0.01010307, 0.00973915,
       0.00940753, 0.00936124, 0.00927971, 0.00927077, 0.00924642,
       0.0091296 , 0.00893104, 0.00871158, 0.00847014, 0.00829788,
       0.00810115, 0.0078265 , 0.00747684, 0.00730121, 0.00731439,
       0.00711373, 0.00700222, 0.00682969, 0.00679187, 0.00685345,
       0.00680081, 0.00670139, 0.0066748 , 0.00667486, 0.00657382,
       0.00658131, 0.00648684, 0.00647207, 0.00634179, 0.00610777,
       0.00590697, 0.00596324, 0.00582966, 0.00555605, 0.00531352,
       0.00519838, 0.00509774, 0.00495738, 0.00465495, 0.00443468,
       0.00442723, 0.00449963, 0.00444945, 0.00443073, 0.00439212,
       0.00442283, 0.00439132, 0.00440861, 0.00424705, 0.00412554,
       0.00409017, 0.00407872, 0.00412035, 0.00416685, 0.00411738,
       0.00415   , 0.004297  , 0.00437403, 0.00437578, 0.00427223,
       0.00430755, 0.00444595, 0.00440173, 0.00413143, 0.00401931,
       0.00397831, 0.00392137, 0.00382088, 0.00364447, 0.00367903,
       0.00381641, 0.00373187, 0.00365768, 0.00361897, 0.00356612,
       0.00345909, 0.0034393 , 0.00356116, 0.00377187, 0.00386057,
       0.00393001, 0.00393223, 0.00402173, 0.0040508 , 0.0038275 ,
       0.00363359, 0.00349167, 0.00330802, 0.0032594 , 0.00313164],
      dtype=float32)

Coordinates: (3)
- latitude
  ()
  float32
  40.75
  units :
  degrees_north
  long_name :
  latitude
```
array(40.75, dtype=float32)
```
- longitude
  ()
  float32
  286.0
  units :
  degrees_east
  long_name :
  longitude
```
array(286., dtype=float32)
```
- time
  (time)
  int64
  1 2 3 4 5 6 ... 361 362 363 364 365
```
array([  1,   2,   3, ..., 363, 364, 365])
```
Attributes: (1)
note :
era5

In [39]:

q.rename(time='day_of_year').plot(hue='qyear', lw=1)
qclim.rename(time='day_of_year').plot(color='k', lw=2)

Out[39]:

[<matplotlib.lines.Line2D at 0x7f81f3d04b50>]

In [40]:

dis = 'HKU1'
R0 = q2R0(q, **climdepens[dis])
R0.rename(time='day_of_year').plot(hue='qyear', lw=1)
R0clim = q2R0(qclim, **climdepens[dis])
R0clim.rename(time='day_of_year').plot(color='k', lw=2)

Out[40]:

[<matplotlib.lines.Line2D at 0x7f81cc33a7d0>]

In [41]:

dis = 'HKU1'
T = 365
ii0 = 1/8e6
ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis)

ds_clim = run_climSIRSpy(R0=R0clim, T=T, ii0=ii0, dis=dis)

In [42]:

tspan = slice(0,150)
ds.ii.sel(t=tspan).plot(hue='qyear', lw=1)
ds.ii.sel(t=tspan).mean('qyear').plot(color='k', ls='--', label='ens_mean')
ds_clim.ii.sel(t=tspan).plot(color='k', lw=2, label='2010-2019_clim_q')
plt.legend()
plt.title(f'{city} {dis}')

Out[42]:

Text(0.5, 1.0, 'NYC HKU1')

In [59]:

ds.ii.max('t').pipe(lambda x: x.sortby(x)).to_dataframe().rename(columns={'ii':'I/N peak'})

Out[59]:

	longitude	latitude	I/N peak
qyear
2012	286.0	40.75	0.114647
2016	286.0	40.75	0.120275
2011	286.0	40.75	0.127920
2010	286.0	40.75	0.130802
2019	286.0	40.75	0.142942
2013	286.0	40.75	0.148758
2015	286.0	40.75	0.152011
2018	286.0	40.75	0.152466
2017	286.0	40.75	0.156920
2014	286.0	40.75	0.162966

In [60]:

ds.ii.argmax('t').pipe(lambda x: x.sortby(x)+1).to_dataframe().rename(columns={'ii':'I/N peak day'})

Out[60]:

	longitude	latitude	I/N peak day
qyear
2015	286.0	40.75	74
2010	286.0	40.75	78
2014	286.0	40.75	78
2011	286.0	40.75	80
2019	286.0	40.75	80
2013	286.0	40.75	83
2017	286.0	40.75	83
2018	286.0	40.75	84
2016	286.0	40.75	88
2012	286.0	40.75	89

Wuhan q: impact of different years¶

In [61]:

city = 'Wuhan'
latx, lonx = regions[city]
qyears = range(2010, 2020)
q = [get_q('era5_single_year', qyear=qyear)
     .sel(latitude=latx, longitude=lonx, method='nearest')
     .isel(time=slice(0, 365))
     .pipe(lambda x: x.assign_coords(time=x.time.dt.dayofyear))
     .load()
     for qyear in qyears]
q = xr.concat(q, dim=pd.Index(qyears, name='qyear'))

qclim = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest').isel(time=slice(0, 365))
qclim

# q.rename(time='day_of_year').plot(hue='qyear')
# qclim.rename(time='day_of_year').plot(color='k', lw=2)

Out[61]:

xarray.DataArray

'q2m'

time: 365

0.0038499625 0.0039108004 0.0039654123 ... 0.0037823406 0.00381002

array([0.00384996, 0.0039108 , 0.00396541, 0.00396307, 0.00395205,
       0.00391889, 0.00387438, 0.00383735, 0.00381875, 0.00382306,
       0.00388626, 0.00391503, 0.00393172, 0.00394593, 0.00401043,
       0.00410113, 0.0041013 , 0.00401649, 0.00395115, 0.00386708,
       0.00382464, 0.0037199 , 0.00365428, 0.00369913, 0.00384173,
       0.00396997, 0.00412374, 0.00422879, 0.00435451, 0.00442368,
       0.00440288, 0.00434718, 0.00433435, 0.00431745, 0.00425412,
       0.0041962 , 0.00411805, 0.00407051, 0.00406043, 0.00404898,
       0.00401627, 0.00405741, 0.00414838, 0.00428756, 0.00441097,
       0.00445234, 0.00450543, 0.00455471, 0.00463519, 0.00471951,
       0.00484009, 0.00505468, 0.00529486, 0.00549488, 0.00568399,
       0.00579704, 0.00573618, 0.00564946, 0.00556006, 0.00553768,
       0.00554225, 0.0054817 , 0.00546671, 0.00564948, 0.00572215,
       0.00566538, 0.00565144, 0.00567487, 0.005853  , 0.00598109,
       0.00603175, 0.00619931, 0.00654893, 0.00683868, 0.00711363,
       0.00713937, 0.00717249, 0.00723806, 0.00728204, 0.00715414,
       0.00694373, 0.00676188, 0.00689906, 0.00715784, 0.00741637,
       0.00765254, 0.00798636, 0.00835469, 0.00857449, 0.00871476,
       0.00869505, 0.00866334, 0.00858696, 0.00845183, 0.00847593,
       0.00848362, 0.00848194, 0.00864285, 0.00870671, 0.00877133,
       0.00895137, 0.00909177, 0.0093564 , 0.00963205, 0.00981403,
       0.00994658, 0.01020619, 0.01034047, 0.01031488, 0.01028643,
       0.01034132, 0.01041152, 0.01048678, 0.01044185, 0.01050339,
       0.01073979, 0.01105092, 0.01112098, 0.01112334, 0.01129035,
       0.01151918, 0.01183238, 0.01194634, 0.01201998, 0.0122097 ,
       0.01240856, 0.01252815, 0.01257031, 0.01253231, 0.01250316,
       0.01261722, 0.01260217, 0.01264037, 0.01283118, 0.01304882,
       0.0132553 , 0.01347077, 0.01340064, 0.01335973, 0.01333084,
       0.01342168, 0.01355568, 0.01359165, 0.01356292, 0.01366283,
       0.01375375, 0.01394483, 0.0140325 , 0.01418591, 0.01430564,
       0.01440433, 0.0144713 , 0.01476181, 0.01499403, 0.01512287,
       0.01530661, 0.0155499 , 0.01588971, 0.01607925, 0.01615192,
       0.01613336, 0.01614031, 0.0161065 , 0.01609515, 0.01615035,
       0.01644783, 0.01683641, 0.01725201, 0.01751891, 0.01776231,
       0.01805865, 0.01829231, 0.0182891 , 0.01816111, 0.01811236,
       0.01826091, 0.01838889, 0.01847037, 0.0185464 , 0.01876002,
       0.01900109, 0.01920494, 0.01924849, 0.01919977, 0.01912993,
       0.01915605, 0.01917076, 0.01920425, 0.01926471, 0.01932807,
       0.01944965, 0.01956671, 0.01959944, 0.01960965, 0.01966114,
       0.01966151, 0.01971671, 0.01978462, 0.01994552, 0.02017635,
       0.02041342, 0.02054929, 0.02067195, 0.0208062 , 0.02090543,
       0.02088418, 0.02081405, 0.02071447, 0.02066008, 0.02058296,
       0.02050532, 0.02037441, 0.02029664, 0.02029199, 0.02036418,
       0.02031861, 0.02023637, 0.02017793, 0.0201592 , 0.02007518,
       0.01994652, 0.01980468, 0.01978154, 0.019794  , 0.0197068 ,
       0.01968314, 0.01970389, 0.01963492, 0.01957898, 0.01951839,
       0.01936146, 0.0192409 , 0.0189236 , 0.01865689, 0.01847629,
       0.01827631, 0.01801163, 0.01778408, 0.01743499, 0.01724254,
       0.01709413, 0.01696677, 0.01669925, 0.01641185, 0.01619044,
       0.01611168, 0.01601595, 0.0157966 , 0.01565691, 0.01556407,
       0.0155923 , 0.01563768, 0.01564743, 0.01565   , 0.01565729,
       0.0155642 , 0.0154402 , 0.01523413, 0.0150273 , 0.01475264,
       0.01438931, 0.01399675, 0.01363965, 0.01338767, 0.01322121,
       0.01312092, 0.01314725, 0.01320297, 0.01303107, 0.01287117,
       0.01269304, 0.01242866, 0.0119828 , 0.01154264, 0.01124177,
       0.01123203, 0.01124135, 0.0113008 , 0.0113511 , 0.01139181,
       0.01136164, 0.01119908, 0.0110256 , 0.0107394 , 0.010419  ,
       0.01018723, 0.01001457, 0.00984796, 0.00980664, 0.00985171,
       0.00994064, 0.01007627, 0.0100907 , 0.0100717 , 0.01000659,
       0.00974984, 0.009393  , 0.00911422, 0.00875764, 0.00844768,
       0.00824416, 0.00810408, 0.0080731 , 0.00809514, 0.00806193,
       0.00806805, 0.00810057, 0.008039  , 0.00797968, 0.00793945,
       0.0078195 , 0.0076629 , 0.00750322, 0.00735444, 0.0073238 ,
       0.00736559, 0.00741939, 0.00744461, 0.00732479, 0.00726426,
       0.00726942, 0.00724636, 0.00709801, 0.00685384, 0.00658375,
       0.00648584, 0.00638247, 0.00619096, 0.00603817, 0.00593836,
       0.0058708 , 0.00586869, 0.00583336, 0.00571238, 0.00558259,
       0.00540299, 0.00520421, 0.00502213, 0.00487062, 0.00472002,
       0.00461028, 0.0044898 , 0.00441945, 0.00440583, 0.00444166,
       0.00437731, 0.00428358, 0.00421722, 0.0042081 , 0.00418991,
       0.00420388, 0.00421938, 0.00429102, 0.0043707 , 0.00439542,
       0.00440566, 0.00434589, 0.00424379, 0.00414334, 0.00407297,
       0.00399193, 0.00388142, 0.00374956, 0.00378234, 0.00381002],
      dtype=float32)

Coordinates: (3)
- latitude
  ()
  float32
  30.5
  units :
  degrees_north
  long_name :
  latitude
```
array(30.5, dtype=float32)
```
- longitude
  ()
  float32
  114.25
  units :
  degrees_east
  long_name :
  longitude
```
array(114.25, dtype=float32)
```
- time
  (time)
  int64
  1 2 3 4 5 6 ... 361 362 363 364 365
```
array([  1,   2,   3, ..., 363, 364, 365])
```
Attributes: (1)
note :
era5

In [27]:

q.rename(time='day_of_year').plot(hue='qyear', lw=1)
qclim.rename(time='day_of_year').plot(color='k', lw=2)

Out[27]:

[<matplotlib.lines.Line2D at 0x7f81ec11ef50>]

In [62]:

dis = 'HKU1'
R0 = q2R0(q, **climdepens[dis])
R0.rename(time='day_of_year').plot(hue='qyear', lw=1)
R0clim = q2R0(qclim, **climdepens[dis])
R0clim.rename(time='day_of_year').plot(color='k', lw=2)

Out[62]:

[<matplotlib.lines.Line2D at 0x7f81cc1b7610>]

In [63]:

dis = 'HKU1'
T = 365
ii0 = 1/8e6
ds = run_climSIRSpy(R0=R0, T=T, ii0=ii0, dis=dis)

ds_clim = run_climSIRSpy(R0=R0clim, T=T, ii0=ii0, dis=dis)

In [37]:

tspan = slice(0,365)
ds.ii.sel(t=tspan).plot(hue='qyear', lw=1)
ds.ii.sel(t=tspan).mean('qyear').plot(color='k', ls='--', label='ens_mean')
ds_clim.ii.sel(t=tspan).plot(color='k', lw=2, label='2010-2019_clim_q')
plt.legend()
plt.title(f'{city} {dis}')

Out[37]:

Text(0.5, 1.0, 'Wuhan HKU1')

In [64]:

ds.ii.max('t').pipe(lambda x: x.sortby(x)).to_dataframe().rename(columns={'ii':'I/N peak'})

Out[64]:

	longitude	latitude	I/N peak
qyear
2014	114.25	30.5	0.073773
2016	114.25	30.5	0.074505
2017	114.25	30.5	0.077028
2019	114.25	30.5	0.077278
2018	114.25	30.5	0.080127
2013	114.25	30.5	0.086065
2015	114.25	30.5	0.086778
2010	114.25	30.5	0.090946
2012	114.25	30.5	0.093187
2011	114.25	30.5	0.110028

In [65]:

ds.ii.argmax('t').pipe(lambda x: x.sortby(x)+1).to_dataframe().rename(columns={'ii':'I/N peak day'})

Out[65]:

	longitude	latitude	I/N peak day
qyear
2011	114.25	30.5	90
2012	114.25	30.5	95
2018	114.25	30.5	101
2013	114.25	30.5	103
2014	114.25	30.5	104
2016	114.25	30.5	104
2019	114.25	30.5	104
2017	114.25	30.5	105
2010	114.25	30.5	106
2015	114.25	30.5	106

idealized q: qmin and qmax choice¶

In [15]:

# idealized q for 'HKU1'
dis = 'HKU1'
T = 2*365

#q is idealized
from numpy import pi as Pi
qpeakmons = ['Jul', 'Oct', 'Jan', 'Apr']
phi0s = [-Pi, -Pi*3/2, 0, -Pi/2,]
phi0s = xr.DataArray(phi0s, dims='qpeakmon', coords=[qpeakmons])
qmin = np.arange(0, 0.023, 0.002)
qmin = xr.DataArray(qmin, dims='qmin', coords=[qmin])
qmax = np.arange(0, 0.023, 0.002)
qmax = xr.DataArray(qmax, dims='qmax', coords=[qmax])
qs = get_q(qmin=qmin, qmax=qmax, phi0=phi0s)

R0s = q2R0(qs, **climdepens[dis])
ds = run_climSIRSpy(R0=R0s, T=T, ii0=1/8e6, dis=dis)

In [31]:

ds.ii.max('t') \
    .where(lambda x: x.qmax>=qmin).transpose() \
    .assign_attrs(long_name='I/N peak') \
    .plot.contourf(col='qpeakmon', col_wrap=2, cmap='plasma', levels=10)

Out[31]:

<xarray.plot.facetgrid.FacetGrid at 0x7fb57de725d0>

In [32]:

# qpeamon mean and anomaly
plt.figure()
ds.ii.max('t') \
    .where(lambda x: x.qmax>=qmin).transpose() \
    .mean('qpeakmon') \
    .assign_attrs(long_name='I/N peak mean over qpeakmon') \
    .plot.contourf(cmap='plasma', levels=21)
plt.tight_layout()

plt.figure()
ds.ii.max('t') \
    .where(lambda x: x.qmax>=qmin).transpose() \
    .pipe(lambda x: x - x.mean('qpeakmon')) \
    .assign_attrs(long_name='I/N peak anomaly') \
    .plot.contourf(col='qpeakmon', col_wrap=2, levels=21)

Out[32]:

<xarray.plot.facetgrid.FacetGrid at 0x7fb57db1b510>

<Figure size 819.2x460.8 with 0 Axes>

I/N peak day of year:

In [37]:

ds.ii.sel(t=slice(0, 365)).argmax('t').pipe(lambda x: x+1) \
    .where(lambda x: x.qmax>=qmin).transpose() \
    .assign_attrs(long_name='I/N peak day of year') \
    .plot.contourf(col='qpeakmon', col_wrap=2, cmap='plasma', levels=11)

Out[37]:

<xarray.plot.facetgrid.FacetGrid at 0x7fb57d4174d0>

In [39]:

# qpeamon mean and anomaly
plt.figure()
ds.ii.sel(t=slice(0, 365)).argmax('t').pipe(lambda x: x+1) \
    .where(lambda x: x.qmax>=qmin).transpose() \
    .mean('qpeakmon') \
    .assign_attrs(long_name='I/N peak day of year: mean over qpeakmon') \
    .plot.contourf(cmap='plasma', levels=11)
plt.tight_layout()

plt.figure()
ds.ii.sel(t=slice(0, 365)).argmax('t').pipe(lambda x: x+1) \
    .where(lambda x: x.qmax>=qmin).transpose() \
    .pipe(lambda x: x - x.mean('qpeakmon')) \
    .assign_attrs(long_name='I/N peak day of year anomaly') \
    .plot.contourf(col='qpeakmon', col_wrap=2, levels=13)

Out[39]:

<xarray.plot.facetgrid.FacetGrid at 0x7fb57ce86990>

<Figure size 819.2x460.8 with 0 Axes>

In [80]:

q = get_q('era5')

In [86]:

R0 = q2R0(q, **climdepens[dis])

In [87]:

R0.mean('time').plot(cmap='plasma')

Out[87]:

<matplotlib.collections.QuadMesh at 0x7f08d0f58250>

In [88]:

dis = 'HKU1'
T = 5
ds = run_climSIRSpy(R0=R0, T=T, ii0=1/8e6, dis=dis)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
/tiger/scratch/gpfs/GEOCLIM/wenchang/wython/wystart.py in <module>
      1 dis = 'HKU1'
      2 T = 5
----> 3 ds = run_climSIRSpy(R0=R0, T=T, ii0=1/8e6, dis=dis)

/tiger/scratch/gpfs/GEOCLIM/wenchang/wython/wystart.py in run_climSIRSpy(ii0, ss0, T, dt, L, D, R0, dis)
     49 
     50     # scipy.integrate.solve_ivp backend, pure Python code
---> 51     sol = solve_ivp(model_climSIRSp,[0, T], y0, t_eval=t,  args=(L, D, R0))
     52     ss = sol.y[0:y0.size//2, :]
     53     ii = sol.y[y0.size//2:, :]

~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/ivp.py in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)
    575     status = None
    576     while status is None:
--> 577         message = solver.step()
    578 
    579         if solver.status == 'finished':

~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/base.py in step(self)
    180         else:
    181             t = self.t
--> 182             success, message = self._step_impl()
    183 
    184             if not success:

~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/rk.py in _step_impl(self)
    144 
    145             y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A,
--> 146                                    self.B, self.C, self.K)
    147             scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol
    148             error_norm = self._estimate_error_norm(self.K, h, scale)

~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/rk.py in rk_step(fun, t, y, f, h, A, B, C, K)
     63     for s, (a, c) in enumerate(zip(A[1:], C[1:]), start=1):
     64         dy = np.dot(K[:s].T, a[:s]) * h
---> 65         K[s] = fun(t + c * h, y + dy)
     66 
     67     y_new = y + h * np.dot(K[:-1].T, B)

~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/base.py in fun(t, y)
    137         def fun(t, y):
    138             self.nfev += 1
--> 139             return self.fun_single(t, y)
    140 
    141         self.fun = fun

~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/base.py in fun_wrapped(t, y)
     19 
     20     def fun_wrapped(t, y):
---> 21         return np.asarray(fun(t, y), dtype=dtype)
     22 
     23     return fun_wrapped, y0

~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/ivp.py in <lambda>(t, x, fun)
    513         # additional parameters.  Pass in the original fun as a keyword
    514         # argument to keep it in the scope of the lambda.
--> 515         fun = lambda t, x, fun=fun: fun(t, x, *args)
    516         jac = options.get('jac')
    517         if callable(jac):

/tiger/scratch/gpfs/GEOCLIM/wenchang/wython/wystart.py in model_climSIRSp(t, y, L, D, R0)
      1 # t instead of y is the first argument required by scipy.integrate.solve_ivp
      2 def model_climSIRSp(t, y, L, D, R0):
----> 3     return model_climSIRS(y=y, t=t, L=L, D=D, R0=R0)

/tiger/scratch/gpfs/GEOCLIM/wenchang/wython/wystart.py in model_climSIRS(y, t, L, D, R0)
     26     dydt = np.array([
     27         (1-ss-ii)/L - _R0*ii*ss/D,
---> 28         _R0*ii*ss/D - ii/D
     29     ]).ravel()
     30 

~/anaconda3/lib/python3.7/site-packages/xarray/core/common.py in _iter(self)
    142     def _iter(self: Any) -> Iterator[Any]:
    143         for n in range(len(self)):
--> 144             yield self[n]
    145 
    146     def __iter__(self: Any) -> Iterator[Any]:

~/anaconda3/lib/python3.7/site-packages/xarray/core/dataarray.py in __getitem__(self, key)
    629         else:
    630             # xarray-style array indexing
--> 631             return self.isel(indexers=self._item_key_to_dict(key))
    632 
    633     def __setitem__(self, key: Any, value: Any) -> None:

~/anaconda3/lib/python3.7/site-packages/xarray/core/dataarray.py in isel(self, indexers, drop, **indexers_kwargs)
   1018             }
   1019             if coord_indexers:
-> 1020                 coord_value = coord_value.isel(coord_indexers)
   1021                 if drop and coord_value.ndim == 0:
   1022                     continue

~/anaconda3/lib/python3.7/site-packages/xarray/core/variable.py in isel(self, indexers, **indexers_kwargs)
   1058 
   1059         key = tuple(indexers.get(dim, slice(None)) for dim in self.dims)
-> 1060         return self[key]
   1061 
   1062     def squeeze(self, dim=None):

~/anaconda3/lib/python3.7/site-packages/xarray/core/variable.py in __getitem__(self, key)
    705         if new_order:
    706             data = duck_array_ops.moveaxis(data, range(len(new_order)), new_order)
--> 707         return self._finalize_indexing_result(dims, data)
    708 
    709     def _finalize_indexing_result(self: VariableType, dims, data) -> VariableType:

~/anaconda3/lib/python3.7/site-packages/xarray/core/variable.py in _finalize_indexing_result(self, dims, data)
   2132         if getattr(data, "ndim", 0) != 1:
   2133             # returns Variable rather than IndexVariable if multi-dimensional
-> 2134             return Variable(dims, data, self._attrs, self._encoding)
   2135         else:
   2136             return type(self)(dims, data, self._attrs, self._encoding, fastpath=True)

~/anaconda3/lib/python3.7/site-packages/xarray/core/variable.py in __init__(self, dims, data, attrs, encoding, fastpath)
    301             unrecognized encoding items.
    302         """
--> 303         self._data = as_compatible_data(data, fastpath=fastpath)
    304         self._dims = self._parse_dimensions(dims)
    305         self._attrs = None

~/anaconda3/lib/python3.7/site-packages/xarray/core/variable.py in as_compatible_data(data, fastpath)
    224     if isinstance(data, np.ndarray):
    225         if data.dtype.kind == "O":
--> 226             data = _possibly_convert_objects(data)
    227         elif data.dtype.kind == "M":
    228             data = np.asarray(data, "datetime64[ns]")

~/anaconda3/lib/python3.7/site-packages/xarray/core/variable.py in _possibly_convert_objects(values)
    160     datetime64 and timedelta64, according to the pandas convention.
    161     """
--> 162     return np.asarray(pd.Series(values.ravel())).reshape(values.shape)
    163 
    164 

~/anaconda3/lib/python3.7/site-packages/pandas/core/series.py in __init__(self, data, index, dtype, name, copy, fastpath)
    282                 if not is_list_like(data):
    283                     data = [data]
--> 284                 index = ibase.default_index(len(data))
    285             elif is_list_like(data):
    286 

~/anaconda3/lib/python3.7/site-packages/pandas/core/indexes/base.py in default_index(n)
   5388     from pandas.core.indexes.range import RangeIndex
   5389 
-> 5390     return RangeIndex(0, n, name=None)
   5391 
   5392 

~/anaconda3/lib/python3.7/site-packages/pandas/core/indexes/range.py in __new__(cls, start, stop, step, dtype, copy, name)
     81     # Constructors
     82 
---> 83     def __new__(
     84         cls, start=None, stop=None, step=None, dtype=None, copy=False, name=None,
     85     ):

KeyboardInterrupt:

xclimSIRS¶

In [89]:

def model_climSIRSwy(ss, ii, t, L, D, R0):
    '''SIRS model right hand side from Baker et al., 2020.

    Input:
    -------
        ss: S/N
        ii: I/N
        t: time (units is 'day')
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        tuple (dsdt, didt)
        dsdt: S/N tendency
        didt: I/N tendency
    '''
    if hasattr(R0, 'time'): # R0 is time-dependent
        _R0 = R0_at_t(R0, t)
        if _R0.ndim < 1:
            _R0 = _R0.item()
    else:
        _R0 = R0
    dsdt = (1-ss-ii)/L - _R0*ii*ss/D
    didt = _R0*ii*ss/D - ii/D

    return dsdt, didt

In [90]:

def rk4_climSIRSwy(ss, ii, t, dt, L, D, R0):
    '''Step change from classic Runge–Kutta (RK4) method for SIRS model of Baker et al., 2020.

    Input:
    -------
        ss: S/N
        ii: I/N
        t: time (units is 'day')
        dt: time step
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        tuple (dss, dii)
        dss: S/N step increase
        dii: I/N step increase
    '''
    model_climSIRS = model_climSIRSwy
    k1s, k1i = model_climSIRS(ss=ss, ii=ii, t=t, L=L, D=D, R0=R0)
    k2s, k2i = model_climSIRS(ss=ss+k1s*dt/2, ii=ii+k1i*dt/2, t=t+dt/2, L=L, D=D, R0=R0)
    k3s, k3i = model_climSIRS(ss=ss+k2s*dt/2, ii=ii+k2i*dt/2, t=t+dt/2, L=L, D=D, R0=R0)
    k4s, k4i = model_climSIRS(ss=ss+k3s*dt, ii=ii+k3i*dt, t=t+dt, L=L, D=D, R0=R0)

    dss = dt*(k1s + 2*k2s + 2*k3s + k4s)/6
    dii = dt*(k1i + 2*k2i + 2*k3i + k4i)/6

    return dss, dii

In [92]:

# def run_climSIRS(ii0=1e-5, ss0=None, T=5*365, dt=1, L=None, D=5, R0=2, dis='HKU1'):
def run_climSIRSwy(ii0=1e-5, ss0=None, T=1*365, dt=1, L=None, D=5, R0=2, dis='HKU1'):
    '''Integrate SIRS model from Baker et al., 2020.

    Input:
    -------
        I0: initial I/N
        T: time stop (start is 1; units is 'day')
        dt: time step
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        ds: solution for the xclimSIRS model.
    '''
    
    t_vec = np.arange(0, T, dt)
    
    if L is None and dis is not None:
        L = Ls[dis]
    if hasattr(R0, 'time'): # R0 time-dependent
        zeros = ii0*L*D*R0.isel(time=0, drop=True)*0 # broadcast array by R0.dims[1:]
    else: # R0 time-independent
        zeros = ii0*L*D*R0*0
    # default initial S/N value determined by initial I/N value
    if ss0 is None:
        ss0 = 1 - ii0
    # broadcast initial values by adding zeros
    ii = ii0 + zeros
    ss = ss0 + zeros
    
    ss_vec = [ss,]
    ii_vec = [ii,]
    print('time step (units: day): ', end='')
    for t in t_vec[1:]:
        if np.mod(t, 365) == 0:
            print(f'365x{t//365}', end='; ')
        dss, dii = rk4_climSIRSwy(ss=ss, ii=ii, t=t, dt=dt, L=L, D=D, R0=R0)
        ss = ss + dss
        ii = ii + dii
        ss_vec.append(ss)
        ii_vec.append(ii)
    print()
    if isinstance(ss, xr.DataArray):
        try:
            da_ss = xr.concat(ss_vec, dim=pd.Index(t_vec, name='t'))
            da_ii = xr.concat(ii_vec, dim=pd.Index(t_vec, name='t'))
        except MemoryError:
            ss_vec = [_ss.expand_dims('dum').chunk({'dum':1}).squeeze() for _ss in ss_vec]
            ii_vec = [_ii.expand_dims('dum').chunk({'dum':1}).squeeze() for _ii in ii_vec]
            da_ss = xr.concat(ss_vec, dim=pd.Index(t_vec, name='t'))
            da_ii = xr.concat(ii_vec, dim=pd.Index(t_vec, name='t'))
    else:
        da_ss = xr.DataArray(ss_vec, dims=['t'], coords=[t_vec])
        da_ii = xr.DataArray(ii_vec, dims=['t'], coords=[t_vec])
    da_ss.attrs['long_name'] = 'S/N'
    da_ii.attrs['long_name'] = 'I/N'
    ds = xr.Dataset(dict(ss=da_ss, ii=da_ii, year=da_ss.t/365))
    ds.t.attrs['units'] = 'day'
    

    return ds

ds = run_climSIRSwy()
# ds.ii.sel(latitude=latx, longitude=lonx, method='nearest').plot()
ds.ii.assign_coords(t=ds.year).rename(t='year').plot()
# da.sel(si='I/N').plot()

time step (units: day):

Out[92]:

[<matplotlib.lines.Line2D at 0x7f0898743650>]

In [140]:

dis = 'HKU1'
T = 15*365
ii0 = 1/8e6
q = get_q()
# q.plot()
R0 = q2R0(q, **climdepens[dis])
plt.figure()
# R0.plot()
# plt.figure()
ds = run_climSIRS(dis=dis, R0=R0, T=T, ii0=ii0)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot()
ds = run_climSIRSpy(dis=dis, R0=R0, T=T, ii0=ii0)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot()
ds = run_climSIRSwy(dis=dis, R0=R0, T=T, ii0=ii0)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot()

plt.legend(['odeint', 'solve_ivp', 'RK4'])

time step (units: day): 365x1; 365x2; 365x3; 365x4; 365x5; 365x6; 365x7; 365x8; 365x9; 365x10; 365x11; 365x12; 365x13; 365x14;

Out[140]:

<matplotlib.legend.Legend at 0x7f08bf938ad0>

In [109]:

dis = 'HKU1'
city = 'NYC'
latx, lonx = regions[city]
q = get_q('era5').sel(latitude=latx, longitude=lonx, method='nearest')
# q.plot()
R0 = q2R0(q, **climdepens[dis])
plt.figure()
# R0.plot()
ds = run_climSIRSwy(dis=dis, R0=R0, T=10*365)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot()
ds = run_climSIRS(dis=dis, R0=R0, T=10*365)
ds.ii.assign_coords(t=ds.year).rename(t='year').plot()

time step (units: day): 365x1; 365x2; 365x3; 365x4; 365x5; 365x6; 365x7; 365x8; 365x9;

Out[109]:

[<matplotlib.lines.Line2D at 0x7f08cd13f610>]

In [62]:

ifiles = '/tigress/wenchang/data/era5/analysis_wy/2m_specific_humidity/daily/clim/q2m.run7clim.2010-2019.nc'
ds = xr.open_dataset(ifiles)
q2m = ds.q2m.load().rename(dayofyear='time')
q2m

Out[62]:

xarray.DataArray

'q2m'

time: 366
latitude: 721
longitude: 1440

0.00047764648 0.00047764648 ... 0.00049648643 0.00049648643

array([[[0.00047765, 0.00047765, 0.00047765, ..., 0.00047765,
         0.00047765, 0.00047765],
        [0.00047971, 0.00047973, 0.00047975, ..., 0.0004796 ,
         0.00047963, 0.00047967],
        [0.00048245, 0.0004825 , 0.00048255, ..., 0.00048217,
         0.00048224, 0.00048233],
        ...,
        [0.00049061, 0.0004906 , 0.00049058, ..., 0.00049068,
         0.00049066, 0.00049064],
        [0.00048719, 0.00048718, 0.00048716, ..., 0.00048722,
         0.00048721, 0.0004872 ],
        [0.00047952, 0.00047952, 0.00047952, ..., 0.00047952,
         0.00047952, 0.00047952]],

       [[0.00048481, 0.00048481, 0.00048481, ..., 0.00048481,
         0.00048481, 0.00048481],
        [0.00048764, 0.00048765, 0.00048766, ..., 0.00048757,
         0.00048759, 0.00048761],
        [0.00049198, 0.00049201, 0.00049203, ..., 0.00049179,
         0.00049184, 0.0004919 ],
        ...,
        [0.00048531, 0.0004853 , 0.00048528, ..., 0.00048539,
         0.00048537, 0.00048534],
        [0.00048243, 0.00048242, 0.0004824 , ..., 0.00048246,
         0.00048245, 0.00048244],
        [0.0004755 , 0.0004755 , 0.0004755 , ..., 0.0004755 ,
         0.0004755 , 0.0004755 ]],

       [[0.00047514, 0.00047514, 0.00047514, ..., 0.00047514,
         0.00047514, 0.00047514],
        [0.00047989, 0.00047989, 0.00047989, ..., 0.00047985,
         0.00047986, 0.00047987],
        [0.0004862 , 0.00048621, 0.00048621, ..., 0.00048606,
         0.00048609, 0.00048613],
        ...,
        [0.00048336, 0.00048335, 0.00048333, ..., 0.00048343,
         0.00048341, 0.00048339],
        [0.00048104, 0.00048103, 0.00048102, ..., 0.00048107,
         0.00048106, 0.00048105],
        [0.00047475, 0.00047475, 0.00047475, ..., 0.00047475,
         0.00047475, 0.00047475]],

       ...,

       [[0.00049022, 0.00049022, 0.00049022, ..., 0.00049022,
         0.00049022, 0.00049022],
        [0.00049472, 0.00049475, 0.00049478, ..., 0.00049463,
         0.00049466, 0.00049469],
        [0.00049776, 0.00049782, 0.00049787, ..., 0.00049751,
         0.00049758, 0.00049765],
        ...,
        [0.00050307, 0.00050308, 0.00050307, ..., 0.00050312,
         0.00050311, 0.00050309],
        [0.00049993, 0.00049993, 0.00049992, ..., 0.00049995,
         0.00049994, 0.00049994],
        [0.00049273, 0.00049273, 0.00049273, ..., 0.00049273,
         0.00049273, 0.00049273]],

       [[0.00048698, 0.00048698, 0.00048698, ..., 0.00048698,
         0.00048698, 0.00048698],
        [0.00048996, 0.00048999, 0.00049002, ..., 0.00048986,
         0.00048989, 0.00048993],
        [0.00049317, 0.00049324, 0.0004933 , ..., 0.00049291,
         0.00049298, 0.00049307],
        ...,
        [0.0004991 , 0.0004991 , 0.00049909, ..., 0.00049915,
         0.00049914, 0.00049912],
        [0.0004961 , 0.00049609, 0.00049608, ..., 0.00049612,
         0.00049612, 0.00049611],
        [0.00048855, 0.00048855, 0.00048855, ..., 0.00048855,
         0.00048855, 0.00048855]],

       [[0.00067848, 0.00067848, 0.00067848, ..., 0.00067848,
         0.00067848, 0.00067848],
        [0.00070442, 0.00070443, 0.00070444, ..., 0.00070436,
         0.00070438, 0.0007044 ],
        [0.00072884, 0.00072889, 0.00072891, ..., 0.00072859,
         0.00072864, 0.00072871],
        ...,
        [0.00051024, 0.00051024, 0.00051024, ..., 0.00051025,
         0.00051025, 0.00051025],
        [0.00050557, 0.00050557, 0.00050557, ..., 0.00050557,
         0.00050557, 0.00050557],
        [0.00049649, 0.00049649, 0.00049649, ..., 0.00049649,
         0.00049649, 0.00049649]]], dtype=float32)

Coordinates: (3)
- latitude
  (latitude)
  float32
  90.0 89.75 89.5 ... -89.75 -90.0
  units :
  degrees_north
  long_name :
  latitude
```
array([ 90.  ,  89.75,  89.5 , ..., -89.5 , -89.75, -90.  ], dtype=float32)
```
- longitude
  (longitude)
  float32
  0.0 0.25 0.5 ... 359.5 359.75
  units :
  degrees_east
  long_name :
  longitude
```
array([0.0000e+00, 2.5000e-01, 5.0000e-01, ..., 3.5925e+02, 3.5950e+02,
       3.5975e+02], dtype=float32)
```
- time
  (time)
  int64
  1 2 3 4 5 6 ... 362 363 364 365 366
```
array([  1,   2,   3, ..., 364, 365, 366])
```
Attributes: (0)

In [64]:

# use ERA5 clim q
city = 'NYC'
latx, lonx = regions[city]
ds = run_climSIRSwy(L=66.25*7, D=5, I0=1e-5, T=10*365, dt=1, 
                    q=q2m.sel(latitude=latx, longitude=lonx, method='nearest'), 
                    R0=2)
# ds.ii.sel(latitude=latx, longitude=lonx, method='nearest').plot()
ds.ii.assign_coords(t=ds.year).rename(t='year').plot()

time step (units: day): 365x1; 365x2; 365x3; 365x4; 365x5; 365x6; 365x7; 365x8; 365x9; 365x10;

Out[64]:

[<matplotlib.lines.Line2D at 0x7f58fb16e310>]

In [65]:

ds

In [244]:

# sensitivity to time step
dss = []
for i,dt in enumerate([2, 1, 1/2]):
    ds = run_climSIRSwy(I0=1e-5, T=5*365, dt=dt, L=66.25*7, D=5, R0=None,
             q=q2m.sel(latitude=latx, longitude=lonx, method='nearest'), a=-227.5)
    ds.ii.plot(color=f'C{i}', label=f'dt = {dt}')
    dss.append(ds)
plt.legend()

time step (units: day): 365x1; 365x3; 365x5; 
time step (units: day): 365x1; 365x2; 365x3; 365x4; 365x5; 
time step (units: day): 365x1.0; 365x2.0; 365x3.0; 365x4.0; 365x5.0;

Out[244]:

<matplotlib.legend.Legend at 0x7ff4f07db290>

In [265]:

# sensitivity to time step
for i,dt in enumerate([2, 1, 1/2]):
#     ds = run_climSIRSwy(L=66.25*7, D=5, I0=1e-5, T=5*365, dt=dt, 
#              q=q2m.sel(latitude=latx, longitude=lonx, method='nearest'), R0=2)
    ds = dss[i]
    ds.ii.sel(t=slice(0,365)).plot(color=f'C{i}', label=f'dt = {dt}')
    dss.append(ds)
plt.legend()
plt.title('RK4', loc='left')

Out[265]:

Text(0.0, 1.0, 'RK4')

In [264]:

# sensitivity to time step
for i,dt in enumerate([2, 1, 1/2]):
#     ds = run_climSIRSwy(L=66.25*7, D=5, I0=1e-5, T=5*365, dt=dt, 
#              q=q2m.sel(latitude=latx, longitude=lonx, method='nearest'), R0=2)
    ds = dss[i]
    ds.ii.sel(t=slice(63,69)).plot(color=f'C{i}', label=f'dt = {dt}')
    dss.append(ds)
plt.legend()
plt.title('RK4', loc='left')

Out[264]:

Text(0.0, 1.0, 'RK4')

In [289]:

ds.ii.sel(latitude=30.5928, longitude=114.3055, method='nearest').plot()

Out[289]:

[<matplotlib.lines.Line2D at 0x7fe622b8ca58>]

In [282]:

ds.ii.sel(longitude=lonx, method='nearest').plot()

Out[282]:

[<matplotlib.lines.Line2D at 0x7fe6231c9c88>]

In [271]:

ds.ii.plot()
da.sel(si='I/N').plot()

Out[271]:

[<matplotlib.lines.Line2D at 0x7fe6232a4d30>]

In [297]:

ds = xr.open_dataset('/tigress/wenchang/analysis/climate_and_disease/climSIRS/HKU1.q2019.nc')
ds

Out[297]:

<xarray.Dataset>
Dimensions:    (latitude: 721, longitude: 1440, t: 3650)
Coordinates:
  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.25 359.5 359.75
  * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
Data variables:
    ss         (t, latitude, longitude) float32 ...
    ii         (t, latitude, longitude) float32 ...

In [303]:

ds2018 = xr.open_dataset('/tigress/wenchang/analysis/climate_and_disease/climSIRS/HKU1.q2018.nc')
ds2018

Out[303]:

<xarray.Dataset>
Dimensions:    (latitude: 721, longitude: 1440, t: 3650)
Coordinates:
  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.25 359.5 359.75
  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
  * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
Data variables:
    ss         (t, latitude, longitude) float32 ...
    ii         (t, latitude, longitude) float32 ...

In [307]:

ds.ii.sel(latitude=latx, longitude=lonx, method='nearest').assign_coords(t=ds.t.values/365).rename(t='year').plot()
ds2018.ii.sel(latitude=latx, longitude=lonx, method='nearest').assign_coords(t=ds.t.values/365).rename(t='year').plot()

plt.xlim(0, 4)

Out[307]:

(0, 4)

In [302]:

ds.sel(longitude=lonx, latitude=latx, method='nearest').to_dataframe().plot(x='ss', y='ii')

Out[302]:

<matplotlib.axes._subplots.AxesSubplot at 0x7fe621b2b240>

In [309]:

da = ds.ii.stack(grid=['latitude', 'longitude'])
da

Out[309]:

<xarray.DataArray 'ii' (t: 3650, grid: 1038240)>
array([[9.9999997e-06, 9.9999997e-06, 9.9999997e-06, ..., 9.9999997e-06,
        9.9999997e-06, 9.9999997e-06],
       [1.3074646e-05, 1.3074646e-05, 1.3074646e-05, ..., 1.3151950e-05,
        1.3151950e-05, 1.3151950e-05],
       [1.7212880e-05, 1.7212880e-05, 1.7212880e-05, ..., 1.7302013e-05,
        1.7302013e-05, 1.7302013e-05],
       ...,
       [7.2539071e-05, 7.2539071e-05, 7.2539071e-05, ..., 5.4199803e-03,
        5.4199803e-03, 5.4199803e-03],
       [7.7993653e-05, 7.7993653e-05, 7.7993653e-05, ..., 5.3909491e-03,
        5.3909491e-03, 5.3909491e-03],
       [8.4108389e-05, 8.4108389e-05, 8.4108389e-05, ..., 5.3762808e-03,
        5.3762808e-03, 5.3762808e-03]], dtype=float32)
Coordinates:
  * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
  * grid       (grid) MultiIndex
  - latitude   (grid) float64 90.0 90.0 90.0 90.0 ... -90.0 -90.0 -90.0 -90.0
  - longitude  (grid) float64 0.0 0.25 0.5 0.75 1.0 ... 359.0 359.2 359.5 359.8
Attributes:
    long_name:  I/N

spread¶

In [145]:

# NYC = (40.7128, -74.0060+360)
dis = 'OC43'
city = 'NYC'
latx, lonx = regions[city]
qyears = range(2010, 2020)
dss = []
for qyear in qyears:
    ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.{qyear}q.nc'
    ds = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')
    dss.append(ds)
ds = xr.concat(dss, dim=pd.Index(qyears, name='qyear'))

# ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.climq.nc'
# ds_clim = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')

ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.run7climq.nc'
ds_run7clim = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')

In [140]:

for qyear in qyears:
    ds.ii.sel(qyear=qyear).plot(color='C0', alpha=0.5, lw=0.5)
ds.ii.mean('qyear').plot(color='C0')
# ds_clim.ii.plot(color='C1')
ds_run7clim.ii.plot(color='C1')
plt.xlim(0, 365*2)
plt.title(f'{city} {dis}')

Out[140]:

Text(0.5, 1.0, 'NYC OC43')

In [147]:

# NYC = (40.7128, -74.0060+360)
dis = 'HKU1'
city = 'NYC'
latx, lonx = regions[city]
qyears = range(2010, 2020)
dss = []
for qyear in qyears:
    ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.{qyear}q.nc'
    ds = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')
    dss.append(ds)
ds = xr.concat(dss, dim=pd.Index(qyears, name='qyear'))

# ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.climq.nc'
# ds_clim = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')

ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.run7climq.nc'
ds_run7clim = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')

In [148]:

for qyear in qyears:
    ds.ii.sel(qyear=qyear).plot(color='C0', alpha=0.5, lw=0.5)
ds.ii.mean('qyear').plot(color='C0')
# ds_clim.ii.plot(color='C1')
ds_run7clim.ii.plot(color='C1')
plt.xlim(0, 365*3)
plt.title(f'{city} {dis}')

Out[148]:

Text(0.5, 1.0, 'NYC HKU1')

In [149]:

for qyear in qyears:
    ds.ii.sel(qyear=qyear).plot(color='C0', alpha=0.5, lw=0.5)
ds.ii.mean('qyear').plot(color='C0')
# ds_clim.ii.plot(color='C1')
ds_run7clim.ii.plot(color='C1')
plt.xlim(0, 365)
plt.title(f'{city} {dis}')

Out[149]:

Text(0.5, 1.0, 'NYC HKU1')

Regions¶

In [150]:

regions = dict(
    London=(51.5074, 0.1278),
    NYC=(40.7128, -74.0060+360),
    LA=(34.0522, -118.2437+360),
    Wuhan=(30.5928, 114.3055),
    Delhi=(28.7041, 77.1025),
    Maracaibo=(10.6427, -71.6125+360),
    Kinshasa=(-4.4419, 15.2663),
    Jakarta=(-6.2088, 106.8456),
    Johannesburg=(-26.2041, 28.0473),
    Buenos_Aires=(-34.6037, -58.3816+360),
    Melbourne=(-37.8136, 144.9631)
    
)
def do_region_sirs(region_name, dis='OC43', do_plot=False):
    latx, lonx = regions[region_name]
    qyears = range(2010, 2020)
    dss = []
    for qyear in qyears:
        ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.{qyear}q.nc'
        ds = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')
        dss.append(ds)
    ds = xr.concat(dss, dim=pd.Index(qyears, name='qyear'))
    
    ifile = f'/tigress/wenchang/analysis/climate_and_disease/climSIRS/{dis}.run7climq.nc'
    ds_clim = xr.open_dataset(ifile).sel(latitude=latx, longitude=lonx, method='nearest')
    
    # fig
    if do_plot:
        for qyear in qyears:
            ds.ii.sel(qyear=qyear).plot(color='C0', alpha=0.5, lw=0.5)
        ds.ii.mean('qyear').plot(color='C0', label='single-year-q ensemble')
        ds_clim.ii.plot(color='C1', label='climq')
        plt.xlim(0, 365*2)
        plt.title(f'{region_name}({latx:.1f}, {lonx:.1f}) {dis}')
        plt.legend()
    
    return ds, ds_clim

In [151]:

#
do_region_sirs(region_name='London', dis='OC43', do_plot=True)

Out[151]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 0.25
     latitude   float32 51.5
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999845 ... 0.41685742 0.4165671
     ii         (qyear, t) float32 1e-05 1.31814895e-05 ... 0.008281448,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 51.5
     longitude  float32 0.25
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [152]:

#
do_region_sirs(region_name='London', dis='HKU1', do_plot=True)

Out[152]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 0.25
     latitude   float32 51.5
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999858 ... 0.6685567 0.66908866
     ii         (qyear, t) float32 1e-05 1.1958771e-05 ... 0.00080141873,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 51.5
     longitude  float32 0.25
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [153]:

#
do_region_sirs(region_name='NYC', dis='OC43', do_plot=True)

Out[153]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 286.0
     latitude   float32 40.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999844 ... 0.376664 0.37640247
     ii         (qyear, t) float32 1e-05 1.3283939e-05 ... 0.009299738,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 40.75
     longitude  float32 286.0
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [154]:

#
do_region_sirs(region_name='NYC', dis='HKU1', do_plot=True)

Out[154]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 286.0
     latitude   float32 40.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.99998546 ... 0.66516936 0.6649421
     ii         (qyear, t) float32 1e-05 1.23389145e-05 ... 0.0040116212,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 40.75
     longitude  float32 286.0
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [155]:

#
do_region_sirs(region_name='LA', dis='OC43', do_plot=True)

Out[155]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 241.75
     latitude   float32 34.0
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999846 ... 0.412506 0.4123191
     ii         (qyear, t) float32 1e-05 1.3059735e-05 ... 0.0077307476,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 34.0
     longitude  float32 241.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [156]:

#
do_region_sirs(region_name='LA', dis='HKU1', do_plot=True)

Out[156]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 241.75
     latitude   float32 34.0
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.99998623 ... 0.6856362 0.68545973
     ii         (qyear, t) float32 1e-05 1.162454e-05 ... 0.0034230368,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 34.0
     longitude  float32 241.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [365]:

#
do_region_sirs(region_name='Wuhan', dis='OC43', do_plot=True)

Out[365]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 114.25
     latitude   float32 30.5
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.99998665 ... 0.8827582 0.8830247
     ii         (qyear, t) float32 1e-05 1.1219222e-05 ... 3.9486004e-06,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 30.5
     longitude  float32 114.25
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [366]:

#
do_region_sirs(region_name='Delhi', dis='OC43', do_plot=True)

Out[366]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 77.0
     latitude   float32 28.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999872 ... 0.895116 0.8946375
     ii         (qyear, t) float32 1e-05 1.0739206e-05 ... 0.0027466281,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 28.75
     longitude  float32 77.0
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [367]:

#
do_region_sirs(region_name='Maracaibo', dis='OC43', do_plot=True)

Out[367]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 288.5
     latitude   float32 10.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.99998766 ... 0.8850004 0.8850432
     ii         (qyear, t) float32 1e-05 1.0278916e-05 ... 0.0010715513,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 10.75
     longitude  float32 288.5
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [368]:

#
do_region_sirs(region_name='Kinshasa', dis='OC43', do_plot=True)

Out[368]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 15.25
     latitude   float32 -4.5
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999877 ... 0.8767824 0.87676954
     ii         (qyear, t) float32 1e-05 1.0252332e-05 ... 0.0014819218,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 -4.5
     longitude  float32 15.25
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [369]:

#
do_region_sirs(region_name='Jakarta', dis='OC43', do_plot=True)

Out[369]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 106.75
     latitude   float32 -6.25
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999877 ... 0.88639516 0.88637114
     ii         (qyear, t) float32 1e-05 1.0246806e-05 ... 0.0014117584,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 -6.25
     longitude  float32 106.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [370]:

#
do_region_sirs(region_name='Johannesburg', dis='OC43', do_plot=True)

Out[370]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 28.0
     latitude   float32 -26.25
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.99998766 ... 0.5247395 0.52582014
     ii         (qyear, t) float32 1e-05 1.0323673e-05 ... 3.4278935e-05,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 -26.25
     longitude  float32 28.0
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [371]:

#
do_region_sirs(region_name='Buenos_Aires', dis='OC43', do_plot=True)

Out[371]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 301.5
     latitude   float32 -34.5
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999876 ... 0.62869155 0.6295276
     ii         (qyear, t) float32 1e-05 1.0336555e-05 ... 7.819701e-05,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 -34.5
     longitude  float32 301.5
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

In [372]:

#
do_region_sirs(region_name='Melbourne', dis='OC43', do_plot=True)

Out[372]:

(<xarray.Dataset>
 Dimensions:    (qyear: 10, t: 3650)
 Coordinates:
     longitude  float32 145.0
     latitude   float32 -37.75
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
   * qyear      (qyear) int64 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
 Data variables:
     ss         (qyear, t) float32 0.99999 0.9999871 ... 0.6695969 0.66961473
     ii         (qyear, t) float32 1e-05 1.0805876e-05 ... 0.003868763,
 <xarray.Dataset>
 Dimensions:    (t: 3650)
 Coordinates:
     latitude   float32 -37.75
     longitude  float32 145.0
   * t          (t) int64 1 2 3 4 5 6 7 8 ... 3644 3645 3646 3647 3648 3649 3650
 Data variables:
     ss         (t) float32 ...
     ii         (t) float32 ...)

sensitivity to methods¶

In [107]:

#ERA5 q
ds_sp = xr.open_dataset('/tigress/wenchang/analysis/climate_and_disease/climSIRS/ClimSIRS.OC43.run7climq.cities.nc')
ds_wy = xr.open_dataset('/tigress/wenchang/analysis/climate_and_disease/climSIRS/xClimSIRS.OC43.run7climq.cities.nc')

In [126]:

cities = ds_sp.city.values
colors = [f'C{i}' for i in range(cities.size)]
peaks = []
peak_days = []
for city,c in zip(cities,colors):
    ds_sp.ii.sel(city=city).plot(color=c, label=f'{city} (odeint)')
    ds_wy.ii.sel(city=city).plot(color='k', ls=':', label=f'{city} (RK4)')
    da_sp = ds_sp.sel(city=city).ii.squeeze().isel(t=ds_sp.sel(city=city).ii.squeeze()==ds_sp.sel(city=city).ii.max())
    da_wy = ds_wy.sel(city=city).ii.squeeze().isel(t=ds_wy.sel(city=city).ii.squeeze()==ds_wy.sel(city=city).ii.max())
    peaks.append([da_sp.item(), da_wy.item()])
    peak_days.append([da_sp.t.astype('int').item(), da_wy.t.item()])
    plt.xlim(0, 365*2)
plt.legend()
plt.title('OC43')
df_peaks = pd.DataFrame(peaks, index=cities, columns=['odeint', 'RK4'])
print('I/N peak')
print(df_peaks)
print()
df_peak_days = pd.DataFrame(peak_days, index=cities, columns=['odeint', 'RK4'])
print('I/N peak day')
print(df_peak_days)

peaks
             odeint       RK4
NYC        0.218072  0.218362
Maracaibo  0.163862  0.163689
Delhi      0.190519  0.190240

peak days
           odeint  RK4
NYC            59   59
Maracaibo      75   76
Delhi          65   65

In [5]:

#MERRA2 q
ds_sp = xr.open_dataset('/tigress/wenchang/analysis/climate_and_disease/climSIRS/ClimSIRS.OC43.merra2wclimq.cities.nc')
ds_wy = xr.open_dataset('/tigress/wenchang/analysis/climate_and_disease/climSIRS/xClimSIRS.OC43.merra2wclimq.cities.nc')

In [7]:

#plot
dis = 'OC43'
qsource = 'MERRA2'
cities = ds_sp.city.values
colors = [f'C{i}' for i in range(cities.size)]
peaks = []
peak_days = []
for city,c in zip(cities,colors):
    ds_sp.ii.sel(city=city).plot(color=c, label=f'{city} (odeint)')
    ds_wy.ii.sel(city=city).plot(color='k', ls=':', label=f'{city} (RK4)')
    da_sp = ds_sp.sel(city=city).ii.squeeze().isel(t=ds_sp.sel(city=city).ii.squeeze()==ds_sp.sel(city=city).ii.max())
    da_wy = ds_wy.sel(city=city).ii.squeeze().isel(t=ds_wy.sel(city=city).ii.squeeze()==ds_wy.sel(city=city).ii.max())
    peaks.append([da_sp.item(), da_wy.item()])
    peak_days.append([da_sp.t.astype('int').item(), da_wy.t.item()])
    plt.xlim(0, 365*2)
plt.legend()
plt.title(f'{dis} {qsource}')
df_peaks = pd.DataFrame(peaks, index=cities, columns=['odeint', 'RK4'])
print('I/N peak')
print(df_peaks)
print()
df_peak_days = pd.DataFrame(peak_days, index=cities, columns=['odeint', 'RK4'])
print('I/N peak day')
print(df_peak_days)

I/N peak
             odeint       RK4
NYC        0.217644  0.217615
Maracaibo  0.169362  0.169373
Delhi      0.214861  0.214796

I/N peak day
           odeint  RK4
NYC            60   60
Maracaibo      75   75
Delhi          60   60

Use solve_ivp to build SIRS¶

In [12]:

ifiles = '/tigress/wenchang/data/era5/analysis_wy/2m_specific_humidity/daily/clim/q2m.run7clim.2010-2019.nc'
ds = xr.open_dataset(ifiles)
q2m = ds.q2m.load().rename(dayofyear='time')
q2m

Out[12]:

xarray.DataArray

'q2m'

time: 366
latitude: 721
longitude: 1440

0.00047764648 0.00047764648 ... 0.00049648643 0.00049648643

array([[[0.00047765, 0.00047765, 0.00047765, ..., 0.00047765,
         0.00047765, 0.00047765],
        [0.00047971, 0.00047973, 0.00047975, ..., 0.0004796 ,
         0.00047963, 0.00047967],
        [0.00048245, 0.0004825 , 0.00048255, ..., 0.00048217,
         0.00048224, 0.00048233],
        ...,
        [0.00049061, 0.0004906 , 0.00049058, ..., 0.00049068,
         0.00049066, 0.00049064],
        [0.00048719, 0.00048718, 0.00048716, ..., 0.00048722,
         0.00048721, 0.0004872 ],
        [0.00047952, 0.00047952, 0.00047952, ..., 0.00047952,
         0.00047952, 0.00047952]],

       [[0.00048481, 0.00048481, 0.00048481, ..., 0.00048481,
         0.00048481, 0.00048481],
        [0.00048764, 0.00048765, 0.00048766, ..., 0.00048757,
         0.00048759, 0.00048761],
        [0.00049198, 0.00049201, 0.00049203, ..., 0.00049179,
         0.00049184, 0.0004919 ],
        ...,
        [0.00048531, 0.0004853 , 0.00048528, ..., 0.00048539,
         0.00048537, 0.00048534],
        [0.00048243, 0.00048242, 0.0004824 , ..., 0.00048246,
         0.00048245, 0.00048244],
        [0.0004755 , 0.0004755 , 0.0004755 , ..., 0.0004755 ,
         0.0004755 , 0.0004755 ]],

       [[0.00047514, 0.00047514, 0.00047514, ..., 0.00047514,
         0.00047514, 0.00047514],
        [0.00047989, 0.00047989, 0.00047989, ..., 0.00047985,
         0.00047986, 0.00047987],
        [0.0004862 , 0.00048621, 0.00048621, ..., 0.00048606,
         0.00048609, 0.00048613],
        ...,
        [0.00048336, 0.00048335, 0.00048333, ..., 0.00048343,
         0.00048341, 0.00048339],
        [0.00048104, 0.00048103, 0.00048102, ..., 0.00048107,
         0.00048106, 0.00048105],
        [0.00047475, 0.00047475, 0.00047475, ..., 0.00047475,
         0.00047475, 0.00047475]],

       ...,

       [[0.00049022, 0.00049022, 0.00049022, ..., 0.00049022,
         0.00049022, 0.00049022],
        [0.00049472, 0.00049475, 0.00049478, ..., 0.00049463,
         0.00049466, 0.00049469],
        [0.00049776, 0.00049782, 0.00049787, ..., 0.00049751,
         0.00049758, 0.00049765],
        ...,
        [0.00050307, 0.00050308, 0.00050307, ..., 0.00050312,
         0.00050311, 0.00050309],
        [0.00049993, 0.00049993, 0.00049992, ..., 0.00049995,
         0.00049994, 0.00049994],
        [0.00049273, 0.00049273, 0.00049273, ..., 0.00049273,
         0.00049273, 0.00049273]],

       [[0.00048698, 0.00048698, 0.00048698, ..., 0.00048698,
         0.00048698, 0.00048698],
        [0.00048996, 0.00048999, 0.00049002, ..., 0.00048986,
         0.00048989, 0.00048993],
        [0.00049317, 0.00049324, 0.0004933 , ..., 0.00049291,
         0.00049298, 0.00049307],
        ...,
        [0.0004991 , 0.0004991 , 0.00049909, ..., 0.00049915,
         0.00049914, 0.00049912],
        [0.0004961 , 0.00049609, 0.00049608, ..., 0.00049612,
         0.00049612, 0.00049611],
        [0.00048855, 0.00048855, 0.00048855, ..., 0.00048855,
         0.00048855, 0.00048855]],

       [[0.00067848, 0.00067848, 0.00067848, ..., 0.00067848,
         0.00067848, 0.00067848],
        [0.00070442, 0.00070443, 0.00070444, ..., 0.00070436,
         0.00070438, 0.0007044 ],
        [0.00072884, 0.00072889, 0.00072891, ..., 0.00072859,
         0.00072864, 0.00072871],
        ...,
        [0.00051024, 0.00051024, 0.00051024, ..., 0.00051025,
         0.00051025, 0.00051025],
        [0.00050557, 0.00050557, 0.00050557, ..., 0.00050557,
         0.00050557, 0.00050557],
        [0.00049649, 0.00049649, 0.00049649, ..., 0.00049649,
         0.00049649, 0.00049649]]], dtype=float32)

Coordinates: (3)
- latitude
  (latitude)
  float32
  90.0 89.75 89.5 ... -89.75 -90.0
  units :
  degrees_north
  long_name :
  latitude
```
array([ 90.  ,  89.75,  89.5 , ..., -89.5 , -89.75, -90.  ], dtype=float32)
```
- longitude
  (longitude)
  float32
  0.0 0.25 0.5 ... 359.5 359.75
  units :
  degrees_east
  long_name :
  longitude
```
array([0.0000e+00, 2.5000e-01, 5.0000e-01, ..., 3.5925e+02, 3.5950e+02,
       3.5975e+02], dtype=float32)
```
- time
  (time)
  int64
  1 2 3 4 5 6 ... 362 363 364 365 366
```
array([  1,   2,   3, ..., 364, 365, 366])
```
Attributes: (0)

In [13]:

def func_R0(t, q, a=-227.5, R0_min=1.4, R0_max=2.5, cycle=True):
    '''R0 as function of t (time) and specific humidity (q).
    q has dims of (n_time, ...).'''
    #R0_min, R0_max = 1.4, 2.5
    #a = -227.5 for HKU1 and -32.5 for OC43
    if cycle:# cycle through the first 365 daily specific humidity
        dayofyear = int( np.mod(t-1, 365) + 1 )
        qt = q.isel(time=dayofyear-1, drop=True)
    else:
        qt = q.isel(time=int(t)-1, drop=True)
    R0t = (R0_max - R0_min) * np.exp(a*qt) + R0_min
    
    return R0t
# func_R0(365, q2m)

In [266]:

def model_climSIRSp(t, y, L, D, R0, q, a):
    '''SIRS model right hand side from Baker et al., 2020.
    Input:
    -------
        y: y[0:y.size//2] is S/N; y[y.size//2:] is I/N
        t: time (units is 'day')
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        dydy: half d(S/N)/dt, half d(I/N)/dt
    '''
    N = y.size
    ss = y[:N//2]
    ii = y[N//2:]
    if q is not None:
        R0 = func_R0(t=t, q=q, a=a, R0_min=1.4, R0_max=2.5)
        if R0.ndim < 1:
            R0 = R0.item()
    dydt = np.array([
        (1-ss-ii)/L - R0*ii*ss/D, 
        R0*ii*ss/D - ii/D
    ]).ravel()
    
    return dydt

In [61]:

from scipy.integrate import solve_ivp
def run_climSIRSp(I0=1e-5, T=5*365, dt=1, L=66.25*7, D=5, R0=2, q=None, a=-227.5):
    '''Integrate SIRS model from Baker et al., 2020.

    Input:
    -------
        I0: initial I/N
        T: time stop (start is 1; units is 'day')
        dt: time step
        L: immunity length, 7x66.25(62.5) days for HKU1 (OC43)
        D: infection period, 5 days
        R0 = beta*D, ignored if q is not None
        q: specific humidity, xr.DataArray, dims ('time', ...)
        a: climate dependence parameter, -227.5(-32.5) for HKU1(OC43)

    Return:
    --------
        ds: solution for the SIRS model.
    '''
    if q is not None and q.ndim > 1:
        grid = q.stack(grid=q.isel(time=0).dims).grid
        zeros = np.zeros(grid.size)
    else:
        zeros = I0*L*D*R0*0
    y0 = np.hstack([1.0-I0+zeros, I0+zeros])
    t = np.linspace(1, T, int(T/dt))
    sol = solve_ivp(model_climSIRS_RH,[1, T], y0, t_eval=t,  args=(L, D, R0, q, a))
    ss = sol.y[0:y0.size//2, :]
    ii = sol.y[y0.size//2:, :]
    da_ss = xr.DataArray(ss, dims=['grid', 't'],coords=[np.arange(1, y0.size//2+1), t])
    da_ss.attrs['long_name'] = 'S/N'
    da_ii = xr.DataArray(ii, dims=['grid', 't'],coords=[np.arange(1, y0.size//2+1), t])
    da_ii.attrs['long_name'] = 'I/N'
    
    if q is not None and q.ndim > 1:
        da_ss = da_ss.assign_coords(grid=grid).unstack()
        da_ii = da_ii.assign_coords(grid=grid).unstack()
    
    ds = xr.Dataset(dict(ss=da_ss, ii=da_ii))
    ds.t.attrs['units'] = 'day'
    
    return ds

latx, lonx = regions['NYC']
if lonx > 180:
    lonx -= 360
# q = q2m.sel(latitude=[latx,], longitude=[lonx,], method='nearest')
# R0 = np.arange(2, 1.3, -0.2)
# ds = pclimSIRS_int(R0=R0, q=None, T=2*365, I0=1/8e6, L=62.5*7, a=-32.5)
# ds.ii.rename(grid='R0').assign_coords(R0=R0).plot(hue='R0')
# ds

In [67]:

city = 'NYC'
latx, lonx = regions[city]
q = q2m.sel(latitude=[latx], longitude=[lonx], method='nearest')
ds = pclimSIRS_int(q=q, T=2*365, I0=1/8e6, L=62.5*7, a=-32.5)
ds.ii.plot()
ds

WY¶

In [6]:

q2m = (
    xr.open_dataset('/tigress/gvecchi/ANALYSIS/DISEASE/COVID-19/merra_weekly_clim.nc')
    .rename(WEEK_QV2M='q2m', TT='time', LAT='latitude', LON='longitude')['q2m']
    )
q2m

Out[6]:

<xarray.DataArray 'q2m' (time: 52, latitude: 361, longitude: 576)>
[10812672 values with dtype=float64]
Coordinates:
  * longitude  (longitude) float64 -180.0 -179.4 -178.8 ... 178.1 178.8 179.4
  * latitude   (latitude) float64 -90.0 -89.5 -89.0 -88.5 ... 89.0 89.5 90.0
  * time       (time) datetime64[ns] 1981-01-04 1981-01-11 ... 1981-12-27
Attributes:
    long_name:               Weekly-mean 2-meter specific humidity from MERRA2
    units:                   kg/kg
    history:                 From /tigress/samueltb/MERRA2/Daily/merged/QV2M/...
    climatology_time_range:  JAN-1980:AUG-2018
    long_name_mod:           regrid: 7 day on T@MOD

In [9]:

np.arange(1, 365, 7)

Out[9]:

array([  1,   8,  15,  22,  29,  36,  43,  50,  57,  64,  71,  78,  85,
        92,  99, 106, 113, 120, 127, 134, 141, 148, 155, 162, 169, 176,
       183, 190, 197, 204, 211, 218, 225, 232, 239, 246, 253, 260, 267,
       274, 281, 288, 295, 302, 309, 316, 323, 330, 337, 344, 351, 358])

In [18]:

q2m.assign_coords(time=q2m.time.dt.dayofyear)[:, 0, 0].plot()
q2m.assign_coords(time=q2m.time.dt.dayofyear)[:, 0, 0].reindex(time=np.arange(1, 366), method='nearest').plot()

Out[18]:

[<matplotlib.lines.Line2D at 0x7f296577fc18>]

In [19]:

q2m.assign_coords(time=q2m.time.dt.dayofyear).reindex(time=np.arange(1, 366), method='nearest')

Out[19]:

<xarray.DataArray 'q2m' (time: 365, latitude: 361, longitude: 576)>
[75896640 values with dtype=float64]
Coordinates:
  * time       (time) int64 1 2 3 4 5 6 7 8 ... 358 359 360 361 362 363 364 365
  * longitude  (longitude) float64 -180.0 -179.4 -178.8 ... 178.1 178.8 179.4
  * latitude   (latitude) float64 -90.0 -89.5 -89.0 -88.5 ... 89.0 89.5 90.0
Attributes:
    long_name:               Weekly-mean 2-meter specific humidity from MERRA2
    units:                   kg/kg
    history:                 From /tigress/samueltb/MERRA2/Daily/merged/QV2M/...
    climatology_time_range:  JAN-1980:AUG-2018
    long_name_mod:           regrid: 7 day on T@MOD

In [373]:

%%html
<script src="https://cdn.rawgit.com/parente/4c3e6936d0d7a46fd071/raw/65b816fb9bdd3c28b4ddf3af602bfd6015486383/code_toggle.js">
</script>

climSIRS¶

How to do ODE in Python?¶

Apply to the SIRS model¶

Non-dimensionalized y and t¶

Non-dimensionalized y only¶

sensitivity to R0¶

sensitivity to D¶

sensitivity to L¶

sensitivity to dt¶

SIRS model: vectorized version¶

climSIRS: SIRS coupled to climate¶

q and R0 annual cycle: idealized case¶

realistic q¶

climSIRS using odeint backend¶

idealized q¶

idealized q in a warming climate¶

climSIRS using solve_ivp backend¶

idealized q¶

idealized q: vectorized form¶

NYC q: impact of different years¶

Wuhan q: impact of different years¶

idealized q: qmin and qmax choice¶

xclimSIRS¶

spread¶

Regions¶

sensitivity to methods¶

Use solve_ivp to build SIRS¶

WY¶

climSIRS using `odeint` backend¶

climSIRS using `solve_ivp` backend¶