analyse_shap_and_features.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import matplotlib\n",
    "import os\n",
    "\n",
    "import ccc\n",
    "import utils\n",
    "import utils_plots\n",
    "import utils_shap\n",
    "\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'size' : 12}\n",
    "matplotlib.rc('font', **font)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_root_tm_path = os.path.join(ccc.MODEL_ROOT_PATH, 'targetmode_1', '2022_02_21__11-11__ALDIS_paper')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_path = os.path.join(model_root_tm_path, \"_model_00018_test_df.pickle\")\n",
    "\n",
    "dd = pd.read_pickle(df_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PLOT_COLS = [f\"{c}_relative_shap\" for c in (ccc.LVL_TRAIN_COLS + [\"hour\", \"dayofyear\"])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap_path = os.path.join(model_root_tm_path, \"_model_00018_shap_parquet_bg_by_lon_lat_no_flash\")\n",
    "dshap = pd.read_parquet(shap_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dshap.rename(columns={utils_shap.colname_meta_infix(col) : col for col in ccc.INDEX_COLS}, inplace=True)\n",
    "dshap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(dshap[\"shap_base_value\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "used_threshold = 0.8708981871604919  # F1 score based threshold of epoch 18 on validation set\n",
    "\n",
    "vc_threshold = utils.getVeryConfidentThreshold(used_threshold)\n",
    "\n",
    "df_joined = utils.joinDataframes(dshap, dd)\n",
    "\n",
    "df_joined.loc[:, \"pred_class\"] = np.where(df_joined[\"output\"] > used_threshold, \"pred_flash\", \"pred_no_flash\")\n",
    "df_joined.loc[:, \"real_class\"] = np.where(df_joined[\"target\"] > 0.5, \"real_flash\", \"real_no_flash\")  # target col only contains 0s and 1s.\n",
    "\n",
    "df_joined.loc[:, 'cat'] = np.select(\n",
    "    [\n",
    "        (df_joined['pred_class'] == \"pred_flash\") & (df_joined['real_class'] == \"real_flash\") & (df_joined[\"output\"] < vc_threshold), \n",
    "        (df_joined['pred_class'] == \"pred_flash\") & (df_joined['real_class'] == \"real_flash\") & (df_joined[\"output\"] >= vc_threshold), \n",
    "        (df_joined['pred_class'] == \"pred_no_flash\") & (df_joined['real_class'] == \"real_flash\"), \n",
    "        (df_joined['pred_class'] == \"pred_flash\") & (df_joined['real_class'] == \"real_no_flash\"), \n",
    "        (df_joined['pred_class'] == \"pred_no_flash\") & (df_joined['real_class'] == \"real_no_flash\"), \n",
    "    ], \n",
    "    [\n",
    "        'TP_LC', \n",
    "        'TP_VC',\n",
    "        'FN',\n",
    "        'FP',\n",
    "        'TN',\n",
    "    ], \n",
    "    default='ERROR'\n",
    ")\n",
    "\n",
    "df_joined.loc[:, 'cluster'] = np.select(\n",
    "    [df_joined['cat'] == \"TP_LC\", df_joined['cat'] == \"TP_VC\", df_joined['cat'] == \"FN\", df_joined['cat'] == \"FP\", df_joined['cat'] == \"TN\",],\n",
    "    [0, 1, 2, 3, 4,],\n",
    "    default=-1\n",
    ")\n",
    "\n",
    "df_joined = df_joined.rename(columns={\"output\": \"pred_score\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sig(x):\n",
    "    return 1 / (1 + np.exp(-x))\n",
    "\n",
    "\n",
    "def invsig(y):\n",
    "    return np.log(y / (1 - y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot_shap_cols = [df_joined[[f\"{varname}_shapval_lvl{lvl}\" for lvl in range(64, 138)]].sum(axis=1) for varname in ccc.LVL_TRAIN_COLS]\n",
    "df_plot_shap = pd.concat(df_plot_shap_cols, axis=1)\n",
    "df_plot_shap.columns = [f\"{c}_shap\" for c in ccc.LVL_TRAIN_COLS]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot_meta_cols = [df_joined[[f\"{varname}_meta_lvl{lvl}\" for lvl in range(64, 138)]].sum(axis=1) for varname in ccc.LVL_TRAIN_COLS]\n",
    "df_plot_meta = pd.concat(df_plot_meta_cols, axis=1)\n",
    "df_plot_meta.columns = [f\"{c}_meta\" for c in ccc.LVL_TRAIN_COLS]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all = pd.concat([df_joined, df_plot_shap, df_plot_meta], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "invsig(used_threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for c in ccc.LVL_TRAIN_COLS:\n",
    "    df_all[f\"{c}_relative_shap\"] = df_all[f\"{c}_shap\"] / (invsig(used_threshold) - df_all[\"shap_base_value\"])\n",
    "    \n",
    "df_all[\"hour_relative_shap\"] = df_all[\"hour_shapval\"] / (invsig(used_threshold) - df_all[\"shap_base_value\"])\n",
    "df_all[\"dayofyear_relative_shap\"] = df_all[\"dayofyear_shapval\"] / (invsig(used_threshold) - df_all[\"shap_base_value\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = [f\"{c}_shap\" for c in ccc.LVL_TRAIN_COLS] + [f\"{c}_relative_shap\" for c in ccc.LVL_TRAIN_COLS + [\"hour\", \"dayofyear\"]] + [\"shap_base_value\", \"cat\", \"cluster\"] + ccc.INDEX_COLS\n",
    "\n",
    "df_TP = df_all[df_all['cat'].isin([\"TP_LC\", \"TP_VC\"])][cols]\n",
    "df_FP = df_all[df_all['cat'].isin([\"FP\"])][cols]\n",
    "df_FN = df_all[df_all['cat'].isin([\"FN\"])][cols]\n",
    "df_TN = df_all[df_all['cat'].isin([\"TN\"])][cols]\n",
    "\n",
    "df_TP_LC = df_TP.query('cat == \"TP_LC\"')\n",
    "df_TP_VC = df_TP.query('cat == \"TP_VC\"')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cloud_tp_sum = df_TP['cswc_relative_shap'] + df_TP['ciwc_relative_shap'] + df_TP['crwc_relative_shap'] + df_TP['clwc_relative_shap']\n",
    "mass_tp_sum = df_TP['q_relative_shap'] + df_TP['t_relative_shap']\n",
    "wind_tp_sum = df_TP['u_relative_shap'] + df_TP['v_relative_shap'] + df_TP['w_relative_shap']\n",
    "\n",
    "df_cloud = df_TP[cloud_tp_sum > 0.5]\n",
    "df_cloudmass = df_TP[(cloud_tp_sum > 0.5) & (mass_tp_sum > wind_tp_sum)]\n",
    "df_cloudwind = df_TP[(cloud_tp_sum > 0.5) & (mass_tp_sum <= wind_tp_sum)]\n",
    "\n",
    "df_mass = df_TP[mass_tp_sum > 0.5]\n",
    "df_wind = df_TP[wind_tp_sum > 0.5]\n",
    "\n",
    "df_nodom = df_TP[(cloud_tp_sum <= 0.5) & (mass_tp_sum <= 0.5) & (wind_tp_sum <= 0.5)]\n",
    "\n",
    "df_cloud_plus_TN = pd.concat([df_cloud, df_TN])\n",
    "df_cloudmass_plus_TN = pd.concat([df_cloudmass, df_TN])\n",
    "df_cloudwind_plus_TN = pd.concat([df_cloudwind, df_TN])\n",
    "\n",
    "df_mass_plus_TN = pd.concat([df_mass, df_TN])\n",
    "df_wind_plus_TN = pd.concat([df_wind, df_TN])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Number of samples in cloud-dominant TPs:\\t\\t{len(df_cloud)}\")\n",
    "print(f\"Number of samples in cloud-mass-dominant TPs:\\t\\t{len(df_cloudmass)}\")\n",
    "print(f\"Number of samples in cloud-wind-dominant TPs:\\t\\t{len(df_cloudwind)}\")\n",
    "print(f\"Number of samples in mass-dominant TPs: \\t\\t {len(df_mass)}\")\n",
    "print(f\"Number of samples in wind-dominant TPs: \\t\\t{len(df_wind)}\")\n",
    "print(f\"Number of samples TPs without dominance: \\t\\t{len(df_nodom)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Relative number of samples in cloud-dominant TPs:\\t\\t{100 * len(df_cloud) / len(df_TP)}\")\n",
    "print(f\"Relative number of samples in cloud-mass-dominant TPs:\\t\\t{100 * len(df_cloudmass) / len(df_TP)}\")\n",
    "print(f\"Relative number of samples in cloud-wind-dominant TPs:\\t\\t{100 * len(df_cloudwind) / len(df_TP)}\")\n",
    "print(f\"Relative number of samples in mass-dominant TPs: \\t\\t {100 * len(df_mass) / len(df_TP)}\")\n",
    "print(f\"Relative number of samples in wind-dominant TPs: \\t\\t{100 * len(df_wind) / len(df_TP)}\")\n",
    "print(f\"Relative number of samples TPs without dominance: \\t\\t{100 * len(df_nodom) / len(df_TP)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Number of samples being cloud and mass dominant at the same time:\\t\\t{len(df_cloud[['latitude']].join(df_mass[['longitude']], how='inner'))}\")\n",
    "print(f\"Number of samples being cloud and wind dominant at the same time:\\t\\t{len(df_cloud[['latitude']].join(df_wind[['longitude']], how='inner'))}\")\n",
    "print(f\"Number of samples being cloud, wind and mass dominant at the same time:\\t\\t{len(df_cloud[['latitude']].join(df_wind[['longitude']], how='inner').join(df_mass[['hour']], how='inner'))}\")\n",
    "print(f\"Number of samples being mass and wind dominant at the same time:\\t\\t{len(df_mass[['latitude']].join(df_wind[['longitude']], how='inner'))}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Used for debugging to quickly reload utils_plots.py\n",
    "import importlib\n",
    "importlib.reload(utils_plots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_clusters = {\n",
    "                0 : 'TP_LC',\n",
    "                1 : 'TP_VC',\n",
    "                2 : 'FN',\n",
    "                3 : 'FP',\n",
    "                4 : 'TN',\n",
    "}\n",
    "\n",
    "\n",
    "def get_color_palette(categories_to_plot=['TP_LC', 'TP_VC', 'FN', 'FP', 'TN']):    \n",
    "    return sns.color_palette([utils_plots.CLUSTER_COLORS[col] for col in categories_to_plot])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.boxplot(df_TP_LC[PLOT_COLS])\n",
    "plt.xticks(rotation=90)\n",
    "\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.boxplot(df_all, y=\"shap_base_value\", hue=\"cluster\", palette=get_color_palette())\n",
    "fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "for key in plot_clusters:\n",
    "    fig.legend_.texts[key].set_text(plot_clusters[key])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_bars_per_variable(df, title=\"\"):\n",
    "    f, axs = plt.subplots(len(PLOT_COLS), 1, constrained_layout=True, sharex=True)\n",
    "    f.set_figheight(50)\n",
    "    f.set_figwidth(15)\n",
    "\n",
    "    if (title != \"\"):\n",
    "        f.suptitle(title)\n",
    "\n",
    "    df.sort_values(by=\"cluster\", inplace=True)\n",
    "\n",
    "    avail_clusters = df[\"cluster\"].unique()\n",
    "    avail_categories = df[\"cat\"].unique()\n",
    "    \n",
    "    for idx, varname in enumerate(PLOT_COLS):\n",
    "        fig = sns.boxplot(df, ax=axs[idx], y=varname, hue=\"cluster\", palette=get_color_palette(avail_categories))\n",
    "        fig.set_ylim(-1, 1)\n",
    "        fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    \n",
    "        for idx, key in enumerate(avail_clusters):\n",
    "            fig.legend_.texts[idx].set_text(plot_clusters[key])\n",
    "\n",
    "    postfix = f\"_{title}\" if title != \"\" else \"\"\n",
    "    \n",
    "    f.savefig(os.path.join(\"tmp\", f\"boxplots{postfix}.png\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bars_per_variable(df_mass_plus_TN, \"mass_shapsum > 0.5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(len(PLOT_COLS), 1, sharex=True)\n",
    "f.set_figheight(50)\n",
    "f.set_figwidth(15)\n",
    "\n",
    "for idx, varname in enumerate(PLOT_COLS):\n",
    "    fig = sns.boxenplot(df_all, ax=axs[idx], y=varname, hue=\"cluster\", palette=get_color_palette())\n",
    "    fig.set_ylim(-1, 1)\n",
    "    fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    for key in plot_clusters:\n",
    "        fig.legend_.texts[key].set_text(plot_clusters[key])\n",
    "        \n",
    "f.savefig(os.path.join(\"tmp\", \"boxenplots.png\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(len(PLOT_COLS), 1, sharex=True)\n",
    "f.set_figheight(50)\n",
    "f.set_figwidth(15)\n",
    "\n",
    "df_filtered = df_all.query(\"cluster in [0.0, 1.0, 4.0]\").copy()\n",
    "df_filtered.loc[:, \"is_TP\"] = (df_filtered[\"cluster\"] <= 1)\n",
    "\n",
    "plot_clusters_filtered = {\n",
    "                            4: 'TN',\n",
    "                            1: 'TP_VC',\n",
    "}\n",
    "\n",
    "palette_filtered = sns.color_palette([utils_plots.CLUSTER_COLORS[plot_clusters_filtered[col]] for col in plot_clusters_filtered] if len(plot_clusters_filtered) > 0 else utils_plots.CLUSTER_COLORS.values())\n",
    "\n",
    "for idx, varname in enumerate(PLOT_COLS):\n",
    "    fig = sns.boxplot(df_filtered, ax=axs[idx], y=varname, hue=\"is_TP\", palette=get_color_palette(plot_clusters_filtered.values()))\n",
    "    fig.set_ylim(-1, 1)\n",
    "    fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    for idx, key in enumerate(plot_clusters_filtered):\n",
    "        fig.legend_.texts[idx].set_text(plot_clusters_filtered[key])\n",
    "        \n",
    "f.savefig(os.path.join(\"tmp\", \"boxplots_TP_vs_TN.png\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in PLOT_COLS:\n",
    "    tpm = np.quantile(df_TP[col], 0.5)\n",
    "    tnq = np.quantile(df_TN[col], 0.75)\n",
    "\n",
    "    print(f\"{col}: TP median = {tpm:0.2f} - TN 3rd quartile = {tnq:0.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in PLOT_COLS:\n",
    "    print(f\"{col}: {stats.ttest_ind(df_TP[col], df_TN[col], equal_var=False, alternative='greater')}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_histograms_per_variable(df, title=\"\"):\n",
    "    size_of_smallest_cl = df['cluster'].value_counts().min()\n",
    "    \n",
    "    f, axs = plt.subplots(len(PLOT_COLS), 1, constrained_layout=True, sharey=True)\n",
    "    f.set_figheight(60)\n",
    "    f.set_figwidth(15)\n",
    "    \n",
    "    if (title != \"\"):\n",
    "        f.suptitle(title)\n",
    "        \n",
    "    postfix = f\"_{title}\" if title != \"\" else \"\"\n",
    "    \n",
    "    for idx, varname in enumerate(PLOT_COLS):\n",
    "        cats = [\"TP_LC\", \"TP_VC\", \"TN\"]\n",
    "        \n",
    "        for cat in cats:\n",
    "            sns.histplot(df.query(f\"cat == '{cat}'\").sample(size_of_smallest_cl), ax=axs[idx], color=utils_plots.CLUSTER_COLORS[cat], x=varname)\n",
    "            \n",
    "        axs[idx].set_xlim([-1, 1])\n",
    "        \n",
    "        axs[idx].legend(labels=cats)\n",
    "    \n",
    "    f.savefig(os.path.join(\"tmp\", f\"histgram_relative{postfix}.png\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_histograms_per_variable(df_cloud_plus_TN, \"cloud\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_histograms_per_variable(df_mass_plus_TN, \"mass\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_histograms_per_variable(df_wind_plus_TN, \"wind\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_of_smallest_cl = min([len(df_TP_LC), len(df_TP_VC), len(df_TN)])\n",
    "\n",
    "ax = sns.histplot(df_TP_LC.sample(size_of_smallest_cl), x=\"shap_base_value\", color=utils_plots.CLUSTER_COLORS[\"TP_LC\"])\n",
    "sns.histplot(df_TP_VC.sample(size_of_smallest_cl), ax=ax, x=\"shap_base_value\", color=utils_plots.CLUSTER_COLORS[\"TP_VC\"])\n",
    "sns.histplot(df_TN.sample(size_of_smallest_cl), ax=ax, x=\"shap_base_value\", color=utils_plots.CLUSTER_COLORS[\"TN\"])\n",
    "    \n",
    "ax.legend(labels=[\"TP less confident\", \"TP very confident\", \"TN\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.get_figure().savefig(\"tmp/basevalues.png\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(df_TN.sample(size_of_smallest_cl), x=\"shap_base_value\", color=utils_plots.CLUSTER_COLORS[\"TN\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all[[\"longitude\", \"latitude\", \"shap_base_value\"]].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_basevalues = df_all.groupby([\"longitude\", \"latitude\"])[\"shap_base_value\"].mean().reset_index()\n",
    "\n",
    "m = Basemap(projection='lcc', resolution='f', lon_0=12.5, lat_0=47.5, llcrnrlon=8, llcrnrlat=45, urcrnrlon=17, urcrnrlat=50)\n",
    "m.drawmapboundary()\n",
    "m.drawcountries(linewidth=2)\n",
    "\n",
    "m.scatter(df_basevalues['longitude'], df_basevalues['latitude'], c=df_basevalues[\"shap_base_value\"], cmap=\"jet\", s=10, latlon=True)\n",
    "plt.colorbar(label='BASE_VALUE', extend=\"max\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_countlightningcells = df_all.query(\"target == 1.0\").groupby([\"longitude\", \"latitude\"])[\"target\"].count().reset_index()\n",
    "df_countlightningcells.rename(columns={\"target\": \"count\"}, inplace=True)\n",
    "\n",
    "m = Basemap(projection='lcc', resolution='f', lon_0=12.5, lat_0=47.5, llcrnrlon=8, llcrnrlat=45, urcrnrlon=17, urcrnrlat=50)\n",
    "m.drawmapboundary()\n",
    "m.drawcountries(linewidth=2)\n",
    "\n",
    "m.scatter(df_countlightningcells['longitude'], df_countlightningcells['latitude'], c=df_countlightningcells[\"count\"], cmap=\"jet\", s=10, latlon=True)\n",
    "plt.colorbar(label='NR_LIGHTNING_CELLS', extend=\"max\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_countvsbase = df_countlightningcells.merge(df_basevalues, how=\"outer\", on=[\"longitude\", \"latitude\"])\n",
    "\n",
    "sns.scatterplot(df_countvsbase, x=\"count\", y=\"shap_base_value\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_map(df, plotcols=PLOT_COLS, vmin=0, vmax=0.5):\n",
    "    for varname in plotcols:\n",
    "        df_var = df.groupby([\"longitude\", \"latitude\"])[varname].mean().reset_index()\n",
    "        \n",
    "        m = Basemap(projection='lcc', resolution='f', lon_0=12.5, lat_0=47.5, llcrnrlon=8, llcrnrlat=45, urcrnrlon=17, urcrnrlat=50)\n",
    "        m.drawmapboundary()\n",
    "        m.drawcountries(linewidth=2)\n",
    "        \n",
    "        m.scatter(df_var['longitude'], df_var['latitude'], c=df_var[varname], vmin=vmin, vmax=vmax, cmap=\"jet\", s=20, latlon=True)\n",
    "        plt.colorbar(label=varname, extend=\"max\")\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_map_count(df, suffix=\"\", df_norm=None, df_elevation=None, fixed_scale=True):\n",
    "    def roundup(x, divisor):\n",
    "        return math.ceil(x / divisor) * divisor\n",
    "    \n",
    "    varname = \"shap_base_value\"  # does not matter; we only count anyways\n",
    "    \n",
    "    df_var = df.groupby([\"longitude\", \"latitude\"])[varname].count().reset_index()\n",
    "\n",
    "    if df_norm is not None:\n",
    "        tmpcol = \"hour\"  # does not matter; we only count\n",
    "        df_norm_count = df_norm.groupby([\"longitude\", \"latitude\"])[tmpcol].count().reset_index()\n",
    "        df_var = df_var.merge(df_norm_count, on=[\"longitude\", \"latitude\"])\n",
    "\n",
    "        df_var[varname] /= df_var[tmpcol]\n",
    "        df_var[varname] *= 100\n",
    "\n",
    "        df_var.drop(tmpcol, axis=1, inplace=True) \n",
    "        \n",
    "    df_var[\"shap_base_value\"] -= 0.01 # hacky solution to ensure each sample is in the correct bin\n",
    "\n",
    "    if fixed_scale:\n",
    "        divisor = 5\n",
    "        \n",
    "        ticks = list(range(3, 30, divisor))\n",
    "        ticklabels = [f'{math.floor(n - 2)} - {math.floor(n + 2)}' for n in ticks]\n",
    "    else:\n",
    "        max_value = df_var[\"shap_base_value\"].max()\n",
    "    \n",
    "        if max_value > 50:\n",
    "            divisor = 20        \n",
    "            max_rounded = roundup(max_value, divisor)\n",
    "            \n",
    "            ticks = list(np.arange(10, max_rounded, divisor))\n",
    "            ticklabels = [f'{math.floor(n - 9)} - {math.floor(n + 10)}' for n in ticks]\n",
    "        elif max_value > 30:\n",
    "            divisor = 10        \n",
    "            max_rounded = roundup(max_value, divisor)\n",
    "            \n",
    "            ticks = list(np.arange(5.5, max_rounded, divisor))\n",
    "            ticklabels = [f'{math.floor(n - 4)} - {math.floor(n + 5)}' for n in ticks]\n",
    "        elif max_value >= 20:\n",
    "            divisor = 5        \n",
    "            max_rounded = roundup(max_value, divisor)\n",
    "            \n",
    "            ticks = list(range(3, max_rounded, divisor))\n",
    "            ticklabels = [f'{math.floor(n - 2)} - {math.floor(n + 2)}' for n in ticks]\n",
    "        elif max_value >= 10:\n",
    "            divisor = 3        \n",
    "            max_rounded = roundup(max_value, divisor)\n",
    "            \n",
    "            ticks = list(np.arange(1.5, max_rounded, divisor))\n",
    "            ticklabels = [f'{math.ceil(n - 1)} - {math.ceil(n + 1)}' for n in ticks]\n",
    "        else:\n",
    "            divisor = 2       \n",
    "            max_rounded = roundup(max_value, divisor)\n",
    "            \n",
    "            ticks = list(np.arange(1, max_rounded, divisor))\n",
    "            ticklabels = [f'{math.ceil(n - 0.5)} - {math.ceil(n + 0.5)}' for n in ticks]\n",
    "\n",
    "    \n",
    "    m = Basemap(projection='lcc', resolution='f', lon_0=12.5, lat_0=47.5, llcrnrlon=(8 - 0.125), llcrnrlat=(45 - 0.125), urcrnrlon=17.125, urcrnrlat=50.125)\n",
    "    m.drawmapboundary()\n",
    "    m.drawcountries(linewidth=2)\n",
    "\n",
    "    if df_elevation is not None:\n",
    "        m.scatter(df_elevation[\"longitude\"], df_elevation[\"latitude\"], c=df_elevation[\"elevation\"], latlon=True, cmap=cmap_map(lambda x: 0.5 + x / 1.7, cm.gray_r), alpha=0.4)\n",
    "\n",
    "    nr_categories = len(ticks) \n",
    "    cmap = mcolors.ListedColormap((\"#F4B8C0\", \"#E198B5\",\"#C87AAD\",\"#AA5FA5\",\"#87489D\",\"#5B3794\")).resampled(nr_categories)\n",
    "    \n",
    "    m.scatter(df_var['longitude'], df_var['latitude'], c=df_var[varname], cmap=cmap, vmin=0, vmax=(nr_categories * divisor), latlon=True, s=20, alpha=1)\n",
    "    \n",
    "    cbar = m.colorbar(ticks=ticks)\n",
    "    cbar.ax.set_yticklabels(ticklabels)\n",
    "\n",
    "    if suffix != \"\":\n",
    "        plt.savefig(os.path.join(\"tmp\", f\"mapplot_{suffix}.png\"), bbox_inches='tight')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map(df_cloud)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map(df_mass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map(df_wind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_elevation = utils.getElevationMap()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_cloud, \"cloud_percent\", df_TP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cloudmass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cloudwind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_cloudmass, \"cloudmass_percent\", df_TP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_cloudwind, \"cloudwind_percent\", df_TP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_mass, \"mass_percent\", df_TP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_wind, \"wind_percent\", df_TP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_TP, \"true_positives\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_cloud, \"cloud\", df_elevation=df_elevation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_cloudmass, \"cloud-mass\", df_elevation=df_elevation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_cloudwind, \"cloud-wind\", df_elevation=df_elevation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_mass, \"mass\", df_elevation=df_elevation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_map_count(df_wind, \"wind\", df_elevation=df_elevation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(len(PLOT_COLS), 1, sharey=True)\n",
    "f.set_figheight(40)\n",
    "f.set_figwidth(15)\n",
    "\n",
    "for idx, varname in enumerate(PLOT_COLS):\n",
    "    sns.histplot(df_TP_LC.sample(size_of_smallest_cl), ax=axs[idx], color=utils_plots.CLUSTER_COLORS[\"TP_LC\"], x=varname)\n",
    "    sns.histplot(df_TP_VC.sample(size_of_smallest_cl), ax=axs[idx], color=utils_plots.CLUSTER_COLORS[\"TP_VC\"], x=varname)\n",
    "    sns.histplot(df_TN.sample(size_of_smallest_cl), ax=axs[idx], color=utils_plots.CLUSTER_COLORS[\"TN\"], x=varname)\n",
    "    axs[idx].set_xlim([-1, 1])\n",
    "    \n",
    "    axs[idx].legend(labels=[\"TP less confident\", \"TP very confident\", \"TN\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axs = plt.subplots(len(ccc.LVL_TRAIN_COLS), 1, sharey=True)\n",
    "f.set_figheight(40)\n",
    "f.set_figwidth(15)\n",
    "\n",
    "for idx, varname in enumerate(ccc.LVL_TRAIN_COLS):\n",
    "    sns.histplot(df_filtered.query(\"cluster == 0\").sample(size_of_smallest_cl), ax=axs[idx], color=utils_plots.CLUSTER_COLORS[\"TP_LC\"], x=f\"{varname}_meta\")\n",
    "    sns.histplot(df_filtered.query(\"cluster == 1\").sample(size_of_smallest_cl), ax=axs[idx], color=utils_plots.CLUSTER_COLORS[\"TP_VC\"], x=f\"{varname}_meta\")\n",
    "    sns.histplot(df_filtered.query(\"cluster == 4\").sample(size_of_smallest_cl), ax=axs[idx], color=utils_plots.CLUSTER_COLORS[\"TN\"], x=f\"{varname}_meta\")\n",
    "    \n",
    "    axs[idx].legend(labels=[\"TP less confident\", \"TP very confident\", \"TN\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_stack(df, nr_samples=100):    \n",
    "    cols = [\n",
    "                'cswc',\n",
    "                'q',\n",
    "                'w',\n",
    "                'ciwc',\n",
    "                'hour',\n",
    "                'clwc',\n",
    "                'crwc',\n",
    "                't',\n",
    "                'u',\n",
    "                'v',\n",
    "    ]\n",
    "    \n",
    "    cols_r = [f\"{c}_relative_shap\" for c in cols]\n",
    "    \n",
    "    df_samp = df.sample(nr_samples)[cols_r]\n",
    "    df_samp.sort_values(by=cols_r, inplace=True, ascending=False)\n",
    "    \n",
    "    plt.figure(figsize=(20, 6))\n",
    "    plt.stackplot(range(nr_samples), df_samp[cols_r].T, labels=cols)\n",
    "    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.hlines(1, -1, nr_samples, color=\"red\", linestyle=\"dashed\")\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_stack_other(df, nr_samples=100):    \n",
    "    cols = [\n",
    "                'cswc',\n",
    "                'q',\n",
    "                'w',\n",
    "                'ciwc'\n",
    "    ]\n",
    "\n",
    "    other_cols =  [\n",
    "                    'clwc',\n",
    "                    'crwc',\n",
    "                    't',\n",
    "                    'u',\n",
    "                    'v',\n",
    "                    'hour',\n",
    "                    'dayofyear'\n",
    "    ]\n",
    "    \n",
    "    \n",
    "    cols_r = [f\"{c}_relative_shap\" for c in cols]\n",
    "    other_cols_r = [f\"{c}_relative_shap\" for c in other_cols]\n",
    "    \n",
    "    df_samp = df.sample(nr_samples)[cols_r + other_cols_r].copy()\n",
    "\n",
    "    df_samp[\"other\"] = df_samp[other_cols_r].sum(axis=1)\n",
    "    df_samp.sort_values(by=cols_r, inplace=True, ascending=False)\n",
    "    \n",
    "    plt.figure(figsize=(20, 6))\n",
    "    plt.stackplot(range(nr_samples), df_samp[cols_r + [\"other\"]].T, labels=cols + [\"other\"])\n",
    "    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.hlines(1, -1, nr_samples, color=\"red\", linestyle=\"dashed\")\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_stack_grouped(df, nr_samples=100):\n",
    "    df_samp = df.sample(nr_samples)\n",
    "\n",
    "    df_samp[\"cloud\"] = df_samp[[f\"{c}_relative_shap\" for c in ['cswc', 'ciwc', 'clwc', 'crwc']]].sum(axis=1)\n",
    "    df_samp[\"mass\"] = df_samp[[f\"{c}_relative_shap\" for c in ['q', 't']]].sum(axis=1)\n",
    "    df_samp[\"wind\"] = df_samp[[f\"{c}_relative_shap\" for c in ['u', 'v', 'w']]].sum(axis=1)\n",
    "    df_samp[\"time\"] = df_samp[[f\"{c}_relative_shap\" for c in ['hour', 'dayofyear']]].sum(axis=1)\n",
    "\n",
    "    cols = [\n",
    "                \"time\",\n",
    "                \"cloud\",\n",
    "                \"wind\",\n",
    "                \"mass\",\n",
    "    ]\n",
    "    \n",
    "    df_samp.sort_values(by=cols, inplace=True, ascending=False)\n",
    "    \n",
    "    plt.figure(figsize=(20, 6))\n",
    "    plt.stackplot(range(nr_samples), df_samp[cols].T, labels=cols)\n",
    "    plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    plt.hlines(1, -1, nr_samples, color=\"red\", linestyle=\"dashed\")\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_barplots_grouped(df):\n",
    "    df_tmp = df.copy()\n",
    "\n",
    "    df_tmp[\"cloud\"] = df_tmp[[f\"{c}_relative_shap\" for c in ['cswc', 'ciwc', 'clwc', 'crwc']]].sum(axis=1)\n",
    "    df_tmp[\"mass\"] = df_tmp[[f\"{c}_relative_shap\" for c in ['q', 't']]].sum(axis=1)\n",
    "    df_tmp[\"wind\"] = df_tmp[[f\"{c}_relative_shap\" for c in ['u', 'v', 'w']]].sum(axis=1)\n",
    "    df_tmp[\"time\"] = df_tmp[[f\"{c}_relative_shap\" for c in ['hour', 'dayofyear']]].sum(axis=1)\n",
    "\n",
    "    cols = [\n",
    "                \"time\",\n",
    "                \"cloud\",\n",
    "                \"mass\",\n",
    "                \"wind\",\n",
    "    ]\n",
    "\n",
    "    f, axs = plt.subplots(len(cols), 1, sharex=True)\n",
    "    f.set_figheight(50)\n",
    "    f.set_figwidth(15)\n",
    "    \n",
    "    for idx, varname in enumerate(cols):\n",
    "        fig = sns.boxplot(df_tmp, ax=axs[idx], y=varname, hue=\"cluster\", palette=get_color_palette())\n",
    "        fig.set_ylim(-1, 1)\n",
    "        fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "    \n",
    "        for key in plot_clusters:\n",
    "            fig.legend_.texts[key].set_text(plot_clusters[key])\n",
    "    \n",
    "    f.savefig(os.path.join(\"tmp\", \"boxplots_grouped.png\"), bbox_inches='tight')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_barplots_grouped(df_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stack_other(df_TP, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stack_grouped(df_TP, 300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stack_grouped(df_TP_LC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stack(df_TP_VC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stack(df_TN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 6))\n",
    "\n",
    "sns.boxplot(df_all.query(\"cluster <= 1.0\"), y=\"hour_relative_shap\", x=\"hour\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "f77099b03e26f5f2a39d729b0e5ff7cce62b3e0b846f69c8115517d984996de7"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}