{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% Scratch Script for EWR Flume Sensitivity Testing\n"
    }
   },
   "outputs": [],
   "source": [
    "# Alexander Rey, 2021\n",
    "\n",
    "import os\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import netCDF4 as nc\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "from scipy.interpolate import LinearNDInterpolator, interp1d\n",
    "\n",
    "from dfm_tools.get_nc import get_netdata, get_ncmodeldata, plot_netmapdata\n",
    "from dfm_tools.get_nc_helpers import get_timesfromnc, get_ncfilelist, get_hisstationlist, get_ncvardimlist, get_timesfromnc\n",
    "import cartopy.crs as ccrs\n",
    "import contextily as ctx\n",
    "from dfm_tools.regulargrid import scatter_to_regulargrid\n",
    "import pathlib as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "pycharm": {
     "name": "#%% Read Model Log\n"
    }
   },
   "outputs": [],
   "source": [
    "# runLog = pd.read_excel(\"Y:/12828.101 English Wabigoon River/06_Models/00_Delft3D/ModelRuns.xlsx\", \"Sensitivity\")\n",
    "pth = pl.Path(\"//srv-ott3.baird.com/\", \"Projects\", \"12828.101 English Wabigoon River\", \"06_Models\", \"00_Delft3D\", \"ModelRuns.xlsx\")\n",
    "runLog = pd.read_excel(pth.as_posix(), \"Sensitivity\")\n",
    "\n",
    "dataPath = \"FlowFM_map.nc\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%% Import using DFM functions\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname found in long_name attribute\n",
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n"
     ]
    }
   ],
   "source": [
    "modelPlot  = [11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 24, 25, 26]\n",
    "\n",
    "dataOUT = dict()\n",
    "\n",
    "for i in modelPlot:\n",
    "    file_nc_hist = pl.Path(runLog['Run Location'][i],'Water_Quality','output','deltashell_his.nc')\n",
    "\n",
    "    # Print file variables\n",
    "    A = get_ncvardimlist(file_nc_hist.as_posix())\n",
    "    # Print file stations\n",
    "    # print(get_hisstationlist(file_nc_hist, 'Hg'))\n",
    "\n",
    "    # Define extraction variables\n",
    "    #dataVars = ['Hg', 'IM1', 'HgS1', 'FrHgIM1', 'DisHg', 'QHgSS']\n",
    "    dataVars = []\n",
    "    dataNames = []\n",
    "    dataLayer = []\n",
    "\n",
    "    for plottingVcount in range(0, 2):\n",
    "\n",
    "        allPlotVars  = runLog['Plotvars ' + str(plottingVcount+1)][i]\n",
    "\n",
    "        if not isinstance(allPlotVars, str):\n",
    "            continue\n",
    "\n",
    "        allPlotVars2 = allPlotVars.split(';')\n",
    "\n",
    "        for d in allPlotVars2:\n",
    "            if 'S1' in d:\n",
    "                dataVars.append(d[3:] + 'S1')\n",
    "                dataLayer.append(5)\n",
    "                dataNames.append(d[3:] + '_S1')\n",
    "            elif 'S2' in d:\n",
    "                dataVars.append(d[3:] + 'S2')\n",
    "                dataLayer.append(5)\n",
    "                dataNames.append(d[3:] + '_S2')\n",
    "            else:\n",
    "                dataVars.append(d[3:])\n",
    "                dataLayer.append(int(d[1]))\n",
    "                dataNames.append(d[3:] + '_L' + d[1])\n",
    "\n",
    "    for vIDX, v in enumerate(dataVars):\n",
    "        # Extract data from NetCDF file\n",
    "        dataIN = get_ncmodeldata(file_nc=file_nc_hist.as_posix(),\n",
    "                                 varname=v, timestep='all',\n",
    "                                 station=dataLayer[vIDX], silent=True)\n",
    "\n",
    "        # Convert to Numpy array\n",
    "        if vIDX == 0:\n",
    "            dataOUT[runLog['Run Name'][i]] = dict()\n",
    "\n",
    "        dataOUT[runLog['Run Name'][i]][dataNames[vIDX]] = np.array(dataIN)\n",
    "\n",
    "    # Get thicknesses\n",
    "    dataOUT[runLog['Run Name'][i]]['ActThS1'] = get_ncmodeldata(file_nc=file_nc_hist.as_posix(),\n",
    "                                                                varname='ActThS1', timestep='all',\n",
    "                                                                station=5, silent=True)\n",
    "\n",
    "    dataOUT[runLog['Run Name'][i]]['ActThS2'] = get_ncmodeldata(file_nc=file_nc_hist.as_posix(),\n",
    "                                                                varname='ActThS2', timestep='all',\n",
    "                                                                station=5, silent=True)\n",
    "\n",
    "    # Extract time\n",
    "    dataOUT[runLog['Run Name'][i]]['Times'] = get_timesfromnc(file_nc_hist)\n",
    "\n",
    "# Clean up input variables\n",
    "del dataIN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 32948 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 33788 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 27028 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 25296 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 20124 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 32772 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 32948 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 31888 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 38896 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 27568 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 38468 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 29908 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 29580 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 25180 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 36056 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 27576 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 26660 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 27460 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 24432 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x432 with 1 Axes>",
      "image/png": "[encoded data, 24424 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "modelPlot  = [11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 24, 25, 26]\n",
    "# modelPlot  = [11]\n",
    "\n",
    "for i in modelPlot:\n",
    "\n",
    "    for plottingVcount in range(0, 2):\n",
    "\n",
    "        allPlotVars  = runLog['Plotvars ' + str(plottingVcount+1)][i]\n",
    "\n",
    "        if not isinstance(allPlotVars, str):\n",
    "            continue\n",
    "\n",
    "        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 6))\n",
    "        fig.patch.set_facecolor('white')\n",
    "        allPlotVars2 = allPlotVars.split(';')\n",
    "\n",
    "        for d in allPlotVars2:\n",
    "            if 'S1' in d:\n",
    "                dataName = (d[3:] + '_S1')\n",
    "                layerDepth = dataOUT[runLog['Run Name'][i]]['ActThS1']\n",
    "            elif 'S2' in d:\n",
    "                dataName = (d[3:] + '_S2')\n",
    "                layerDepth = dataOUT[runLog['Run Name'][i]]['ActThS2']\n",
    "            else:\n",
    "                dataName = (d[3:] + '_L' + d[1])\n",
    "\n",
    "            if 'S1' in d or 'S2' in d:\n",
    "                axes.plot(dataOUT[runLog['Run Name'][i]]['Times'],\n",
    "                      dataOUT[runLog['Run Name'][i]][dataName]/layerDepth,\n",
    "                      label=dataName)\n",
    "                B = dataOUT[runLog['Run Name'][i]][dataName]\n",
    "            else:\n",
    "                axes.plot(dataOUT[runLog['Run Name'][i]]['Times'],\n",
    "                      dataOUT[runLog['Run Name'][i]][dataName],\n",
    "                      label=dataName)\n",
    "\n",
    "\n",
    "        fig.suptitle(runLog['Short Description'][i])\n",
    "        plotCountInt = plottingVcount+1\n",
    "        plt.title(runLog['Plot Title ' + str(plotCountInt)][i])\n",
    "        plt.ylabel(runLog['PlotUnits ' + str(plotCountInt)][i])\n",
    "        axes.legend(loc='upper right')\n",
    "        fig.autofmt_xdate()\n",
    "\n",
    "\n",
    "        fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/Modelling/'+\n",
    "                        'Figures/ProcessTest/Round1/' +\n",
    "                        runLog['Short Description'][i] + '_' +\n",
    "                        runLog['Plot Title ' + str(plotCountInt)][i] +\n",
    "                        '.png', bbox_inches='tight', dpi=200)\n",
    "\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Plot Sensitivity Data\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "pycharm": {
     "name": "#%% Plot Sensitivity Data\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 124772 characters]\n",
      "text/plain": [
       "<Figure size 648x1080 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(9, 15))\n",
    "fig.patch.set_facecolor('white')\n",
    "fig.tight_layout(pad=3)\n",
    "ax = axes.flat\n",
    "\n",
    "cmap = cm.get_cmap('Paired')\n",
    "\n",
    "dataVars        = ['Hg', 'IM1', 'HgS1', 'FrHgIM1', 'DisHg', 'QHgSS']\n",
    "dataVarsLabel   = ['Hg (g/m3)', 'IM1 (gDM/m3)', 'HgS1 (gDM/m2)', 'FrHgIM1 (-)', 'DisHg (gHg/m3)', 'QHgSS (mgHg/kgDM)']\n",
    "dataVarsLegend  = ['Mercury', ' Inorganic matter in water', 'Inorganic matter in layer S1',\n",
    "                   'Fraction Hg adsorbed to inorganic matter', 'Free dissolved Hg in water column',\n",
    "                   'Overall quality suspended solid for Hg']\n",
    "for vIDX, v in enumerate(dataVars):\n",
    "\n",
    "    # Loop through sensitivity runs\n",
    "    for i in modelPlot:\n",
    "        if i == 0:\n",
    "            ax[vIDX].plot(dataTime, dataOUT[v][:,i], color='k', linewidth=3)\n",
    "        else:\n",
    "            ax[vIDX].plot(dataTime, dataOUT[v][:,i], color=cmap.colors[i-1], linewidth=1)\n",
    "\n",
    "    # Plot settings\n",
    "    ax[vIDX].set_ylabel(dataVarsLabel[vIDX])\n",
    "    ax[vIDX].set_title(dataVarsLegend[vIDX])\n",
    "\n",
    "    if vIDX == 0:\n",
    "        ax[vIDX].legend(runLog['Short Description'], loc='upper right')\n",
    "\n",
    "plt.gcf().autofmt_xdate()\n",
    "# plt.gcf().align_ylabels()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/Modelling/FlumeTest/SensTest/Figures/SensTestParams.pdf',\n",
    "             bbox_inches='tight', dpi=200)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PyCharm (NTC_DFM)",
   "language": "python",
   "name": "pycharm-b2185061"
  },
  "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}