Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[plotter] Feat: Implement $F_d$, $F_k$ and other improvements #44

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

[plotter] Feat: Implement $F_d$, $F_k$ and other improvements #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0385019
plotter: Plot prettification
mixayloff-dimaaylov Mar 3, 2023
9ccb0c7
plotter: feat: add F_d, d, F_k calculation
mixayloff-dimaaylov Mar 3, 2023
62b04bb
perf: Replace conda with mamba
mixayloff-dimaaylov Mar 3, 2023
5b863e6
plotter: feat: l_s calculation
mixayloff-dimaaylov Dec 19, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/jupyter/ClickHouseExample.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
"# Install a conda packages in the current Jupyter kernel\n",
"import sys\n",
"\n",
"!conda install --yes --prefix {sys.prefix} -c conda-forge clickhouse-driver clickhouse-sqlalchemy ipython-sql"
"!mamba install --yes --prefix {sys.prefix} -c conda-forge clickhouse-driver clickhouse-sqlalchemy ipython-sql"
]
},
{
Expand Down
136 changes: 111 additions & 25 deletions docs/jupyter/plotter.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@
"# Install a conda packages in the current Jupyter kernel\n",
"import sys\n",
"\n",
"!conda install --yes --prefix {sys.prefix} -c conda-forge clickhouse-driver clickhouse-sqlalchemy ipywidgets"
"!mamba install --yes --prefix {sys.prefix} -c conda-forge clickhouse-driver clickhouse-sqlalchemy ipywidgets"
]
},
{
Expand Down Expand Up @@ -167,7 +167,50 @@
"\n",
"\n",
"def sigPhi(sigNT, f):\n",
" return 1e16 * 80.8 * pi * sigNT / (C * f)"
" return 1e16 * 80.8 * pi * sigNT / (C * f)\n",
"\n",
"\n",
"def F_d(avgNT, f_0, alpha):\n",
" return np.sqrt(np.absolute((C * f_0 ** 3) /\n",
" (80.8 * pi * (avgNT * 1e16) / np.sin(np.radians(alpha)))))\n",
"\n",
"\n",
"def d(h_max, l_s, f_0, alpha):\n",
" return (4 * (h_max / np.sin(np.radians(alpha))) ** 2 * (C ** 2) /\n",
" ((pi ** 2) * (f_0 ** 2) * (l_s ** 4)))\n",
"\n",
"\n",
"def F_k(sigPhi, h_max, l_s, f_0, alpha):\n",
" d_ = d(h_max, l_s, f_0, alpha)\n",
" return (np.sqrt(2) * f_0 /\n",
" (sigPhi * np.sqrt(1.0 + d_ / 2)))\n",
"\n",
"\n",
"def autocorr(x, t):\n",
" if t <= 0:\n",
" return 1\n",
"\n",
" return np.corrcoef(x[:-t], x[t:])[0, 1]\n",
"\n",
"\n",
"def R_delNT(delNT, sigNT):\n",
" n = 50\n",
"\n",
" res = np.zeros((n-1, n))\n",
" for delnt, signt in zip(\n",
" np.lib.stride_tricks.sliding_window_view(delNT, n),\n",
" np.lib.stride_tricks.sliding_window_view(sigNT, n)):\n",
" acc = np.empty((1, n))\n",
" for i in range(0, n):\n",
" acc[0, i] = autocorr(delnt, i)/signt[i]\n",
" res = np.append(res, acc, axis=0)\n",
"\n",
" return res\n",
"\n",
"\n",
"def l_s(delNT, sigNT):\n",
" k = (R_delNT(delNT, sigNT) < np.exp(-1)).argmax(axis=1)\n",
" return 1000 * (k + 1) * 0.02"
]
},
{
Expand Down Expand Up @@ -375,10 +418,18 @@
" df_range['sigNT'] = pd.Series(sigNT(df_range.delNT)).shift(59, fill_value=0.0)\n",
" df_range['sigPhi'] = sigPhi(df_range.sigNT, df_range.f2)\n",
"\n",
" df_range['l_s'] = l_s(df_range.delNT, df_range.sigNT)\n",
"\n",
" h_max = 300000\n",
" df_range['F_d'] = F_d(df_range.avgNT, df_range.f1, df_satxyz2.elevation)\n",
" df_range['d'] = d(h_max, df_range.l_s, df_range.f1, df_satxyz2.elevation)\n",
" df_range['F_k'] = F_k(df_range.sigPhi, h_max, df_range.l_s, df_range.f1, df_satxyz2.elevation)\n",
"\n",
" # For export\n",
" df_range['ism_tec'] = df_ismrawtec.tec\n",
" df_range['ism_primaryfreq'] = df_ismrawtec.primaryfreq\n",
" df_range['ism_secondaryfreq'] = df_ismrawtec.secondaryfreq\n",
" df_range['elevation'] = df_satxyz2.elevation\n",
"\n",
" return df_range\n",
"\n",
Expand All @@ -403,44 +454,79 @@
" locator = mdates.AutoDateLocator()\n",
" formatter = mdates.ConciseDateFormatter(locator)\n",
"\n",
" gfig, gax = plt.subplots()\n",
" gax.xaxis.set_major_locator(locator)\n",
" gax.xaxis.set_major_formatter(formatter)\n",
"\n",
" def dumpplot(xs, ys, vname):\n",
" def dumpplot(gax, xs, ys, yname, ylabel):\n",
" fig, ax = plt.subplots()\n",
"\n",
" ax.xaxis.set_major_locator(locator)\n",
" ax.xaxis.set_major_formatter(formatter)\n",
"\n",
" ax.set_title(f\"{vname} {track_name_human}\")\n",
" ax.set_title(f\"${yname}$ {track_name_human}\")\n",
" ax.set_xlabel(\"Datetime\")\n",
" ax.plot(xs, ys, label=vname)\n",
" ax.set_ylabel(f\"${ylabel}$\")\n",
" # Cuttoff filter splashes\n",
" ax.plot(xs[200:], ys[200:], label=f\"${yname}$\")\n",
" ax.grid()\n",
" # Rotate and align the tick labels so they look better.\n",
" fig.autofmt_xdate()\n",
" fig.legend()\n",
"\n",
" plt.title(f\"{vname} {track_name_human}\")\n",
" plt.savefig(f\"{_path}/{track_name} {vname}.png\")\n",
" plt.title(f\"${yname}$ {track_name_human}\")\n",
" plt.savefig(f\"{_path}/{track_name} {yname}.png\")\n",
" plt.close(fig)\n",
"\n",
" gax.plot(xs, ys, label=vname)\n",
" gax.plot(xs[200:], ys[200:], label=f\"${yname}$\")\n",
" gax.set_xlabel(\"Datetime\")\n",
"\n",
" dumpplot(sat.time, sat.NTpsr, \"NT(P1-P2)\")\n",
" dumpplot(sat.time, sat.NTadr, \"NT(adr1-adr2)\")\n",
" dumpplot(sat.time, sat.ism_tec, \"ISMRAWTEC's TEC\")\n",
" dumpplot(sat.time, sat.avgNT, \"avgNT\")\n",
" dumpplot(sat.time, sat.delNT, \"delNT\")\n",
" dumpplot(sat.time, sat.sigNT, \"sigNT\")\n",
" dumpplot(sat.time, sat.sigPhi, \"sigPhi\")\n",
"\n",
" gax.legend()\n",
" gax.grid()\n",
" plt.title(f\"ПЭСы {track_name_human} ({_date})\")\n",
" # Rotate and align the tick labels so they look better.\n",
" gfig.autofmt_xdate()\n",
" def init_plot():\n",
" fig, ax = plt.subplots()\n",
" ax.xaxis.set_major_locator(locator)\n",
" ax.xaxis.set_major_formatter(formatter)\n",
" return fig, ax\n",
"\n",
" def plot_finalize(gax, title):\n",
" gax.legend()\n",
" gax.grid()\n",
" plt.title(f\"{title} {track_name_human} ({_date})\")\n",
" # Rotate and align the tick labels so they look better.\n",
" gfig.autofmt_xdate()\n",
"\n",
" # First plot with TECU\n",
" gfig, gax = init_plot()\n",
" dumpplot(gax, sat.time, sat.NTpsr, \"N_T (P_1 - P_2)\", \"TECU\")\n",
" dumpplot(gax, sat.time, sat.NTadr, \"N_T (adr_1 - adr_2)\", \"TECU\")\n",
" dumpplot(gax, sat.time, sat.ism_tec, \"ISMRAWTEC's TEC\", \"TECU\")\n",
" dumpplot(gax, sat.time, sat.avgNT, \"\\overline{{N_T}}\", \"TECU\")\n",
" dumpplot(gax, sat.time, sat.delNT, \"\\Delta N_T\", \"TECU\")\n",
" dumpplot(gax, sat.time, sat.sigNT, \"\\sigma N_T\", \"TECU\")\n",
" dumpplot(gax, sat.time, sat.sigPhi, \"\\sigma \\\\varphi\", \"TECU\")\n",
" gax.set_ylabel(\"TECU\")\n",
" plot_finalize(gax, \"ПЭСы\")\n",
"\n",
" # Other plots:\n",
" gfig, gax = init_plot()\n",
" dumpplot(gax, sat.time, sat.l_s, \"l_s\", \"m\")\n",
" gax.set_ylabel(\"безразмерная\")\n",
" plot_finalize(gax, \"$l_s$\")\n",
"\n",
" gfig, gax = init_plot()\n",
" dumpplot(gax, sat.time, sat.F_d, \"F_d\", \"Hz\")\n",
" gax.set_ylabel(\"Hz\")\n",
" plot_finalize(gax, \"$F_d$\")\n",
"\n",
" gfig, gax = init_plot()\n",
" dumpplot(gax, sat.time, sat.F_k, \"F_k\", \"Hz\")\n",
" gax.set_ylabel(\"Hz\")\n",
" plot_finalize(gax, \"$F_k$\")\n",
"\n",
" gfig, gax = init_plot()\n",
" dumpplot(gax, sat.time, sat.d, \"d\", \"безразмерная\")\n",
" gax.set_ylabel(\"безразмерная\")\n",
" plot_finalize(gax, \"$d$\")\n",
"\n",
" gfig, gax = init_plot()\n",
" dumpplot(gax, sat.time, sat.elevation, \"elevation\", \"Deg\")\n",
" gax.set_ylabel(\"Deg\")\n",
" plot_finalize(gax, \"Угол возвышения\")\n",
"\n",
"\n",
"# Ref: https://stackoverflow.com/questions/15411967\n",
Expand Down
134 changes: 110 additions & 24 deletions util/plotter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
# Install a conda packages in the current Jupyter kernel
import sys

get_ipython().system('conda install --yes --prefix {sys.prefix} -c conda-forge clickhouse-driver clickhouse-sqlalchemy ipywidgets')
get_ipython().system('mamba install --yes --prefix {sys.prefix} -c conda-forge clickhouse-driver clickhouse-sqlalchemy ipywidgets')


# ## Исходная программа
Expand Down Expand Up @@ -131,6 +131,49 @@ def sigPhi(sigNT, f):
return 1e16 * 80.8 * pi * sigNT / (C * f)


def F_d(avgNT, f_0, alpha):
return np.sqrt(np.absolute((C * f_0 ** 3) /
(80.8 * pi * (avgNT * 1e16) / np.sin(np.radians(alpha)))))


def d(h_max, l_s, f_0, alpha):
return (4 * (h_max / np.sin(np.radians(alpha))) ** 2 * (C ** 2) /
((pi ** 2) * (f_0 ** 2) * (l_s ** 4)))


def F_k(sigPhi, h_max, l_s, f_0, alpha):
d_ = d(h_max, l_s, f_0, alpha)
return (np.sqrt(2) * f_0 /
(sigPhi * np.sqrt(1.0 + d_ / 2)))


def autocorr(x, t):
if t <= 0:
return 1

return np.corrcoef(x[:-t], x[t:])[0, 1]


def R_delNT(delNT, sigNT):
n = 50

res = np.zeros((n-1, n))
for delnt, signt in zip(
np.lib.stride_tricks.sliding_window_view(delNT, n),
np.lib.stride_tricks.sliding_window_view(sigNT, n)):
acc = np.empty((1, n))
for i in range(0, n):
acc[0, i] = autocorr(delnt, i)/signt[i]
res = np.append(res, acc, axis=0)

return res


def l_s(delNT, sigNT):
k = (R_delNT(delNT, sigNT) < np.exp(-1)).argmax(axis=1)
return 1000 * (k + 1) * 0.02


# ### Работа с выгрузками

# In[ ]:
Expand Down Expand Up @@ -314,10 +357,18 @@ def perf_cal(values):
df_range['sigNT'] = pd.Series(sigNT(df_range.delNT)).shift(59, fill_value=0.0)
df_range['sigPhi'] = sigPhi(df_range.sigNT, df_range.f2)

df_range['l_s'] = l_s(df_range.delNT, df_range.sigNT)

h_max = 300000
df_range['F_d'] = F_d(df_range.avgNT, df_range.f1, df_satxyz2.elevation)
df_range['d'] = d(h_max, df_range.l_s, df_range.f1, df_satxyz2.elevation)
df_range['F_k'] = F_k(df_range.sigPhi, h_max, df_range.l_s, df_range.f1, df_satxyz2.elevation)

# For export
df_range['ism_tec'] = df_ismrawtec.tec
df_range['ism_primaryfreq'] = df_ismrawtec.primaryfreq
df_range['ism_secondaryfreq'] = df_ismrawtec.secondaryfreq
df_range['elevation'] = df_satxyz2.elevation

return df_range

Expand All @@ -342,44 +393,79 @@ def plot_build(sat):
locator = mdates.AutoDateLocator()
formatter = mdates.ConciseDateFormatter(locator)

gfig, gax = plt.subplots()
gax.xaxis.set_major_locator(locator)
gax.xaxis.set_major_formatter(formatter)

def dumpplot(xs, ys, vname):
def dumpplot(gax, xs, ys, yname, ylabel):
fig, ax = plt.subplots()

ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)

