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.
649 KiB
649 KiB
Seasonal variation of the simulated ionospheric potential contributions averaged over different latitudes¶
A decomposition into contributions, comparison with temperature patterns, a new parameterisation
The source code of Figures 2.3, 2.5 and 2.6
Importing libraries¶
In [1]:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from matplotlib import cm, colormaps, colors, transforms
import datetime as dt
import numpy as np
Helper functions, variables and classes¶
In [2]:
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))
In [3]:
# a helper function to show the peak-to-peak amplitude
def draw_arrows(ax=None, miny=0, maxy=1, exty=1,
arrow_pos=12, text_pos=12.2, ampl=0, inward=False):
arrow_props = {
"patchA": None, "patchB": None,
"shrinkA": 0.0, "shrinkB": 0.0,
"connectionstyle": "arc3,rad=0.", "fc": "black",
"linewidth": 0.5
}
coords = {
"xycoords": "data",
"textcoords": "data",
"annotation_clip": False
}
if inward:
ax.annotate("", xy=(arrow_pos, miny), xytext=(arrow_pos, maxy),
arrowprops={"arrowstyle": "-", **arrow_props}, **coords)
for y, yt in zip([maxy + exty, miny - exty], [maxy, miny]):
ax.annotate("", xy=(arrow_pos, y), xytext=(arrow_pos, yt),
arrowprops={
"arrowstyle": "<|-,head_length=0.8,head_width=0.3",
**arrow_props
}, **coords)
else:
ax.annotate("", xy=(arrow_pos, miny), xytext=(arrow_pos, maxy),
arrowprops={
"arrowstyle": "<|-|>,head_length=0.8,head_width=0.3",
**arrow_props
}, **coords)
ax.text(
text_pos, (miny + maxy) / 2,
f"{ampl * 100:.0f}%",
fontsize="large", ha="left", va="center", rotation=270
)
In [4]:
month_name = ["J", "F", "M", "A", "M", "J",
"J", "A", "S", "O", "N", "D"]
month_name_3 = ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
In [5]:
# area factors for different latitudes
area_factor = (
np.cos(np.arange(180) * np.pi / 180) -
np.cos(np.arange(1, 181) * np.pi / 180)
) / 2
In [6]:
class LatitudeBand():
"""
Represents a latitude band object with specific
latitude range and model information
Attributes:
slice (slice): slice for the latitude range
model (str): model associated with the latitude band
color (str): color associated with the latitude band
"""
def __init__(self, expression, model=None):
"""
Initializes a LatitudeBand object with the given latitude
expression and model.
Args:
expression (str): latitude range expression like "0-8S", "8N-8S"
model (str, optional): "WRF" or "INMCM", defaults to None
"""
# splitting the expression into two parts
# and replacing "0" with "0S" if necessary
s1, s2 = [s if s != "0" else "0S" for s in expression.split("-")]
assert all(s[-1] in "NS" and s[:-1].isdigit() for s in [s1, s2])
# calculating the boundaries of the latitude range
# taking into account the direction
N1 = int(s1[:-1]) * (-1 if s1[-1] == "S" else 1)
N2 = int(s2[:-1]) * (-1 if s2[-1] == "S" else 1)
# formatting latitude values for display using thin spaces (" ")
p1 = "0°" if N1 == 0 else f"{abs(N1):d}° {'S' if N1 < 0 else 'N'}"
p2 = "0°" if N2 == 0 else f"{abs(N2):d}° {'S' if N2 < 0 else 'N'}"
# set the colour for each latitude band
Z = (N1 + N2) / 2
cmap = colormaps.get_cmap("tab20")
if 0 < Z < 9:
self.color = cmap(7 / 40)
elif -9 < Z < 0:
self.color = cmap(5 / 40)
elif -18 < Z < -9:
self.color = cmap(39 / 40)
elif -90 < Z < -18:
self.color = cmap(37 / 40)
elif 9 < Z < 18:
self.color = cmap(35 / 40)
elif 18 < Z < 90:
self.color = cmap(33 / 40)
# ensuring that the model is valid, setting the step in degrees
# and the central index based on the model
assert model in ["WRF", "INMCM"]
degrees = 1 if model == "WRF" else 1.5
central_idx = 90 if model == "WRF" else 60
# checking that latitude values are divisible
# by the appropriate degrees
assert all(N % degrees == 0 for N in [N1, N2])
# converting latitude values to indices and sorting them
N1, N2 = sorted([int(N // degrees + central_idx) for N in [N1, N2]])
# creating a slice for the latitude range
self.slice = slice(N1, N2, 1)
# saving the latitude band information
self.model = model
self.title = f"{p1}–{p2}" # like ‘9° S–18° S’
self.label = self.title.replace(" ", "")
def __repr__(self):
return self.label
Loading precalculated arrays¶
In [7]:
# the number of simulated days used for analysis
wrf_N_days = 4992
inm_N_days = 41 * 365
In [8]:
# dates corresponding to the indices (0 axis) of the data arrays
# note: in the case of the WRF the dates correspond to real dates
# note: in the case of the INMCM the are 41 365-day years
wrf_dt_indices = np.array(
[dt.date(1980, 1, 1) + dt.timedelta(i * 3) for i in range(wrf_N_days)]
)
inm_dt_indices = np.array(
[dt.date(2022, 1, 1) + dt.timedelta(i % 365) for i in range(inm_N_days)]
)
In [9]:
# dictionaries where the processed data are saved
# dictionary keys represent CAPE threshold values in J/kg
# 500 J/kg corresponds to WRF simulations with a temperature threshold of 25 °C
# 600, 800, 1000, 1200 J/kg correspond to usual WRF and INMCM simulations
# monthly average air temperature at the height of 2 m for different latitudes
# with the shape (180, 12)
wrf_LATxMON_t2 = np.load("./data/WRF/WRF_T2_LATxMON.npy")
# monthly average WRF-based IP contributions for different latitudes
# with the shape (180, 12)
wrf_LATxMON_ip = {
key: np.load(f"./data/WRF/WRF_IP_{parameters}_LATxMON.npy")
for key, parameters in zip([500, 600, 800, 1000, 1200],
["500_T2_25", "600", "800", "1000", "1200"])
}
# monthly average WRF-based INMCM contributions for different latitudes
# with the shape (180, 12)
inm_LATxMON_ip = {
key: np.load(f"./data/INMCM/INMCM_IP_{parameters}_LATxMON.npy")
for key, parameters in zip([600, 800, 1000, 1200],
["600", "800", "1000", "1200"])
}
# a dictionary containing arrays with the dimensions (4992, 24)
# axis 0 indicates the index of the day (see `wrf_dt_indices`)
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, 600, 800, 1000, 1200],
["500_T2_25", "600", "800", "1000", "1200"])
}
Auxiliary calculations of some values mentioned in the text¶
In [10]:
# the number of days for each month
wrf_numdays = np.load("./data/WRF/WRF_NUMDAYS_MON.npy")
inm_numdays = np.load("./data/INMCM/INMCM_NUMDAYS_MON.npy")
# to calculate the annual mean value, we use
# `(wrf_data_LATxMON[:, :].sum(axis=0) * days).sum() / days.sum()`
# rather than
# `wrf_data_LATxMON[:, :].sum(axis=0).mean()`,
# since
# `((a1+a2+a3)/3 + (b1+b2)/2)/2 != (a1+a2+a3+b1+b2)/5`
In [11]:
# calculating average contributions of the regions to
# the IP (as fractions of the total IP)
# for specific latitude ranges and parameterisations
# selecting the latitude ranges for the analysis
band_names = ["9S-9N", "18S-18N", "30S-30N"]
wrf_bands = [LatitudeBand(t, model="WRF") for t in band_names]
inm_bands = [LatitudeBand(t, model="INMCM") for t in band_names]
print("Model | CAPE | T2, °C | Band | % ")
# iterating over the latitude ranges selected above
for i, (wrf_band, inm_band) in enumerate(zip(wrf_bands, inm_bands)):
print("-" * 42)
# we will analyse all parameterisations only for the 9N-9S range;
# for the other ranges we will consider only the parameterization
# with CAPE = 1000 J/kg
capes = [1000]
if i == 0:
capes = [500, 600, 800, 1000, 1200]
# iterating over different parameterisations
for cape in capes:
# WRF
# calculating the seasonal variation of IP contributions
# from each range -> array (12,)
wrf_ip_MON = wrf_LATxMON_ip[cape][wrf_band.slice].sum(axis=0)
# calculating the contribution of each range
# to the total IP as a percentage of the mean -> single value
wrf_portion = np.average(wrf_ip_MON, weights=wrf_numdays) / 240e3 * 100
# note: The averaging over months takes into account the number of
# days taken from the original series for each month, as a direct
# averaging may lead to a slight error
p = "25 " if cape == 500 else "- "
print(f"WRF | {cape:4d} | {p} |{wrf_band.label:>10} |"
f" {wrf_portion:>4.2f}%")
# INMCM
if cape != 500:
# calculating the seasonal variation of IP contributions
# from each range -> array (12,)
inm_ip_MIN = inm_LATxMON_ip[cape][inm_band.slice].sum(axis=0)
# calculating the contribution of each range
# to the total IP as a percentage of the mean -> single value
inm_portion = np.average(inm_ip_MIN, weights=inm_numdays) / 240e3 * 100
print(f"INMCM | {cape:4d} | - |{inm_band.label:>10} |"
f" {inm_portion:>4.2f}%")
In [12]:
# calculating
# (i) the positions of the extrema of the seasonal variation of
# contributions to the IP from individual latitude ranges,
# (ii) the peak-to-peak amplitudes of the seasonal variation of
# contributions to the IP from individual latitude ranges
# for different parameterisations
# selecting the latitude ranges for the analysis
band_names = ["90S-9S", "9S-9N", "9N-90N"]
wrf_bands = [LatitudeBand(t, model="WRF") for t in band_names]
inmcm_bands = [LatitudeBand(t, model="INMCM") for t in band_names]
print("Model | CAPE | T2,°C | Band | Min | Max | Peak-to-peak, %")
# iterating over the latitude ranges selected above
for i, (wrf_band, inm_band) in enumerate(zip(wrf_bands, inmcm_bands)):
print("-" * 61)
for cape in [500, 600, 800, 1000, 1200]:
#WRF
# calculating the seasonal variation of IP contributions
# from each range -> array (12,)
seas_var = wrf_LATxMON_ip[cape][wrf_band.slice].sum(axis=0)
# finding the positions of the extrema
month_min = month_name_3[np.argmin(seas_var)]
month_max = month_name_3[np.argmax(seas_var)]
# computing the peak-to-peak amplitude of the seasonal variation
pk_pk_ampl = (seas_var.max() - seas_var.min())/seas_var.mean() * 100
p = "25 " if cape == 500 else "- "
print(f"WRF | {cape:4} | {p:<5} |{inm_band.label:>9} |"
f" {month_min} | {month_max} | {pk_pk_ampl:>6.2f}%")
# INMCM
if cape != 500:
# calculating the seasonal variation of IP contributions
# from each range -> array (12,)
seas_var = inm_LATxMON_ip[cape][inm_band.slice].sum(axis=0)
# finding the positions of the extrema
month_min = month_name_3[np.argmin(seas_var)]
month_max = month_name_3[np.argmax(seas_var)]
# computing the peak-to-peak amplitude of the seasonal variation
pk_pk_ampl = (seas_var.max() - seas_var.min())/seas_var.mean() * 100
print(f"INMCM | {cape:4} | - |{inm_band.label:>9} |"
f" {month_min} | {month_max} | {pk_pk_ampl:>6.2f}%")
In [13]:
# calculating maximum values of contributions to the IP
# from different latitude ranges; this calculation is needed to
# manually select the colour bar boundaries
# selecting the latitude ranges for the analysis
band_names = ["18N-90N", "9N-18N", "0-9N", "0-9S", "9S-18S", "18S-90S"]
wrf_bands = [LatitudeBand(t, model="WRF") for t in band_names]
inmcm_bands = [LatitudeBand(t, model="INMCM") for t in band_names]
print("CAPE = 1000 J/kg\n")
print("Model | Band | Max IP contrib.")
# iterating over the latitude ranges selected above
for i, (wrf_band, inm_band) in enumerate(zip(wrf_bands, inmcm_bands)):
print("-" * 35)
# calculating maximum values (over months), V
wrf_max = wrf_LATxMON_ip[1000][wrf_band.slice].sum(axis=0).max()
inm_max = inm_LATxMON_ip[1000][inm_band.slice].sum(axis=0).max()
print(f"WRF | {wrf_band.label:>9} | {wrf_max:>9.2f} V")
print(f"INMCM | {inm_band.label:>9} | {inm_max:>9.2f} V")
In [14]:
# calculating maximum and minimum values of the average
# air temperature at the height of 2 m in each latitude range;
# this calculation is needed to manually select the colour bar boundaries
# selecting the latitude ranges for the analysis
band_names = ["18N-30N", "9N-18N", "0-9N", "0-9S", "9S-18S", "18S-30S"]
wrf_bands = [LatitudeBand(t, model="WRF") for t in band_names]
print("Model | Band | Min T2 | Max T2")
print("-" * 37)
# iterating over the latitude ranges selected above
for wrf_band in wrf_bands:
# when averaging, we take the area factor into account
tmp_t2 = np.average(wrf_LATxMON_t2[wrf_band.slice],
axis=0,
weights=area_factor[wrf_band.slice])
print(f"WRF | {wrf_band.label:>9} | {tmp_t2.min():>3.2f}°C | {tmp_t2.max():>3.2f}°C")
Figure 2.5¶
Decomposition of the simulated IP variation into contributions
In [15]:
# selecting the latitude ranges for the analysis
band_names = ["0-9N", "0-9S", "9S-18S", "18S-90S", "9N-18N", "18N-90N"]
wrf_bands = [LatitudeBand(t, model="WRF") for t in band_names]
inm_bands = [LatitudeBand(t, model="INMCM") for t in band_names]
bar_data = np.zeros((8, len(band_names), 12))
# (subplots, strips, months)
for j, cape in enumerate([600, 800, 1000, 1200]):
ax_idx = j * 2
for i, band in enumerate(wrf_bands):
bar_data[ax_idx, i] = wrf_LATxMON_ip[cape][band.slice].sum(axis=0)
ax_idx = j * 2 + 1
for i, band in enumerate(inm_bands):
bar_data[ax_idx, i] = inm_LATxMON_ip[cape][band.slice].sum(axis=0)
In [16]:
fig = plt.figure(figsize=(10, 16), constrained_layout=False)
ax = [None for _ in range(8)]
for n in range(8):
ax[n] = fig.add_subplot(5, 4, (2*n + 1, 2*n + 2))
low = [0e3] * 8
high = [300e3] * 8
step = [50e3] * 8
coeff = [1e3] * 8
caption = ["WRF, 1980–2020, $\\varepsilon_0 = 0.6$ kJ/kg",
"INMCM, 41 years, $\\varepsilon_0 = 0.6$ kJ/kg",
"WRF, 1980–2020, $\\varepsilon_0 = 0.8$ kJ/kg",
"INMCM, 41 years, $\\varepsilon_0 = 0.8$ kJ/kg",
"WRF, 1980–2020, $\\varepsilon_0 = 1$ kJ/kg",
"INMCM, 41 years, $\\varepsilon_0 = 1$ kJ/kg",
"WRF, 1980–2020, $\\varepsilon_0 = 1.2$ kJ/kg",
"INMCM, 41 years, $\\varepsilon_0 = 1.2$ kJ/kg"]
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")
ax[n].set_ylabel("Monthly mean\nionospheric potential, kV",
fontsize="large")
ax[n].set_title(caption[n], fontsize="large")
fig.align_ylabels([ax[0], ax[2], ax[4]])
fig.align_ylabels([ax[1], ax[3], ax[5]])
for n in range(8):
bottom_data = np.zeros(12)
for k in range(len(band_names)):
data = bar_data[n, k]
ax[n].bar(np.arange(12), data, width=0.8,
color=wrf_bands[k].color,
bottom=bottom_data, label=wrf_bands[k].title,
)
bottom_data += data
for n in range(8):
ax[n].text(-0.28, 1.05, chr(ord("a") + 4 * (n % 2) + n // 2),
fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=ax[n].transAxes)
fig.subplots_adjust(hspace=0.3, wspace=1.)
leg_width_rel = (ax[7].get_position().x1 - ax[6].get_position().x0) \
/ (ax[6].get_position().x1 - ax[6].get_position().x0)
leg = ax[6].legend(bbox_to_anchor=(0., -1./3, leg_width_rel, 1./6), ncols=6,
loc="lower center", borderaxespad=0.,
mode="expand", fontsize="large",
framealpha=1, edgecolor="0.",
handlelength=1., handletextpad=0.5)
leg.get_frame().set_linewidth(0.5)
fig.savefig("figures/ip_decomposed.eps", bbox_inches="tight")
Figure 2.3¶
Comparison of the seasonal patterns in IP contributions and temperatures
In [17]:
# selecting the latitude ranges for the analysis
band_names = ["18N-90N", "9N-18N", "0-9N", "0-9S", "9S-18S", "18S-90S"]
wrf_bands = [LatitudeBand(t, model="WRF") for t in band_names]
inm_bands = [LatitudeBand(t, model="INMCM") for t in band_names]
wrf_t2_bands = [LatitudeBand(t, model="WRF")
for t in ["18N-30N"] + band_names[1:-1] + ["18S-30S"]]
In [18]:
wrf_slice = LatitudeBand("30S-30N", model="WRF").slice
inm_slice = LatitudeBand("30S-30N", model="INMCM").slice
WRF_IP_LxM = wrf_LATxMON_ip[1000][wrf_slice] / 360
INM_IP_LxM = inm_LATxMON_ip[1000][inm_slice] / 180
T2_LxM = wrf_LATxMON_t2[wrf_slice]
In [19]:
strip_map = np.full((180, 360, 4), 1.)
for band in wrf_t2_bands:
for i in range(band.slice.start, band.slice.stop, band.slice.step):
for j in range(360):
strip_map[i, j] = band.color
In [20]:
wrf_ip_sum_over_bands = np.zeros((6, 12))
inm_ip_sum_over_bands = np.zeros((6, 12))
wrf_t2_avg_over_bands = np.zeros((6, 12))
for j in range(6):
wrf_ip_sum_over_bands[j] = wrf_LATxMON_ip[1000][wrf_bands[j].slice].sum(axis=0)
inm_ip_sum_over_bands[j] = inm_LATxMON_ip[1000][inm_bands[j].slice].sum(axis=0)
wrf_t2_avg_over_bands[j] = np.sum(
wrf_LATxMON_t2[wrf_t2_bands[j].slice]
* area_factor[wrf_t2_bands[j].slice, np.newaxis],
axis=0,
) / np.sum(area_factor[wrf_t2_bands[j].slice, np.newaxis])
In [21]:
fig_width = 10
fig_height = 12
width = 3.
wsp = 1.
tot_width = 3 * width + 4 * wsp
widthM = 0.6 * tot_width
height = np.array([2., 2., 3., 3., 2., 1.])
hsp = 0.8
hsep1 = 2.
hsep2 = 3.5
assert len(height) == 6
heightD = 4.
heightM = widthM/tot_width * 60/360 * fig_width/fig_height
tot_height = (sum(height) + heightD + 8 * hsp + hsep1 + hsep2) / (1 - heightM)
heightM *= tot_height
sizeT = np.array([6., 3., 2., 2., 3., 4.])
hspT = hsp
assert len(sizeT) == 6
heightT = sizeT * (sum(height) + 5 * hsp - 5 * hspT) / sum(sizeT)
minT = np.array([16., 24., 24., 24., 22., 18.])
fig = plt.figure(figsize=(fig_width, fig_height), constrained_layout=False)
ax = [None for _ in range(12)]
axT = [None for _ in range(6)]
for n in range(6):
ax[2*n] = fig.add_axes(
(
(width + 2 * wsp) / tot_width,
(sum(height[n+1:]) + heightM + (6 - n) * hsp + hsep1) / tot_height,
width / tot_width,
height[n] / tot_height
)
)
ax[2*n + 1] = fig.add_axes(
(
(2 * width + 3 * wsp) / tot_width,
(sum(height[n+1:]) + heightM + (6 - n) * hsp + hsep1) / tot_height,
width / tot_width,
height[n] / tot_height
)
)
axT[n] = fig.add_axes(
(
wsp / tot_width,
(sum(heightT[n+1:]) + (5 - n) * hspT
+ heightM + hsp + hsep1) / tot_height,
width / tot_width,
heightT[n] / tot_height
)
)
axD = [None for _ in range(3)]
for n in range(3):
axD[n] = fig.add_axes(
(
(n * width + (n + 1) * wsp) / tot_width,
(sum(height) + heightM + 7 * hsp + hsep1 + hsep2) / tot_height,
width / tot_width,
heightD / tot_height
)
)
axM = fig.add_axes(
(
(1 - widthM / tot_width) / 2,
hsp / tot_height,
widthM / tot_width,
heightM / tot_height
),
projection=ccrs.PlateCarree(central_longitude=0)
)
axM.coastlines("110m", linewidth=0.5)
caption = ["WRF, 1980–2020, $\\varepsilon_0 = 1$ kJ/kg",
"INMCM, 41 years, $\\varepsilon_0 = 1$ kJ/kg"]
captionT = "WRF, 1980–2020"
captionD = ["WRF, 1980–2020",
"WRF, 1980–2020, $\\varepsilon_0 = 1$ kJ/kg",
"INMCM, 41 years, $\\varepsilon_0 = 1$ kJ/kg"]
for n in range(12):
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")
low = 0
high = 50e3 * height[n // 2]
step = 50e3
coeff = 1e3
ax[n].set_ylim((low, high))
ax[n].set_yticks(np.arange(low, high + step / 2, step))
ax[n].set_yticklabels((np.arange(low, high + step / 1,
step) / coeff).astype(int),
fontsize="large")
if n < 2:
ax[n].set_title(caption[n], fontsize="large")
for n in range(2):
lbl = ax[n].set_ylabel("Monthly mean contributions "
"to the ionospheric potential, kV",
fontsize="large", labelpad=0.)
lbl.set_position((0., 1 - (sum(height) + 5 * hsp) / 2 / height[0]))
for n in range(6):
for axis in ["top", "bottom", "left", "right"]:
axT[n].spines[axis].set_linewidth(0.5)
axT[n].tick_params(length=6, width=0.5, axis="y")
axT[n].tick_params(length=0, width=0.5, axis="x")
axT[n].grid(color="0.", linewidth=0.5, axis="y")
axT[n].set_xlim((-0.5, 11.5))
axT[n].set_xticks(np.arange(12))
axT[n].set_xticklabels(month_name, fontsize="large", va="top")
low = minT[n]
step = 2
high = minT[n] + step * sizeT[n]
coeff = 1
axT[n].set_ylim((low, high))
axT[n].set_yticks(np.arange(low, high + step / 2, step))
axT[n].set_yticklabels((np.arange(low, high + step / 2,
step) / coeff).astype(int),
fontsize="large")
if n == 0:
axT[n].set_title(captionT, fontsize="large")
lbl = axT[0].set_ylabel("Monthly mean surface air temperature, °C",
fontsize="large")
lbl.set_position((0., 1 - (sum(heightT) + 5 * hspT) / 2 / heightT[0]))
for n in range(3):
for axis in ["top", "bottom", "left", "right"]:
axD[n].spines[axis].set_linewidth(0.5)
axD[n].tick_params(length=6, width=0.5)
axD[n].tick_params(length=0, width=0.5, which="minor")
axD[n].grid(color="0.", linewidth=0.5, axis="y")
axD[n].set_xlim((0, 12))
axD[n].set_xticks(np.arange(0, 12.5, 1))
axD[n].set_xticks(np.arange(0.5, 12, 1), minor=True)
axD[n].set_xticklabels([])
axD[n].set_xticklabels(month_name, fontsize="large", va="top",
minor=True)
axD[n].set_ylim((-30, 30))
axD[n].set_yticks(np.arange(-30, 31, 10))
axD[n].set_yticklabels(["0" if v == 0 else f"{-v}S" if v < 0 else f"{v}N"
for v in range(-30, 31, 10)],
fontsize="large")
axD[n].set_title(captionD[n], fontsize="large")
cax = [None for _ in range(3)]
for n in range(3):
cax[n] = fig.add_axes(
(
axD[n].get_position().x0,
axD[n].get_position().y0 - 0.035,
axD[n].get_position().x1 - axD[n].get_position().x0,
0.01
)
)
min_val = [10, 0, 0]
max_val = [30, 50, 120]
step_val = [5, 10, 20]
cbar_cmap = ["inferno", "plasma_r", "plasma_r"]
cbar_norm = [None] * 3
for n in range(3):
cbar_norm[n] = colors.Normalize(vmin=min_val[n], vmax=max_val[n])
cbar = fig.colorbar(cm.ScalarMappable(norm=cbar_norm[n],
cmap=cbar_cmap[n]),
cax=cax[n], orientation="horizontal")
cbar.outline.set_linewidth(0.5)
cbar.ax.tick_params(length=6, width=0.5)
cbar.set_ticks(np.arange(min_val[n], max_val[n] + 1, step_val[n]))
cbar.ax.set_xticklabels(["0" if v == 0 else f"−{-v}" if v < 0 else f"{v}"
for v in range(min_val[n], max_val[n] + 1,
step_val[n])],
fontsize="large")
if n == 0:
cbar.set_label("Mean surface air temperature, °C",
fontsize="large")
else:
cbar.set_label("Mean grid cell contributions\n"
"to the ionospheric potential, V",
fontsize="large")
for axis in ["geo"]:
axM.spines[axis].set_linewidth(0.5)
axM.tick_params(top=True, right=True, which="both")
axM.tick_params(length=6, width=0.5)
axM.tick_params(length=3, width=0.5, which="minor")
axM.grid(False)
axM.set_xticks(np.arange(-180, 181, 30))
axM.set_xticks(np.arange(-180, 181, 10), minor=True)
axM.set_xticklabels(["180°", "", "120° W", "", "60° W", "", "0°",
"", "60° E", "", "120° E", "", "180°"],
fontsize="large", va="top")
# thin spaces (" ") are used between ‘°’ signs and letters
axM.set_yticks(np.arange(-90, 91, 30))
axM.set_yticks(np.arange(-90, 91, 10), minor=True)
axM.set_yticklabels(["", "", "30° S", "0°", "30° N", "", ""],
fontsize="large")
# thin spaces (" ") are used between ‘°’ signs and letters
axM.set_title("Partition of the Earth’s surface into regions",
fontsize="large", pad=10)
for j in range(6):
ax[2*j].bar(range(12), wrf_ip_sum_over_bands[j],
width=0.8, color=wrf_bands[j].color)
ax[2*j+1].bar(range(12), inm_ip_sum_over_bands[j],
width=0.8, color=inm_bands[j].color)
axT[j].bar(range(12), wrf_t2_avg_over_bands[j],
width=0.8, color=wrf_t2_bands[j].color)
offset = lambda p: transforms.ScaledTranslation(0, p/72,
fig.dpi_scale_trans)
for AX in [ax[2*j], ax[2*j+1]]:
AX.text(0.01 if j < 2 else 0.99, 1,
wrf_bands[j].title,
fontsize="large",
ha="left" if j < 2 else "right", va="top",
transform=AX.transAxes + offset(-2))
axT[j].text(0.01 if j < 2 else 0.99, 1,
wrf_t2_bands[j].title,
fontsize="large",
ha="left" if j < 2 else "right", va="top",
transform=axT[j].transAxes + offset(-2))
axD[0].imshow(T2_LxM,
cmap=cbar_cmap[0], norm=cbar_norm[0],
extent=[0, 12, 30, -30], aspect="auto", interpolation="none",
rasterized=True)
axD[1].imshow(WRF_IP_LxM,
cmap=cbar_cmap[1], norm=cbar_norm[1],
extent=[0, 12, 30, -30], aspect="auto", interpolation="none",
rasterized=True)
axD[2].imshow(INM_IP_LxM,
cmap=cbar_cmap[2], norm=cbar_norm[2],
extent=[0, 12, 30, -30], aspect="auto", interpolation="none",
rasterized=True)
axM.imshow(np.flip(strip_map[60:-60], axis=0), extent=[-180, 180, -30, 30])
for n in range(3):
axD[n].text(-0.28 * width/width, 1 + 0.1 * heightD/heightD,
chr(ord("a") + 2 * n), fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=axD[n].transAxes)
axT[0].text(-0.28 * width/width, 1 + 0.1 * heightD/heightT[0],
"b", fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=axT[0].transAxes)
for n in range(2):
ax[n].text(-0.28 * width/width, 1 + 0.1 * heightD/height[0],
chr(ord("d") + 2 * n), fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=ax[n].transAxes)
axM.text(-0.28 * width/widthM, 1 + 0.1 * heightD/heightM,
"g", fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=axM.transAxes)
fig.savefig("figures/ip_separate.eps", bbox_inches="tight")
Figure 2.6¶
New IP parameterisation
In [22]:
vostok_diurnal = np.load(
"./data/Vostok/vostok_diurnal_2006_2020_adjusted.npy"
)
vostok_seasonal = np.load(
"./data/Vostok/vostok_2006_2020_results_adjusted.npz"
)["mean"][0]
vostok_seasonal_counter = np.load(
"./data/Vostok/vostok_2006_2020_results_adjusted.npz"
)["counter"][0]
In [23]:
wrf_diurnal = wrf_hourly_total_ip[1000].mean(axis=0) / 240e3
wrf_new_diurnal = wrf_hourly_total_ip[500].mean(axis=0) / 240e3
vostok_diurnal /= vostok_diurnal.mean()
In [24]:
# making cycles
vostok_diurnal = [
(vostok_diurnal[0] + vostok_diurnal[-1]) / 2,
*vostok_diurnal,
(vostok_diurnal[0] + vostok_diurnal[-1]) / 2]
wrf_new_diurnal = [
*wrf_new_diurnal,
wrf_new_diurnal[0]
]
wrf_diurnal = [
*wrf_diurnal,
wrf_diurnal[0]
]
In [25]:
# selecting the latitude ranges for the analysis
band_names = ["0-9N", "0-9S", "9S-18S", "18S-90S", "9N-18N", "18N-90N"]
wrf_bands = [LatitudeBand(t, model="WRF") for t in band_names]
inmcm_bands = [LatitudeBand(t, model="INMCM") for t in band_names]
bar_data = np.zeros((2, len(band_names), 12))
# (subplots, strips, months)
for i, band in enumerate(wrf_bands):
bar_data[0, i] = wrf_LATxMON_ip[1000][band.slice].sum(axis=0)
bar_data[1, i] = wrf_LATxMON_ip[500][band.slice].sum(axis=0)
In [26]:
fig = plt.figure(figsize=(10, 15), constrained_layout=False)
ax = [None for _ in range(6)]
for n in range(6):
ax[n] = fig.add_subplot(4, 4, (2*n + 1, 2*n + 2))
low = [0.8, 0] + [0.8, 0e3] * 2
high = [1.2, 180] + [1.2, 300e3] * 2
step = [0.1, 30] + [0.1, 50e3] * 2
coeff = [1] * 2 + [1, 1e3] * 2
caption = ["Vostok station, 2006–2020 (adjusted)",
"Vostok station, 2006–2020 (adjusted)",
"WRF, 1980–2020, $\\varepsilon_0 = 1$ kJ/kg",
"WRF, 1980–2020, $\\varepsilon_0 = 1$ kJ/kg",
"WRF, 1980–2020, $\\varepsilon_0 = 0.5$ kJ/kg, $T_0 = 25$ °C",
"WRF, 1980–2020, $\\varepsilon_0 = 0.5$ kJ/kg, $T_0 = 25$ °C"]
for n in range(0, 5, 2):
for axis in ["top", "bottom", "left", "right"]:
ax[n].spines[axis].set_linewidth(0.5)
ax[n].tick_params(length=6, width=0.5)
ax[n].grid(color="0.", linewidth=0.5)
ax[n].set_xlim((0, 24))
ax[n].set_xticks(np.arange(0, 25, 2))
ax[n].set_xticklabels(np.arange(0, 25, 2), fontsize="large")
ax[n].set_xlabel("UTC, hours", fontsize="large")
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([f"{x:.1f}" for x in
np.arange(low[n], high[n] + step[n] / 2,
step[n]) / coeff[n]],
fontsize="large")
if n == 0:
ax[n].set_ylabel("Fair-weather pot. gradient\n"
"as a fraction of the mean",
fontsize="large")
else:
ax[n].set_ylabel("Ionospheric potential\nas a fraction of the mean",
fontsize="large")
ax[n].set_title(caption[n], fontsize="large")
for n in range(1, 6, 2):
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 == 1:
ax[n].set_ylabel("Monthly mean fair-weather\npotential gradient, V/m",
fontsize="large")
else:
ax[n].set_ylabel("Monthly mean\nionospheric potential, kV",
fontsize="large")
ax[n].set_title(caption[n], fontsize="large")
fig.align_ylabels([ax[0], ax[2], ax[4]])
fig.align_ylabels([ax[1], ax[3], ax[5]])
# diurnal variation curves
ax[0].plot(np.append(np.insert(np.arange(0.5, 24, 1), 0, 0), 24),
vostok_diurnal,
linewidth=1.5, color="purple", clip_on=False,
zorder=4)
ax[2].plot(np.arange(0, 25),
wrf_diurnal,
linewidth=1.5, color="red", clip_on=False,
zorder=4)
ax[4].plot(np.arange(0, 25),
wrf_new_diurnal,
linewidth=1.5, color="red", clip_on=False,
zorder=4)
# arrows and amplitudes
for ax_number, data_diurnal in zip(
[0, 2, 4], [vostok_diurnal, wrf_diurnal, wrf_new_diurnal]
):
ampl = np.max(data_diurnal) - np.min(data_diurnal)
draw_arrows(ax=ax[ax_number],
miny=np.min(data_diurnal),
maxy=np.max(data_diurnal),
arrow_pos=25,
text_pos=25.2,
ampl=ampl)
# seasonal variation of the PG
ax[1].bar(np.arange(12), vostok_seasonal, width=0.8, color="orangered")
# arrows for the seasonal variation of the PG
vostok_ampl = (vostok_seasonal.max() - vostok_seasonal.min()) / \
np.average(vostok_seasonal, weights=vostok_seasonal_counter)
draw_arrows(ax=ax[1],
miny=np.min(vostok_seasonal),
maxy=np.max(vostok_seasonal),
exty=high[1] / 10,
ampl=vostok_ampl,
inward=True)
# seasonal variation of the IP (decomposed into contributions)
for n in range(2):
bottom_data = np.zeros(12)
for k in range(len(band_names)):
data = bar_data[n, k]
ax[2*n+3].bar(np.arange(12), data,
width=0.8, color=wrf_bands[k].color,
bottom=bottom_data, label=wrf_bands[k].title)
bottom_data += data
# arrows for the seasonal variations of the IP
variation_1000 = wrf_LATxMON_ip[1000].sum(axis=0)
variation_500 = wrf_LATxMON_ip[500].sum(axis=0)
for ax_number, seasonal_variation in zip([3, 5],
[variation_1000, variation_500]):
seasonal_ampl = (seasonal_variation.max() - seasonal_variation.min()) / \
np.average(seasonal_variation, weights=wrf_numdays)
draw_arrows(ax=ax[ax_number],
miny=seasonal_variation.min(),
maxy=seasonal_variation.max(),
exty=high[ax_number] / 10,
ampl=seasonal_ampl,
inward=True)
for n in range(0, 5, 2):
ax[n].text(-0.28, 1.05, chr(ord("a") + n),
fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=ax[n].transAxes)
for n in range(1, 6, 2):
ax[n].text(-0.30, 1.05, chr(ord("a") + n),
fontsize="x-large",
fontweight="semibold", ha="left", va="bottom",
transform=ax[n].transAxes)
fig.subplots_adjust(hspace=0.45, wspace=1.6)
leg_width_rel = (ax[5].get_position().x1 - ax[4].get_position().x0) \
/ (ax[5].get_position().x1 - ax[5].get_position().x0)
leg = ax[5].legend(bbox_to_anchor=(1. - leg_width_rel, -1./2,
leg_width_rel, 1./6), ncols=6,
loc="lower center", borderaxespad=0.,
mode="expand", fontsize="large",
framealpha=1, edgecolor="0.",
handlelength=1., handletextpad=0.5)
leg.get_frame().set_linewidth(0.5)
fig.savefig("figures/new_parameterisation.eps", bbox_inches="tight")
In [ ]: