LPT_analyze_stats.py

# %%
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import param

plt.rcParams["font.family"] = "Arial"
plt.rcParams["font.serif"] = "Computer Modern Roman"
plt.rcParams["font.size"] = 12

# %%
number_particles = 400000
is_2D_or_3D = "3D"
nu = "10-6"
simu_name = "one_month_release_constant_seed_by_site"
mesh_name = "gbr_styx_with_rivers"

# %%
p = param.parameters(mesh_name)
base_dir = f"{p.output_directory}/LPT_sediments/{is_2D_or_3D}/{simu_name}_nu{nu}_{str(number_particles)}/"
out_dir = f"{base_dir}/stats/"
fig_dir = f"{out_dir}/fig/"
os.makedirs(fig_dir, exist_ok=True)

sed_sizes = [1, 2, 4, 8, 16, 32, 64, 125, 250, 500, 1000]
n_sites = 8

inside = ["site 1", "site 2", "site 3", "site 4"]
outside = ["site 5", "site 6", "site 7", "site 8"]
# %%
full_summary_df = pd.DataFrame()
reach_summary_df = pd.DataFrame()
sed_summary_df = pd.DataFrame()

inside_full_summary_df = pd.DataFrame()
outside_full_summary_df = pd.DataFrame()

inside_reach_df = pd.DataFrame()
inside_sed_df = pd.DataFrame()
inside_total_df = pd.DataFrame()

outside_reach_df = pd.DataFrame()
outside_sed_df = pd.DataFrame()
outside_total_df = pd.DataFrame()

time_inside_df = pd.DataFrame()
time_outside_df = pd.DataFrame()
time_total_df = pd.DataFrame()

# for sed in sediments:
for min_size, max_size in zip(sed_sizes[:-1], sed_sizes[1:]):
    sed = f"{min_size}_{max_size}"

    stats_dir = base_dir + sed + "/stats/"

    reach_df = pd.DataFrame()
    sed_df = pd.DataFrame()
    total_df = pd.DataFrame()

    site_df = pd.DataFrame()

    time_sed_df = pd.read_csv(f"{stats_dir}stats_{sed}_sedimentation.csv", index_col=0)

    time_inside_df = time_inside_df.append(time_sed_df.loc["Inside"].rename(sed))
    time_outside_df = time_outside_df.append(time_sed_df.loc["Outside"].rename(sed))
    time_total_df = time_total_df.append(time_sed_df.mean().rename(sed))

    for i in range(n_sites):
        tmp_df = pd.read_csv(
            f"{stats_dir}stats_{sed}_release_point_{i}.csv",
            index_col=0,
            parse_dates=True,
        )
        reach_df[f"site {i+1}"] = tmp_df["n_ever_reached_zoi"]
        sed_df[f"site {i+1}"] = tmp_df["n_ever_sedimented_on_zoi"]
        total_df[f"site {i+1}"] = tmp_df["n_from_site"]

        site_df[f"site {i+1}"] = tmp_df.iloc[-1]

    site_df["inside"] = site_df[inside].sum(axis=1)
    site_df["outside"] = site_df[outside].sum(axis=1)

    reach_df["total"] = reach_df.sum(axis=1)
    reach_df["inside"] = reach_df[inside].sum(axis=1)
    reach_df["outside"] = reach_df[outside].sum(axis=1)

    sed_df["total"] = sed_df.sum(axis=1)
    sed_df["inside"] = sed_df[inside].sum(axis=1)
    sed_df["outside"] = sed_df[outside].sum(axis=1)

    total_df["total"] = total_df.sum(axis=1)
    total_df["inside"] = total_df[inside].sum(axis=1)
    total_df["outside"] = total_df[outside].sum(axis=1)

    inside_reach_df[sed] = reach_df["inside"]
    inside_sed_df[sed] = sed_df["inside"]
    inside_total_df[sed] = total_df["inside"]

    outside_reach_df[sed] = reach_df["outside"]
    outside_sed_df[sed] = sed_df["outside"]
    outside_total_df[sed] = total_df["outside"]

    full_summary_df[sed] = site_df.sum(axis=1)
    inside_full_summary_df[sed] = site_df[inside].sum(axis=1)
    outside_full_summary_df[sed] = site_df[outside].sum(axis=1)

    reach_summary_df[sed] = (
        site_df.loc["n_ever_reached_zoi"] / site_df.loc["n_from_site"] * 100
    )
    sed_summary_df[sed] = (
        site_df.loc["n_ever_sedimented_on_zoi"] / site_df.loc["n_from_site"] * 100
    )

    reach_df.to_csv(f"{stats_dir}stats_{sed}_reach_summary.csv")
    sed_df.to_csv(f"{stats_dir}stats_{sed}_sed_summary.csv")
    total_df.to_csv(f"{stats_dir}stats_{sed}_total_summary.csv")


full_summary_df = full_summary_df.T
full_summary_df["prop_ever_reached"] = (
    full_summary_df["n_ever_reached_zoi"] / full_summary_df["n_from_site"] * 100
)
full_summary_df["prop_ever_sedimented"] = (
    full_summary_df["n_ever_sedimented_on_zoi"] / full_summary_df["n_from_site"] * 100
)
full_summary_df.to_csv(f"{out_dir}full_summary.csv")

inside_full_summary_df = inside_full_summary_df.T
inside_full_summary_df["prop_ever_reached"] = (
    inside_full_summary_df["n_ever_reached_zoi"]
    / inside_full_summary_df["n_from_site"]
    * 100
)
inside_full_summary_df["prop_ever_sedimented"] = (
    inside_full_summary_df["n_ever_sedimented_on_zoi"]
    / inside_full_summary_df["n_from_site"]
    * 100
)
inside_full_summary_df.to_csv(f"{out_dir}inside_full_summary.csv")

outside_full_summary_df = outside_full_summary_df.T
outside_full_summary_df["prop_ever_reached"] = (
    outside_full_summary_df["n_ever_reached_zoi"]
    / outside_full_summary_df["n_from_site"]
    * 100
)
outside_full_summary_df["prop_ever_sedimented"] = (
    outside_full_summary_df["n_ever_sedimented_on_zoi"]
    / outside_full_summary_df["n_from_site"]
    * 100
)
outside_full_summary_df.to_csv(f"{out_dir}outside_full_summary.csv")

reach_summary_df.loc["mean"] = reach_summary_df.mean()
sed_summary_df.loc["mean"] = sed_summary_df.mean()

reach_summary_df.to_csv(f"{out_dir}reach_summary.csv")
sed_summary_df.to_csv(f"{out_dir}sed_summary.csv")

prop_inside_reach_df = inside_reach_df / inside_total_df * 100
prop_inside_sed_df = inside_sed_df / inside_total_df * 100
prop_outside_reach_df = outside_reach_df / outside_total_df * 100
prop_outside_sed_df = outside_sed_df / outside_total_df * 100

time_inside_df.to_csv((f"{out_dir}time_inside_summary.csv"))
time_outside_df.to_csv((f"{out_dir}time_outside_summary.csv"))
time_total_df.to_csv((f"{out_dir}time_total_summary.csv"))

# %% Figures param
resample_period = "12H"
sed_sel_for_plot = sed_sizes[:7]
n_lines = len(sed_sel_for_plot) - 1
colors = plt.cm.gist_heat(np.linspace(0, 0.85, n_lines))

# %%
for stat in ["reach", "sedimented"]:

    if stat == "reach":
        df_inside = prop_inside_reach_df.resample(resample_period).mean()
        df_outside = prop_outside_reach_df.resample(resample_period).mean()
    else:
        df_inside = prop_inside_sed_df.resample(resample_period).mean()
        df_outside = prop_outside_sed_df.resample(resample_period).mean()

    fig, ax = plt.subplots(nrows=2, sharex=True, sharey=True)

    for i, (min_size, max_size) in enumerate(
        zip(sed_sel_for_plot[:-1], sed_sel_for_plot[1:])
    ):
        sed = f"{min_size}_{max_size}"
        df_inside[sed].plot(c=colors[i], ax=ax[0], label=f"{min_size}-{max_size} µm")
        df_outside[sed].plot(c=colors[i], ax=ax[1], label="")

    for a in ax:
        a.axvspan("2021-01-05 00:00:00", "2021-02-05 00:00:00", facecolor="0.85")
        a.set_xlim("2021-01-01 00:00:00", "2021-04-01 00:00:00")
        a.set_ylim(0, 100)
        a.grid(ls=":")

    ax[0].set_ylabel("Inside release")
    ax[1].set_ylabel("Outside release")

    ax[1].tick_params(axis="x", which="major", pad=-8)

    plt.subplots_adjust(hspace=0.1)

    fig.legend(bbox_to_anchor=(0, 0.85, 1, 0.2), loc="upper center", ncol=3)
    fig.text(
        -0.01, 0.5, "Proportion of sediments [%]", va="center", rotation="vertical"
    )

    plt.savefig(f"{fig_dir}/prop_sed_{stat}.pdf", dpi=1000, bbox_inches="tight")
    plt.savefig(f"{fig_dir}/prop_sed_{stat}.jpg", dpi=1000, bbox_inches="tight")
    plt.savefig(
        f"{fig_dir}/prop_sed_{stat}.png",
        dpi=1000,
        bbox_inches="tight",
        transparent=True,
    )
    plt.show()
    plt.close()

# %%
fig = plt.figure(constrained_layout=True, figsize=(8, 4))
subfigs = fig.subfigures(1, 2)

# Left side
axsLeft = subfigs[0].subplots(nrows=2, sharex=True, sharey=True)

df_inside = prop_inside_reach_df.resample(resample_period).mean()
df_outside = prop_outside_reach_df.resample(resample_period).mean()

for i, (min_size, max_size) in enumerate(
    zip(sed_sel_for_plot[:-1], sed_sel_for_plot[1:])
):
    sed = f"{min_size}_{max_size}"
    df_inside[sed].plot(c=colors[i], ax=axsLeft[0], label=f"{min_size}-{max_size} µm")
    df_outside[sed].plot(c=colors[i], ax=axsLeft[1], label="")

for a in axsLeft:
    a.axvspan("2021-01-05 00:00:00", "2021-02-05 00:00:00", facecolor="0.85")
    a.set_xlim("2021-01-01 00:00:00", "2021-04-01 00:00:00")
    a.set_ylim(0, 100)
    a.grid(ls=":")

axsLeft[0].set_ylabel("Inside release")
axsLeft[1].set_ylabel("Outside release")

axsLeft[1].tick_params(axis="x", which="minor", label1On=False)
axsLeft[1].tick_params(axis="x", which="major", pad=-20)

subfigs[0].suptitle(
    "Percentage of particles that ever\n" + r"$\bf{reached}$ the dugong sanctuary",
    fontsize="medium",
)

# Right side
axsRight = subfigs[1].subplots(nrows=2, sharex=True, sharey=True)

df_inside = prop_inside_sed_df.resample(resample_period).mean()
df_outside = prop_outside_sed_df.resample(resample_period).mean()

for i, (min_size, max_size) in enumerate(
    zip(sed_sel_for_plot[:-1], sed_sel_for_plot[1:])
):
    sed = f"{min_size}_{max_size}"
    df_inside[sed].plot(c=colors[i], ax=axsRight[0], label="")
    df_outside[sed].plot(c=colors[i], ax=axsRight[1], label="")

for a in axsRight:
    a.axvspan("2021-01-05 00:00:00", "2021-02-05 00:00:00", facecolor="0.85")
    a.set_xlim("2021-01-01 00:00:00", "2021-04-01 00:00:00")
    a.set_ylim(0, 100)
    a.grid(ls=":")

axsRight[1].tick_params(axis="x", which="minor", label1On=False)
axsRight[1].tick_params(axis="x", which="major", pad=-20)

subfigs[1].suptitle(
    "Percentage of particles that ever\n" + r"$\bf{settled}$ on the dugong sanctuary",
    fontsize="medium",
)

fig.legend(
    bbox_to_anchor=(0, -0.1, 1, 0),
    loc="lower center",
    ncol=n_lines,
    fancybox=False,
    columnspacing=1.6,
)
fig.suptitle(" ", fontsize="xx-small")
plt.savefig(f"{fig_dir}/prop_sed_both.pdf", dpi=1000, bbox_inches="tight")
plt.savefig(f"{fig_dir}/prop_sed_both.jpg", dpi=1000, bbox_inches="tight")
plt.savefig(
    f"{fig_dir}/prop_sed_both.png", dpi=1000, bbox_inches="tight", transparent=True
)
plt.show()
plt.close()

# %%
names = [f"{mi}-{ma} µm" for mi, ma in zip(sed_sizes[:-1], sed_sizes[1:])]
names

col1 = "#005AB5"
col2 = "#DC3220"

fig, ax1 = plt.subplots(figsize=(5, 4))

ax2 = ax1.twinx()

ax1.plot(
    time_total_df.index,
    time_total_df.nb_time_sed - 1,
    c=col1,
    linestyle=":",
    marker="o",
)
ax2.plot(
    time_total_df.index,
    time_total_df.percentage_time_sed,
    c=col2,
    linestyle=":",
    marker="o",
)

ax1.set_xticks(np.arange(len(time_total_df)))
ax1.set_xticklabels(labels=names, rotation=45, ha="right")
ax1.grid(which="major", axis="x", linestyle=":")

ax1.set_ylabel("Mean number of resuspensions", color=col1)
ax1.yaxis.label.set_color(col1)
ax1.tick_params(axis="y", colors=col1)

ax2.set_ylabel("Mean percentage of time settled", color=col2, rotation=-90, labelpad=10)
ax2.yaxis.label.set_color(col2)
ax2.tick_params(axis="y", colors=col2)
ax2.spines["left"].set_color(col1)
ax2.spines["right"].set_color(col2)

# To add horizontal lines and align scales
ax1.set_ylim(-50, 1050)
ax2.grid(which="major", axis="y", linestyle=":")
ax2.set_ylim(-5, 105)

plt.savefig(f"{fig_dir}/time_sed.pdf", dpi=1000, bbox_inches="tight")
plt.savefig(f"{fig_dir}/time_sed.jpg", dpi=1000, bbox_inches="tight")
plt.savefig(f"{fig_dir}/time_sed.png", dpi=1000, bbox_inches="tight", transparent=True)
plt.show()

# %%