ax.set_title(f"{vname} {track_name_human}")
ax.set_title(f"${yname}$ {track_name_human}")
ax.set_xlabel("Datetime")
ax.plot(xs, ys, label=vname)
ax.set_ylabel(f"${ylabel}$")
# Cuttoff filter splashes
ax.plot(xs[200:], ys[200:], label=f"${yname}$")
ax.grid()
# Rotate and align the tick labels so they look better.
fig.autofmt_xdate()
fig.legend()

plt.title(f"{vname} {track_name_human}")
plt.savefig(f"{_path}/{track_name} {vname}.png")
plt.title(f"${yname}$ {track_name_human}")
plt.savefig(f"{_path}/{track_name} {yname}.png")
plt.close(fig)

gax.plot(xs, ys, label=vname)
gax.plot(xs[200:], ys[200:], label=f"${yname}$")
gax.set_xlabel("Datetime")

dumpplot(sat.time, sat.NTpsr, "NT(P1-P2)")
dumpplot(sat.time, sat.NTadr, "NT(adr1-adr2)")
dumpplot(sat.time, sat.ism_tec, "ISMRAWTEC's TEC")
dumpplot(sat.time, sat.avgNT, "avgNT")
dumpplot(sat.time, sat.delNT, "delNT")
dumpplot(sat.time, sat.sigNT, "sigNT")
dumpplot(sat.time, sat.sigPhi, "sigPhi")

gax.legend()
gax.grid()
plt.title(f"ПЭСы {track_name_human} ({_date})")
# Rotate and align the tick labels so they look better.
gfig.autofmt_xdate()
def init_plot():
fig, ax = plt.subplots()
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)
return fig, ax

def plot_finalize(gax, title):
gax.legend()
gax.grid()
plt.title(f"{title} {track_name_human} ({_date})")
# Rotate and align the tick labels so they look better.
gfig.autofmt_xdate()

# First plot with TECU
gfig, gax = init_plot()
dumpplot(gax, sat.time, sat.NTpsr, "N_T (P_1 - P_2)", "TECU")
dumpplot(gax, sat.time, sat.NTadr, "N_T (adr_1 - adr_2)", "TECU")
dumpplot(gax, sat.time, sat.ism_tec, "ISMRAWTEC's TEC", "TECU")
dumpplot(gax, sat.time, sat.avgNT, "\overline{{N_T}}", "TECU")
dumpplot(gax, sat.time, sat.delNT, "\Delta N_T", "TECU")
dumpplot(gax, sat.time, sat.sigNT, "\sigma N_T", "TECU")
dumpplot(gax, sat.time, sat.sigPhi, "\sigma \\varphi", "TECU")
gax.set_ylabel("TECU")
plot_finalize(gax, "ПЭСы")

# Other plots:
gfig, gax = init_plot()
dumpplot(gax, sat.time, sat.l_s, "l_s", "m")
gax.set_ylabel("безразмерная")
plot_finalize(gax, "$l_s$")

gfig, gax = init_plot()
dumpplot(gax, sat.time, sat.F_d, "F_d", "Hz")
gax.set_ylabel("Hz")
plot_finalize(gax, "$F_d$")

gfig, gax = init_plot()
dumpplot(gax, sat.time, sat.F_k, "F_k", "Hz")
gax.set_ylabel("Hz")
plot_finalize(gax, "$F_k$")

gfig, gax = init_plot()
dumpplot(gax, sat.time, sat.d, "d", "безразмерная")
gax.set_ylabel("безразмерная")
plot_finalize(gax, "$d$")

gfig, gax = init_plot()
dumpplot(gax, sat.time, sat.elevation, "elevation", "Deg")
gax.set_ylabel("Deg")
plot_finalize(gax, "Угол возвышения")


# Ref: https://stackoverflow.com/questions/15411967
Expand Down