You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
314 KiB
314 KiB
Analysis of simulated IP grouped by model parameters and years¶
Import libraries¶
In [1]:
import datetime as dt
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
Helper functions, variables and classes¶
In [2]:
month_name = ["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"]
In [3]:
def std_error(avg_val, avg_sqr, counter):
"""
Estimate the standard error from the average value
and the average value of the square.
:param avg_val: the average value
:param avg_sqr: the average square value
:param counter: the size of the sample
:return: the standard error
"""
return np.sqrt((avg_sqr - avg_val**2) / (counter - 1))
Loading precalculated arrays¶
In [4]:
# numbers of simulated days for analysis
wrf_N_days = 4992
inm_N_days = 3650
In [5]:
# dates corresponding to the indices (0 axis) of the data arrays
# note: for WRF dates correspond to real dates
wrf_dt_indicies = np.array(
[dt.date(1980, 1, 1) + dt.timedelta(i * 3) for i in range(wrf_N_days)]
)
inm_dt_indicies = np.array(
[dt.date(2022, 1, 1) + dt.timedelta(i % 365) for i in range(inm_N_days)]
)
In [6]:
# dictionaries where processed data is saved. The dictionary keys represent the
# threshold value of CAPE - integer numbers from the list 500, 800, 1000, 1200,
# where the key 500 relates only to temperature-based modeling and is present
# only in dictionaries wrf_<...>, while keys 800, 1000, 1200 correspond to
# classical modeling using the WRF or INMCM model.
# dict contains arrays with dimensions (4992, 24)
# where axis 0 contains the day index (see `wrf_dt_indicies`)
wrf_hourly_total_ip = {
key: np.load(f"./data/WRF/WRF_HOURLY_TOTAL_IP_{parameters}.npy")[:wrf_N_days]
for key, parameters in zip([500, 800, 1000, 1200],
["500_T2_25", "800", "1000", "1200"])
}
# dict contains arrays with dimensions (3650, 24)
# where axis 0 contains the day index (see `inm_dt_indicies`)
# note: years in `inm_dt_indicies` is not real
inm_hourly_total_ip = {
key: np.load(f"./data/INMCM/INMCM_HOURLY_TOTAL_IP_{parameters}.npy")[:inm_N_days]
for key, parameters in zip([800, 1000, 1200],
["800", "1000", "1200"])
}
Figure 2.1¶
In [7]:
# calculate seasonal variation parameters for different parametrizations
# 7 sets for different axes of the figure
# each set with 12 values for each month
# monthly mean values of IP
data = np.zeros((7, 12))
# count of days per month
data_counter = np.zeros((7, 12), dtype=int)
# the sum of squares of IP daily values
data_sqr = np.zeros((7, 12))
for j, cape_thres in enumerate([800, 1000, 1200]):
for m in range(12):
# the indices `ax_idx` from 0 to 6 correspond to the axes:
# 0 - WRF, 1980-2020, CAPE=800,
# 1 - INMCM, 10 years, CAPE=800,
# 2 - WRF, 1980-2020, CAPE=1000,
# 3 - INMCM, 10 years, CAPE=1000,
# 4 - WRF, 1980-2020, CAPE=1200,
# 5 - INMCM, 10 years, CAPE=1200,
# 6 - Vostok station, 2006-2020
# calculate axes index for WRF axes (0, 2, 4)
ax_idx = j * 2
# filtering day indices belonging to a specific month
wrf_inds = [i for i, date in enumerate(wrf_dt_indicies)
if date.month == m + 1]
# slice IP for specific CAPE and day indices
ip = wrf_hourly_total_ip[cape_thres][wrf_inds]
# calculate seasonal variation parameters
data[ax_idx, m] = ip.mean()
data_counter[ax_idx, m] = len(ip)
data_sqr[ax_idx, m] = np.sum(ip.mean(axis=-1) ** 2) / len(ip)
# calculate axes index for INMCM axes (1, 3, 5)
ax_idx = j * 2 + 1
# filtering day indices belonging to a specific month
inmcm_inds =[i for i, date in enumerate(inm_dt_indicies)
if date.month == m + 1]
# slice IP for specific CAPE and day indices
ip = inm_hourly_total_ip[cape_thres][inmcm_inds]
# calculate seasonal variation parameters
data[ax_idx, m] = ip.mean()
data_counter[ax_idx, m] = len(ip)
data_sqr[ax_idx, m] = np.sum(ip.mean(axis=-1) ** 2) / len(ip)
# the last set is loaded from the processed data of Vostok station in script 2
# the data is loaded as a dictionary with keys `mean`, `counter`, `sqr`
# index 0: 2006-2020, index 1: 2006-2012, index 2: 2013-2020
vostok_results = np.load("./data/Vostok/vostok_2006_2020_results.npz")
data[-1] = vostok_results["mean"][0]
data_counter[-1] = vostok_results["counter"][0]
data_sqr[-1] = vostok_results["sqr"][0]
In [8]:
fig = plt.figure(figsize=(10, 14), constrained_layout=False)
ax = [None for _ in range(7)]
for n in range(6):
ax[n] = fig.add_subplot(4, 4, (2*n + 1, 2*n + 2))
ax[6] = fig.add_subplot(4, 4, (14, 15))
low = [200e3] * 6 + [100]
high = [280e3] * 6 + [180]
step = [20e3] * 6 + [20]
coeff = [1e3] * 6 + [1]
caption = ["WRF, 1980–2020, $\\varepsilon_0 = 0.8$ kJ/kg",
"INMCM, 10 years, $\\varepsilon_0 = 0.8$ kJ/kg",
"WRF, 1980–2020, $\\varepsilon_0 = 1$ kJ/kg",
"INMCM, 10 years, $\\varepsilon_0 = 1$ kJ/kg",
"WRF, 1980–2020, $\\varepsilon_0 = 1.2$ kJ/kg",
"INMCM, 10 years, $\\varepsilon_0 = 1.2$ kJ/kg",
"Vostok station, 2006–2020"]
col = ["royalblue"] * 6 + ["orangered"]
for n in range(7):
for axis in ["top", "bottom", "left", "right"]:
ax[n].spines[axis].set_linewidth(0.5)
ax[n].tick_params(length=6, width=0.5, axis="y")
ax[n].tick_params(length=0, width=0.5, axis="x")
ax[n].grid(color="0.", linewidth=0.5, axis="y")
ax[n].set_xlim((-0.5, 11.5))
ax[n].set_xticks(np.arange(12))
ax[n].set_xticklabels(month_name, fontsize="large", va="top")
ax[n].set_ylim((low[n], high[n]))
ax[n].set_yticks(np.arange(low[n], high[n] + step[n] / 2, step[n]))
ax[n].set_yticklabels((np.arange(low[n], high[n] + step[n] / 2,
step[n]) / coeff[n]).astype(int),
fontsize="large")
if n < 6:
ax[n].set_ylabel("Monthly mean\nionospheric potential, kV",
fontsize="large")
else:
ax[n].set_ylabel("Monthly mean fair-weather\npotential gradient, V/m",
fontsize="large")
ax[n].set_title(caption[n], fontsize="large")
ax[n].annotate("", xy=(12, np.min(data[n])), xycoords="data",
xytext=(12, np.max(data[n])), textcoords="data",
annotation_clip=False,
arrowprops=dict(
arrowstyle="<|-|>,head_length=0.8,head_width=0.3",
patchA=None, patchB=None, shrinkA=0., shrinkB=0.,
connectionstyle="arc3,rad=0.", fc="black",
linewidth=0.5
))
# ampl = (np.max(data[n]) - np.min(data[n])) / np.mean(data[n])
ampl = (np.max(data[n]) - np.min(data[n])) / \
np.sum(data[n] * data_counter[n]) * np.sum(data_counter[n])
ax[n].text(12.2, (np.min(data[n]) + np.max(data[n])) / 2,
f"{ampl * 100:.0f}%",
fontsize="large", ha="left", va="center", rotation=270)
fig.align_ylabels([ax[0], ax[2], ax[4]])
fig.align_ylabels([ax[1], ax[3], ax[5]])
for n in range(7):
ax[n].bar(np.arange(12), data[n],
yerr=std_error(data[n],
data_sqr[n],
data_counter[n]),
width=0.8, color=col[n])
for n in range(6):
ax[n].text(-0.3, 1.05, chr(ord("a") + 3 * (n % 2) + n // 2),
fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=ax[n].transAxes)
ax[6].text(-0.3, 1.05, chr(ord("a") + 6), fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=ax[6].transAxes)
fig.subplots_adjust(hspace=0.3, wspace=1.6)
for m in range(12):
ax[6].annotate(f"{data_counter[6, m]}",
xy=(m-0.15, ax[6].get_ylim()[0] + 3),
rotation=270, ha="center", va="bottom",
fontsize="large", color="0.")
fig.savefig("./figures_two_parts/ip_pg_total.eps", bbox_inches="tight")
Figure 2.5¶
In [9]:
# calculate seasonal variation parameters for different temporal parts
# 8 sets for different axes of the figure
# each set with 12 values for each month
# monthly mean values of IP
data = np.zeros((7, 12))
# count of days per month
data_counter = np.zeros((7, 12), dtype=int)
# the sum of squares of IP daily values
data_sqr = np.zeros((7, 12))
data = np.zeros((8, 12))
data_counter = np.zeros((8, 12), dtype=int)
data_sqr = np.zeros((8, 12))
# to construct this figure we divide the datasets into equal ranges of years
# the dictionary keys below denote the axis number on the figure
wrf_ranges = {
0: range(1981, 1990 + 1),
1: range(1991, 2000 + 1),
2: range(2001, 2010 + 1),
3: range(2011, 2020 + 1),
}
inm_ranges = {
4: range(0, 5),
5: range(5, 10)
}
for m in range(12):
for ax_idx in range(6):
if ax_idx in [0, 1, 2, 3]:
# filtering day indices belonging to a specific month
wrf_inds = [i for i, date in enumerate(wrf_dt_indicies)
if date.month == m + 1
and date.year in wrf_ranges[ax_idx]
]
# slice IP for CAPE = 1000 and aforementioned day indices
ip = wrf_hourly_total_ip[1000][wrf_inds]
# calculate seasonal variation parameters
data[ax_idx, m] = ip.mean()
data_counter[ax_idx, m] = len(ip)
data_sqr[ax_idx, m] = np.sum(ip.mean(axis=-1) ** 2) / len(ip)
if ax_idx in [4, 5]:
inmcm_inds = [i for i, date in enumerate(inm_dt_indicies)
if date.month == m + 1
and i//365 in inm_ranges[ax_idx]
]
# slice IP for CAPE = 1000 and aforementioned day indices
ip = inm_hourly_total_ip[1000][inmcm_inds]
# calculate seasonal variation parameters
data[ax_idx, m] = ip.mean()
data_counter[ax_idx, m] = len(ip)
data_sqr[ax_idx, m] = np.sum(ip.mean(axis=-1) ** 2) / len(ip)
# the last sets is loaded from the processed data of Vostok station in script 2
# the data is loaded as a dictionary with keys `mean`, `counter`, `sqr`
# index 0: 2006-2020, index 1: 2006-2012, index 2: 2013-2020
vostok_results = np.load("./data/Vostok/vostok_2006_2020_results.npz")
# part 2006-2012
data[6] = vostok_results["mean"][1]
data_counter[6] = vostok_results["counter"][1]
data_sqr[6] = vostok_results["sqr"][1]
# part 2013-2020
data[7] = vostok_results["mean"][2]
data_counter[7] = vostok_results["counter"][2]
data_sqr[7] = vostok_results["sqr"][2]
In [10]:
fig = plt.figure(figsize=(10, 14), constrained_layout=False)
ax = [None for _ in range(8)]
for n in range(8):
ax[n] = fig.add_subplot(4, 4, (2*n + 1, 2*n + 2))
low = [200e3] * 6 + [80] * 2
high = [280e3] * 6 + [180] * 2
step = [20e3] * 6 + [20] * 2
coeff = [1e3] * 6 + [1] * 2
caption = ["WRF, 1981–1990, $\\varepsilon_0 = 1$ kJ/kg",
"WRF, 1991–2000, $\\varepsilon_0 = 1$ kJ/kg",
"WRF, 2001–2010, $\\varepsilon_0 = 1$ kJ/kg",
"WRF, 2011–2020, $\\varepsilon_0 = 1$ kJ/kg",
"INMCM, 5 years (1–5), $\\varepsilon_0 = 1$ kJ/kg",
"INMCM, 5 years (6–10), $\\varepsilon_0 = 1$ kJ/kg",
"Vostok station, 2006–2012",
"Vostok station, 2013–2020"]
col = ["royalblue"] * 6 + ["orangered"] * 2
for n in range(8):
for axis in ["top", "bottom", "left", "right"]:
ax[n].spines[axis].set_linewidth(0.5)
ax[n].tick_params(length=6, width=0.5, axis="y")
ax[n].tick_params(length=0, width=0.5, axis="x")
ax[n].grid(color="0.", linewidth=0.5, axis="y")
ax[n].set_xlim((-0.5, 11.5))
ax[n].set_xticks(np.arange(12))
ax[n].set_xticklabels(month_name, fontsize="large", va="top")
ax[n].set_ylim((low[n], high[n]))
ax[n].set_yticks(np.arange(low[n], high[n] + step[n] / 2, step[n]))
ax[n].set_yticklabels((np.arange(low[n], high[n] + step[n] / 2,
step[n]) / coeff[n]).astype(int),
fontsize="large")
if n <= 5:
ax[n].set_ylabel("Monthly mean\nionospheric potential, kV",
fontsize="large")
else:
ax[n].set_ylabel("Monthly mean fair-weather\npotential gradient, V/m",
fontsize="large")
ax[n].set_title(caption[n], fontsize="large")
ax[n].annotate("", xy=(12, np.min(data[n])), xycoords="data",
xytext=(12, np.max(data[n])), textcoords="data",
annotation_clip=False,
arrowprops=dict(
arrowstyle="<|-|>,head_length=0.8,head_width=0.3",
patchA=None, patchB=None, shrinkA=0., shrinkB=0.,
connectionstyle="arc3,rad=0.", fc="black",
linewidth=0.5
))
# ampl = (np.max(data[n]) - np.min(data[n])) / np.mean(data[n])
ampl = (np.max(data[n]) - np.min(data[n])) / \
np.sum(data[n] * data_counter[n]) * np.sum(data_counter[n])
ax[n].text(12.2, (np.min(data[n]) + np.max(data[n])) / 2,
f"{ampl * 100:.0f}%",
fontsize="large", ha="left", va="center", rotation=270)
fig.align_ylabels([ax[0], ax[2], ax[4], ax[6]])
fig.align_ylabels([ax[1], ax[3], ax[5], ax[7]])
for n in range(8):
ax[n].bar(np.arange(12), data[n],
yerr=std_error(data[n],
data_sqr[n],
data_counter[n]),
width=0.8, color=col[n])
for n in range(8):
ax[n].text(-0.3, 1.05, chr(ord("a") + n), fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=ax[n].transAxes)
fig.subplots_adjust(hspace=0.3, wspace=1.6)
for n in range(6, 8):
for m in range(12):
ax[n].annotate(f"{data_counter[n, m]}",
xy=(m-0.15, ax[n].get_ylim()[0] + 3),
rotation=270, ha="center", va="bottom",
fontsize="large", color="0.")
fig.savefig("./figures_two_parts/ip_pg_partial.eps", bbox_inches="tight")
In [ ]: