{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d5cf2b3d-0f34-4d13-9ba3-dd2e3fb66337",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot 2D PDFs at four stations to compare: T, V, q, WDA, R3H\n",
    "\n",
    "from docx import Document\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mticker\n",
    "from itertools import chain\n",
    "from matplotlib.colors import LinearSegmentedColormap, Normalize\n",
    "import matplotlib.colors as mcolors\n",
    "from matplotlib.colors import BoundaryNorm\n",
    "import cmocean\n",
    "\n",
    "y_labs = ['T (˚C)', 'V (m/s)', 'q (g/kg)', 'WDA (degree)', 'Rain (mm/3hrs)'] \n",
    "x_labs = \"Local Time\"\n",
    "labs = ['(a)', '(b)', '(c)', '(d)', '(e)', '(f)', '(g)', '(h)', '(i)', '(j)', '(k)', '(l)', '(m)', '(n)', '(o)', '(p)',\n",
    "    '(q)', '(r)']\n",
    "\n",
    "cmap = LinearSegmentedColormap.from_list(\"white_to_dark\", [\"white\", \"black\"])\n",
    "norm = BoundaryNorm([0, 1, np.inf], ncolors=cmap.N, clip=True)\n",
    "\n",
    "count_std = np.zeros(6) \n",
    "\n",
    "y1 = [20, 0,  13,   0,  0]\n",
    "y2 = [25, 10, 25, 360, 30]\n",
    "yd = [1, 2, 3, 60, 5]\n",
    "\n",
    "vsh = ['T', 'V', 'q', 'WDA', 'R3h']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "51ce5ba3-dcbb-4870-a512-b2c338c0e4f1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/cdz/Desktop/Khafid_etal/Manuscript/2Dpdf_clusters3/pdf2d_data_cluster1_R3h.npz\n",
      "/Users/cdz/Desktop/Khafid_etal/Manuscript/2Dpdf_clusters3/pdf2d_data_cluster2_R3h.npz\n",
      "/Users/cdz/Desktop/Khafid_etal/Manuscript/2Dpdf_clusters3/pdf2d_data_cluster3_R3h.npz\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/MAAAGcCAYAAACVwfAFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHx0lEQVR4nO3deXhU5d3/8c9kD2ShAUISCTGWTQURAREUtyoaK4LYonUpKFq5WBRQ3JDNWlAfRbQIilaQxwfRWlHrngoCSqmC8KioCBqWCnmCiIkEss75/eEvpw6TQDJM5r5P8n5dVy6ZM9vHe+75zvmeM3OOz3EcRwAAAAAAwDOiTAcAAAAAAAANQzMPAAAAAIDH0MwDAAAAAOAxNPMAAAAAAHgMzTwAAAAAAB5DMw8AAAAAgMfQzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHtOsmvm9e/cqPT1d27Ztq9ftb731Vt10002NGwoADPt5bXzvvffk8/n0ww8/1Hn71157TT179pTf749cSACIsIasN7LOCMCEZtXMz5o1S4MGDdKxxx5br9vfdtttWrhwoQoKCho3GAAY1NDaePHFF8vn82nJkiWNGwwADGpIbWSdEYAJzaaZP3jwoP7yl7/o+uuvr/d90tPTNXDgQD3++OONmAwAzAmlNkrStddeqz//+c+NlAoAzGpobWSdEYAJzaaZf/PNNxUTE6N+/fpJkqqrqzVy5Ejl5uYqMTFRXbp00SOPPBJ0v0suuUTPPfdcpOMCQEQcWhtrfPDBB+rRo4cSEhLUt29fffrppwHXX3LJJfrwww/1zTffRDIuAEREbbVx06ZN+vWvf62UlBQlJydrwIAB+vrrr93rWWcEEGnNpplftWqVevfu7V72+/1q3769XnjhBX3++eeaOnWq7rrrLr3wwgsB9zv11FO1c+dObd++PdKRAaDRHVoba0yaNEkPPvigPvroI6Wnp+uSSy5RZWWle31OTo7S09O1evXqSMYFgIg4tDZ+++23OvPMM5WQkKDly5dr/fr1uu6661RVVeXehnVGAJEWYzpApGzbtk1ZWVnu5djYWM2YMcO9nJubqzVr1uiFF17QsGHD3OXHHHOMe/+cnJzIBQaACDi0NtaYNm2azj//fEnSM888o/bt22vZsmVB9bG+BxQFAC85tDY+9thjSk1N1dKlSxUbGytJ6ty5c8B9WGcEEGnNppk/ePCgEhISApY9/vjjeuqpp7R9+3YdPHhQFRUVOvnkkwNuk5iYKEk6cOBApKICQMTUVhslBXy1NC0tTV26dNEXX3wRcJvExERqI4Am6dDauHHjRg0YMMBt5GvDOiOASGs2X7Nv06aN9u3b515+4YUXNGHCBF133XV65513tHHjRl177bWqqKgIuN/3338vSWrbtm1E8wJAJBxaGw/H5/MFXP7++++pjQCapENrY02jfjisMwKItGbTzPfs2VOff/65e3n16tXq37+/Ro8erZ49e6pjx44BBzGp8dlnnyk2NlYnnnhiJOMCQEQcWhtrrF271v33vn379NVXX6lr167usrKyMn399dfq2bNnRHICQCQdWhtPOukkrV69OuDYIYdinRFApDWbZv6CCy7Qpk2b3K2sHTt21Lp16/T222/rq6++0pQpU/TRRx8F3W/16tUaMGBAvbbIAoDXHFoba9xzzz1699139dlnn2nEiBFq06aNhgwZ4l6/du1axcfHBx0FHwCagkNr49ixY1VSUqIrrrhC69at05YtW/Tf//3f2rx5s3sf1hkBRFqzaea7d++u3r17u0erHzVqlIYOHarLL79cffv21d69ezV69Oig+z333HO64YYbIh0XACLi0NpY47777tPNN9+sXr16affu3Xr11VcVFxfnXv/cc8/pqquuUosWLSIdGQAa3aG1sXXr1lq+fLn279+vs846S7169dKTTz4Z8Bt61hkBRJrPcRzHdIhIeeONN3Trrbfqs88+U1TUkbdjvP7665o0aZI++eQTxcQ0m2MFAmhmGlob9+zZo65du2rdunXKzc2NQEIAiLyG1EbWGQGY0KyqzUUXXaQtW7bo22+/VXZ29hFvX1paqoULF1KUATRpDa2NBQUFmjdvHo08gCatIbWRdUYAJjSrPfMAAAAAADQFzeY38wAAAAAANBU08wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAxzT5Q276/X7t2rVLycnJ8vl8puMAaIIcx9GPP/6orKysep3azQbURgCNjdoIAMHCWRubfDO/a9euep1qCQCO1s6dO9W+fXvTMeqF2gggUqiNABAsHLWxyTfzycnJkn4arJSUFMNpADRFJSUlys7OduuNF1AbATQ2aiMABAtnbWzyzXzNV6RSUlIoygAalZe+kkltBBAp1EYACBaO2uiNHzABAAAAAAAXzTwAAAAAAB5DMw8AAAAAgMfQzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTGmA0SK3++X3+83HQOH4fP5TEdw2TRXbMoiSVFR9mwDtGXO2PYaNURlZaUqKytNx5DjOKYjuGyZVzViYuz5qLbpdbKJbXPGltfJy7URdrNljtvIpvdddXW16QgB4uLiTEcIO3vWygEAAAAAQL3QzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAeQzMPAAAAAIDH0MwDAAAAAOAxMaYDRIrjOHIcx3QMHIbP5zMdwUp+v990hAC8TsG8XFt8Pp8Vr6lNYxgdHW06QgCbaoANcwVHZsucsSUHmh7balF1dbXpCK7KykrTEVyxsbGmIwSwZV0jnDnYMw8AAAAAgMfQzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAeQzMPAAAAAIDHxJgOECk+n08+n890DDmOYzpCABvGpIZNY2PTuERHR5uOgCbMcRwr3ntRUfZsW/b7/aYjBLCpBtgwV2rYNGdsGhc0LTbNLZvWjWwaF9vYVBttY8scDmcOXm0AAAAAADyGZh4AAAAAAI+hmQcAAAAAwGNo5gEAAAAA8BiaeQAAAAAAPIZmHgAAAAAAj6GZBwAAAADAY2jmAQAAAADwGKPN/KxZs9SnTx8lJycrPT1dQ4YM0ebNmwNuM2LECPl8voC/0047zVBiAGhc1EUACEZtBIBgRpv5lStXasyYMVq7dq3y8/NVVVWlgQMHqrS0NOB2F154oXbv3u3+vfHGG4YSA0Djoi4CQDBqIwAEizH55G+99VbA5YULFyo9PV3r16/XmWee6S6Pj49XRkbGUT1XdXW1qqurj+oxwiEqyq5fNjiOYzqCy+fzmY7g8vv9piO4bJi3P2fTHLZl/oZzvkSyLko/ZbdhvtvyWkpSdHS06QgBbBobsniDLXM4nDkiXRsdx7FijtlQn2vYMq8ku9YZJbvGpqyszHQEV2xsrOkITZ49a+WSiouLJUlpaWkBy9977z2lp6erc+fOuuGGG1RUVFTnY5SXl6ukpCTgDwC8Khx1UaI2AmhaqI0AIPkcGzY76qctoIMHD9a+ffu0evVqd/nzzz+vpKQk5eTkqKCgQFOmTFFVVZXWr1+v+Pj4oMeZPn26ZsyYEbT8u+++U0pKSqP+P9SHTXs1bWPTVlabtoSzZ95+JSUlatOmjYqLi8NaZ8JVF6W6a2NRUZEVtdGSjyJJdu1hkezKY9PrRC2yX0lJiVq1auXJ2vjDDz9YURttWh+xqRahbof+9MSkFi1amI4QwJZeo6SkRKmpqWGpjdY082PGjNHrr7+u999/X+3bt6/zdrt371ZOTo6WLl2qoUOHBl1fXl6u8vJy93JJSYmys7Np5j3AljeYZNeHJ828/RqrmQ9XXZTqro0088FsW2G1KY9NrxO1yH6N1cxHojbSzAezqRahbjTzdbOl1whnM2/0N/M1xo0bp1dffVWrVq06bFGWpMzMTOXk5GjLli21Xh8fH1/n1lcA8Ipw1kWJ2gigaaA2AsB/GG3mHcfRuHHjtGzZMr333nvKzc094n327t2rnTt3KjMzMwIJASCyqIsAEIzaCADBjH5HbcyYMXr22We1ZMkSJScnq7CwUIWFhTp48KAkaf/+/br11lv1z3/+U9u2bdN7772nQYMGqU2bNrr00ktNRgeARkFdBIBg1EYACGZ0z/z8+fMlSWeffXbA8oULF2rEiBGKjo7Wp59+qsWLF+uHH35QZmamzjnnHD3//PNKTk42kBgAGhd1EQCCURsBIJjxr9kfTmJiot5+++0IpQEA86iLABCM2ggAwTgULAAAAAAAHkMzDwAAAACAx1hxarpIiIqKsuKctDado1ey53yLkl1ZbHqdbMpiGxve05I9ORAevOe8gfNv1626utp0BEl2vUYN5ff7rchv29xC7Ww6t7tNp1q04T30c03x/cQaKAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAeQzMPAAAAAIDH0MwDAAAAAOAxNPMAAAAAAHgMzTwAAAAAAB5DMw8AAAAAgMfEmA4QKY7jyHEc0zEUFWXX9hMbxqSGTVl8Pp/pCC6/3286QgCbxsYWtr1GXlRdXW06gismptl8NDaYTXXaptfJthpgy+tkS45QREVFWbfOBnu1aNHCdAQr2fTZLtlTk8KZgyoFAAAAAIDH0MwDAAAAAOAxNPMAAAAAAHgMzTwAAAAAAB5DMw8AAAAAgMfQzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx8SYDhApjuPIcRzTMeT3+01HCGDDmNSIirJn25Jtr5NNbBobm+aMV1VVVamqqsp0DKtey+rqatMRUA821SKfz2c6QgBbxsaWHKHw+XxWvK42rafZxLa5ZdNnWGVlpekIrri4ONMRAjTF95M9Mw8AAAAAANQLzTwAAAAAAB5DMw8AAAAAgMfQzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHtPgU9Nt3rxZzz33nFavXq1t27bpwIEDatu2rXr27KkLLrhAl112meLj4xsjKwAAAAAAUAP2zG/YsEHnn3++evTooVWrVqlPnz4aP368/vjHP+rqq6+W4ziaPHmysrKydP/996u8vLwxcwMAAAAA0GzVe8/8kCFDNGnSJD3//PNKS0ur83b//Oc/9fDDD+uhhx7SXXfdFZaQAAAAAADgP+q9Z37Lli0aO3bsYRt5SerXr59eeOEF3XrrrUd8zFmzZqlPnz5KTk5Wenq6hgwZos2bNwfcxnEcTZ8+XVlZWUpMTNTZZ5+tTZs21Tc2AHgKdREAglEbASBYvffMx8XFHfb6H374Qa1atar37SVp5cqVGjNmjPr06aOqqipNnjxZAwcO1Oeff66WLVtKkh544AHNnj1bixYtUufOnXXvvffq/PPP1+bNm5WcnFzf+Nbw+/2mIwSIiuIYiLZzHMd0hAA2zRlb3k/hzBHpuhgVFWXVa2oDn89nOkIAm/LYVo9sYdu42PKeDmeOSNfG6upqVVdXhy1/qGx5LSW7apFN4yLJirlSgyx1s23ehIPPCeET6P7779exxx6ryy+/XJI0bNgw/e1vf1NGRobeeOMN9ejRI6Qwe/bsUXp6ulauXKkzzzxTjuMoKytL48eP1+233y5JKi8vV7t27XT//ffrxhtvPOJjlpSUKDU1Vd99951SUlJCyhVOfODXzaYPCZuKT2VlpekIAaKjo01HcNkyZ0pKSpSenq7i4uKw15nGqIs1mVNTU1VYWGhFbbSJTXNcsiuPTZ9hNn1+2caWDZ0lJSVq3bq1J2vj999/b0VttGme2/KZK9lViyTWG+tSn527kWTL+6mkpEStWrUKS20M6f/oiSeeUHZ2tiQpPz9f+fn5evPNN5WXl6dJkyaFHKa4uFiS3K/yFxQUqLCwUAMHDnRvEx8fr7POOktr1qwJ+XkAwCuoiwAQjNoIACGcmk6Sdu/e7Tbzr732moYNG6aBAwfq2GOPVd++fUMK4jiOJk6cqDPOOEPdunWTJBUWFkqS2rVrF3Dbdu3aafv27bU+Tnl5ecCR9EtKSkLKAwCmhasuStRGAE0HtREAfhLSnvlf/OIX2rlzpyTprbfe0nnnnSfpp+Ia6tdMxo4dq08++UTPPfdc0HWHfq3HcZw6v+oza9Yspaamun81Gx0AwGvCVRclaiOApoPaCAA/CamZHzp0qK688kqdf/752rt3r/Ly8iRJGzduVMeOHRv8eOPGjdOrr76qFStWqH379u7yjIwMSf/Z2lqjqKgoaMtrjTvvvFPFxcXuX81GBwDwknDWRYnaCKBpoDYCwH+E1Mw//PDDGjdunE444QTl5+crKSlJ0k9fvx89enS9H8dxHI0dO1YvvfSSli9frtzc3IDrc3NzlZGRofz8fHdZRUWFVq5cqf79+9f6mPHx8UpJSQn4AwCvaIy6KFEbAXgbtREAgjX4N/OVlZX6wx/+oClTpui4444LuG78+PENeqwxY8ZoyZIleuWVV5ScnOxuTU1NTVViYqJ8Pp/Gjx+vmTNnqlOnTurUqZNmzpypFi1a6Morr2xodACwHnURAIJRGwEgWIOb+djYWC1btkxTpkw56iefP3++JOnss88OWL5w4UKNGDFCknTbbbfp4MGDGj16tPbt26e+ffvqnXfe8eQ55gHgSKiLABCM2ggAwUI6z/y1116r7t27a+LEiY2RKaw4z/zh2XK+Rcmu85dyvtC62XTOa1vmTGOeZ76xcJ75utk0xyW78tj0GWbT55dtmsN55hsL55mvmy2fuZJdtUhivbEunGe+duE8z3xIp6br2LGj/vjHP2rNmjXq1auXWrZsGXD9TTfddFShAAAAAABA3UJq5p966im1atVK69ev1/r16wOu8/l8Vjbzfr/fmi3VNrFpy6ZNW3xtmis2vUZAY6mqqjIdwVo21UZ4gy2fYbbkCIXjOFZ8/vL+9wabPsNs2jNv2/yNjY01HUFSeGtjSM18QUFB2AIAAAAAAICGseOHAwAAAAAAoN5C2jNfXV2tRYsW6d1331VRUVHQVwWWL18elnAAAAAAACBYSM38zTffrEWLFunXv/61unXrZt3vIQAAAAAAaMpCauaXLl2qF154QRdddFG48wAAAAAAgCMI6TfzcXFx6tixY7izAAAAAACAegipmb/lllv0yCOPWHHKDgAAAAAAmpt6f81+6NChAZeXL1+uN998UyeeeGLQOfteeuml8KQDAAAAAABB6t3Mp6amBly+9NJLwx4GAAAAAAAcWb2b+YULFzZmDgAAAAAAUE8hHc3ei/x+v/x+v+kYio6ONh3BWhyDoXY2zNufq66uNh3BZcv7yctz15baeOjPtUyyaY5LdtWAmBh7Vhu8/L5rbLaMjS05QmFLbUTtOC123WxZN5Ls+syQ7BmbcOZo8AHwnnrqKQ0fPtzdU//888/r+OOP13HHHadp06aFLRgAAAAAAKhdgzaXzJkzR3fffbcuuOACTZ48Wbt27dLDDz+sCRMmyO/366GHHtIxxxyjP/zhD42VFwAAAACAZq9BzfwTTzyhBQsW6Morr9SGDRt06qmn6vHHH9fIkSMlSe3bt9djjz1GMw8AAAAAQCNq0Nfst2/frjPOOEOS1LNnT0VHR+u0005zrx8wYIC+/vrr8CYEAAAAAAABGtTMt2jRQqWlpe7ltm3bKikpKeA2VVVV4UkGAAAAAABq1aBmvmvXrvrkk0/cyzt37lROTo57+csvv9Sxxx4btnAAAAAAACBYg34zf//996tly5Z1Xr9jxw7deOONRx0KAAAAAADUrUHNfNu2bdW5c+c6rx89evRRBwIAAAAAAIfXoK/Z9+zZU8cff7xuv/12rVmzprEyAQAAAACAw2jQnvm9e/cqPz9fr7zyioYOHSrHcXTxxRdr8ODBGjhwoBISEhor51GLiopSVFSDtl00Cr/fbzpCABvGxEY+n890BGs5jmM6gsuW18mWHKGwpTbadPBU217P6Oho0xFcFRUVpiO4YmNjTUdw2fAe+jlb6rQtOUIRGxtrxRzz8hg2JtvqtA1zpUZZWZnpCC7besPq6mrTESSFN0eDPn0SEhI0aNAgPfXUU9q9e7eWLVumtm3b6o477lDr1q01ePBgPf300yoqKgpbQAAAAAAAECjkTck+n0/9+/fXfffdp88//1wbN27UmWeeqUWLFik7O1uPPfZYOHMCAAAAAID/r0Ffsz+cTp066ZZbbtEtt9yivXv36vvvvw/XQwMAAAAAgJ8JuZn/9ttv9cEHH6ioqCjgd+A+n0/jxo1T69atwxIQAAAAAAAECqmZX7hwoUaNGqW4uDi1bt064CAUNc08AAAAAABoHCE181OnTtXUqVN15513WncEVwAAAAAAmrqQOvEDBw7oiiuuoJEHAAAAAMCAkLrxkSNH6q9//Wu4swAAAAAAgHoI6Wv2s2bN0sUXX6y33npL3bt3V2xsbMD1s2fPDks4AAAAAAAQLKRmfubMmXr77bfVpUsXSQo6AB4AAAAAAGg8ITXzs2fP1tNPP60RI0aEOU7j8fv9AafQw08cxzEdwUo2zZXKykrTEQLYdKwMWzYeVlVVmY4Qsv3791sxjja955KTk01HCGBbDUCw6upq0xEQZgcOHFBMTMhncA6bhIQE0xGsFB0dbTpCANana1daWmo6QoC4uDjTESSF93M9pLXy+Ph4nX766Uf95KtWrdKgQYOUlZUln8+nl19+OeD6ESNGyOfzBfyddtppR/28AGAzaiMABKM2AkCgkJr5m2++WX/+85+P+slLS0vVo0cPzZ07t87bXHjhhdq9e7f798Ybbxz18wKAzaiNABCM2ggAgUL6/tCHH36o5cuX67XXXtOJJ54YdAC8l156qV6Pk5eXp7y8vMPeJj4+XhkZGaHEBABPojYCQDBqIwAECqmZb9WqlYYOHRruLLV67733lJ6erlatWumss87Sn/70J6Wnp9d5+/LycpWXl7uXS0pKIhETACKK2ggAwaiNAJqTkJr5hQsXhjtHrfLy8vTb3/5WOTk5Kigo0JQpU3Tuuedq/fr1io+Pr/U+s2bN0owZMyKSDwBMoDYCQDBqI4DmxudYcvhFn8+nZcuWaciQIXXeZvfu3crJydHSpUvr/GZAbVtYs7OzVVhYqJSUlHDH9jyOTF47m45MfODAAdMRAtg0Z2w5KmlJSYmysrJUXFwc9jrT2LVx69atVhy9naPZ182m95xNDv2Jn0mWrEq5bPkMKykpUbt27TxZG3ft2mXFeiNHs6+dbUezt+msNvv37zcdwWXb55dN643hqo0h7Znfu3evpk6dqhUrVqioqChoJez7778/qlB1yczMVE5OjrZs2VLnbeLj4+vc+goATRG1EQCCURsBNHUhNfNXX321vv76a40cOVLt2rWL2B7VvXv3aufOncrMzIzI8wGAF1AbASAYtRFAUxdSM//+++/r/fffV48ePY7qyffv36+tW7e6lwsKCrRx40alpaUpLS1N06dP12WXXabMzExt27ZNd911l9q0aaNLL730qJ4XAGxGbQSAYNRGAAgUUjPftWtXHTx48KiffN26dTrnnHPcyxMnTpQkDR8+XPPnz9enn36qxYsX64cfflBmZqbOOeccPf/889b9nhEAwonaCADBqI0AECikA+B99NFHuuOOOzR16lR169Yt6CA0NhwwpEZJSYlSU1OtOQCebQfJsenAFDZlqaioMB3BVVlZaTpCAJsOVGjTgUwyMzMb5SBPjaWmNu7YscOKzDYdAM+mA6tJdo2NLe85SYqJCWl/RLNgy+dGSUmJMjIyPFkb9+7da0Vmm+a5LQdWlOxaZ5TsOlixTa+TbQdwtOXzvaSkRK1atTJ3ALyaJz/33HMDljuOI5/PZ9UkAgAAAACgqQmpmb/qqqsUFxenJUuWRPQAeAAAAAAAIMRm/rPPPtOGDRvUpUuXcOcBAAAAAABHENIPTnr37q2dO3eGOwsAAAAAAKiHkPbMjxs3TjfffLMmTZqk7t27Bx1M4KSTTgpLOAAAAAAAECykZv7yyy+XJF133XXuMp/PxwHwAAAAAACIgJCa+YKCgnDnAAAAAAAA9RRSM5+TkxPuHAAAAAAAoJ7qfQC8f/7zn/V+0NLSUm3atCmkQAAAAAAA4PDqvWf+97//vY499ljdcMMNuuiii5SUlBR0m88//1zPPvusFi5cqAceeEAnnnhiWMMejaqqKlVVVZmOYR2fz2c6gsumYy3YlKWiosJ0hADR0dGmI7hsyWLTfGmoiooKlZeXm44RdCBVkyorK01HCBAVFdKJZxqF3+83HcFl0+sUExPSFx0bjeM4piNIsieHl9m07mrLZy4Or6yszHQEl23rR4mJiaYjSFJY17vq/enz+eef64knntDUqVN11VVXqXPnzsrKylJCQoL27dunL7/8UqWlpRo6dKjy8/PVrVu3sIUEAAAAAAD/Ue9mPjY2VmPHjtXYsWP18ccfa/Xq1dq2bZsOHjyoHj16aMKECTrnnHOUlpbWmHkBAAAAAGj2Qvpe2CmnnKJTTjkl3FkAAAAAAEA92PNDPAAAAAAAUC808wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMSEdzV6S3n33Xb377rsqKiqS3+8PuO7pp58+6mAAAAAAAKB2ITXzM2bM0D333KPevXsrMzNTPp8v3LkAAAAAAEAdQmrmH3/8cS1atEjXXHNNuPM0eWz4qFtUlD2/+nAcx3QEa9k0h215nWzJEYro6GjFxIT8Ja2wOXDggOkIrvj4eNMRAtiUx6b3vw3ztsah31A0zZY8tuQIxYEDB6yYYy1atDAdwWXT62nb5251dbXpCK64uDjTEVxJSUmmI1gpnLUlpO6poqJC/fv3D1sIAAAAAABQfyE189dff72WLFkS7iwAAAAAAKAeQtrHX1ZWpgULFugf//iHTjrpJMXGxgZcP3v27LCEAwAAAAAAwUJq5j/55BOdfPLJkqTPPvss4Dqbfk8HAAAAAEBTFFIzv2LFinDnAAAAAAAA9WTP4cMBAAAAAEC91HvP/NChQ7Vo0SKlpKRo6NChh73tSy+9dNTBAAAAAABA7erdzKemprq/h09NTW20QAAAAAAA4PDq3cwvXLiw1n8DAAAAAIDI4jfzAAAAAAB4TEhHs5ekF198US+88IJ27NihioqKgOs+/vjjow4WblVVVaqqqjIdQ9HR0aYjWMvv95uO4LIpiw3z1la2nArTy69RWVmZYmNjTcdQVJQ925Ztev9LUnV1tekILl4nb7BlbGzJEYrk5GQlJyebjmEVL7+ejc2mOm3LupEkVVZWmo4QwJb1tYMHD4btsUL6VH700Ud17bXXKj09XRs2bNCpp56q1q1b65tvvlFeXl7YwgEAAAAAgGAhNfPz5s3TggULNHfuXMXFxem2225Tfn6+brrpJhUXF4c7IwAAAAAA+JmQmvkdO3aof//+kqTExET9+OOPkqRrrrlGzz33XL0fZ9WqVRo0aJCysrLk8/n08ssvB1zvOI6mT5+urKwsJSYm6uyzz9amTZtCiQwAnkFtBIBg1EYACBRSM5+RkaG9e/dKknJycrR27VpJUkFBgRzHqffjlJaWqkePHpo7d26t1z/wwAOaPXu25s6dq48++kgZGRk6//zz3Y0HANAUURsBIBi1EQAChXQAvHPPPVd///vfdcopp2jkyJGaMGGCXnzxRa1bt05Dhw6t9+Pk5eXV+Rt7x3E0Z84cTZ482X3MZ555Ru3atdOSJUt04403hhIdAKxHbQSAYNRGAAgUUjO/YMEC94iWo0aNUlpamt5//30NGjRIo0aNCkuwgoICFRYWauDAge6y+Ph4nXXWWVqzZk2dRbm8vFzl5eXu5ZKSkrDkAQAbUBsBIBi1EUBzFNLX7KOiohQT85/tAMOGDdOjjz6qm266SXv27AlLsMLCQklSu3btApa3a9fOva42s2bNUmpqqvuXnZ0dljwAYANqIwAEozYCaI7CdsLYwsJCjRs3Th07dgzXQ0oKPlei4ziHPX/inXfeqeLiYvdv586dYc0DADagNgJAMGojgOakQc38Dz/8oKuuukpt27ZVVlaWHn30Ufn9fk2dOlXHHXec1q5dq6effjoswTIyMiQpaGtqUVFR0FbXn4uPj1dKSkrAHwA0FdRGAAhGbQTQHDWomb/rrru0atUqDR8+XGlpaZowYYIuvvhivf/++3rzzTf10Ucf6Xe/+11YguXm5iojI0P5+fnusoqKCq1cudI9LR4ANDfURgAIRm0E0Bw16AB4r7/+uhYuXKjzzjtPo0ePVseOHdW5c2fNmTMnpCffv3+/tm7d6l4uKCjQxo0blZaWpg4dOmj8+PGaOXOmOnXqpE6dOmnmzJlq0aKFrrzyygY/l+M4DTptXmOpqqoyHSGADWNS4+fHYTCtsrLSdATXwYMHTUcIEB8fbzqCKyoqbL8UOio1BwQNl0jWRlscOHDAdASXTXNcsqs2hnuuHw2bXqfDfY0b4RPJ2lhZWWnFuoAtn3OSXfPchtfm535+EEXTvv/+e9MRXMnJyaYjBEhNTTUdQVJ4+68GrSHs2rVLJ5xwgiTpuOOOU0JCgq6//vqQn3zdunU655xz3MsTJ06UJA0fPlyLFi3SbbfdpoMHD2r06NHat2+f+vbtq3feece6iQEA4URtBIBg1EYACORzGrBpIDo6WoWFhWrbtq2kn7a2fPLJJ8rNzW20gEerpKREqampKigo4HdQtWDPfO3KyspMR3D9+OOPpiMEsGlvWEJCgukIkn56jX75y1+quLjYM3WmpjZ+8cUXVqzo2vSes2mOS1KLFi1MR3CxZ752Nu2xlOx5nUpKSpSdne3J2rhnzx4rMrNnvnYVFRWmIwQoLi42HcHFnvm62bJnPpy1sUHdk+M4GjFihPsBWlZWplGjRqlly5YBt3vppZeOKhQAAAAAAKhbg5r54cOHB1y++uqrwxoGAAAAAAAcWYOa+YULFzZWDgAAAAAAUE/2/BAHAAAAAADUC808AAAAAAAeQzMPAAAAAIDH0MwDAAAAAOAx9pzYu5GVl5dbcR5jm86lLknR0dGmI1ipqqrKdARr2XL+YsmeLLbkCIXf77civw0Zath0LmVJqqysNB3BZdPY2DRnYmNjTUcIUFpaajqCJFmx3hWq8vJylZeXm45hFZvmuU21SJKiouzZP2rT6+Q4jukIAaqrq01HkBTeHPbMPAAAAAAAUC808wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAeQzMPAAAAAIDH0MwDAAAAAOAxNPMAAAAAAHhMjOkAkVJRUaGKigrTMVRZWWk6QoC4uDjTEVxVVVWmI7iqq6tNR3CVlJSYjhAgOjradARXy5YtTUeQJP3444+mI4SsurraivleWFhoOoIrKSnJdIQArVu3Nh3BZdNnWGxsrOkILpvGRZL8fr/pCJLsyRGKqqoqK9ZLbBrD+Ph40xFcNo2LJBUXF5uO4Nq+fbvpCK6ePXuajhDg4MGDpiNIksrKysL2WOyZBwAAAADAY2jmAQAAAADwGJp5AAAAAAA8hmYeAAAAAACPoZkHAAAAAMBjaOYBAAAAAPAYmnkAAAAAADyGZh4AAAAAAI+hmQcAAAAAwGNo5gEAAAAA8JgY0wEipaKiQuXl5aZjWMdxHNMRXLGxsaYjuGyaKwcPHjQdIUBiYqLpCC5bXqeKigrTEUL2/fffW5E/JsaejyPb3nM//vij6Qgum16nyspK0xFcNo2LJEVF2bGvxpYcoSgrK7NivcSmz1yb2LT+Ktk115OSkkxHcO3Zs8d0hADJycmmI0gK7/y1Z+YBAAAAAIB6oZkHAAAAAMBjaOYBAAAAAPAYq5v56dOny+fzBfxlZGSYjgUARlEbASAYtRFAc2PXEVtqceKJJ+of//iHezk6OtpgGgCwA7URAIJRGwE0J9Y38zExMWxVBYBDUBsBIBi1EUBzYvXX7CVpy5YtysrKUm5urq644gp98803piMBgHHURgAIRm0E0JxYvWe+b9++Wrx4sTp37qz/+7//07333qv+/ftr06ZNat26da33KS8vDzj3dElJSaTiAkBEUBsBIBi1EUBzY/We+by8PF122WXq3r27zjvvPL3++uuSpGeeeabO+8yaNUupqanuX3Z2dqTiAkBEUBsBIBi1EUBzY/We+UO1bNlS3bt315YtW+q8zZ133qmJEye6l0tKSpSdna2ysjLFxJj/33Ucx3SEAJWVlaYjuOLj401HcP3444+mI7j2799vOkKAiooK0xFciYmJpiNIkkpLS40+/9HUxgMHDsjn80Ui5mHt2bPHdASXLfOqRlSUPdvdbarTNmVJSkoyHSHAvn37TEeQZP7z62hrow3vPRsy1EhISDAdwVVVVWU6QgC/3286gsumdfutW7eajhCgY8eOpiNICm9tNN/dNkB5ebm++OILDRgwoM7bxMfHW/UBDwCNjdoIAMGojQCaOns299Xi1ltv1cqVK1VQUKB//etf+s1vfqOSkhINHz7cdDQAMIbaCADBqI0Amhur98z/+9//1u9+9zt99913atu2rU477TStXbtWOTk5pqMBgDHURgAIRm0E0NxY3cwvXbrUdAQAsA61EQCCURsBNDdWf80eAAAAAAAEo5kHAAAAAMBjaOYBAAAAAPAYmnkAAAAAADyGZh4AAAAAAI+hmQcAAAAAwGOsPjVdOG3dulUtWrQwHUMJCQmmIwTw+XymI7jS0tJMR3Dt2LHDdASXTVkkKTEx0XQEV0ZGhukIkqQDBw6YjhCyzZs3W/GaHjx40HQEV2VlpekIAdq0aWM6gsumLNXV1aYjuMrKykxHCLBnzx7TESRJpaWlpiOEbPfu3WrZsqXpGEpOTjYdwVVcXGw6gsum978krVmzxnQE16ZNm0xHcKWmppqOEGDXrl2mI0gK7zoPe+YBAAAAAPAYmnkAAAAAADyGZh4AAAAAAI+hmQcAAAAAwGNo5gEAAAAA8BiaeQAAAAAAPIZmHgAAAAAAj6GZBwAAAADAY2jmAQAAAADwGJp5AAAAAAA8JsZ0gEjZsGGD4uPjTcdQTIxdQ96qVSvTEVzJycmmI7i2b99uOoLr008/NR0hQHV1tekIrrZt25qOIEmqqKgwHSFkH374oeLi4kzH0FdffWU6gsuGz4qf++Uvf2k6gqtNmzamI7jy8vJMR3CVl5ebjhBg69atpiNIkg4ePGg6Qsg2b96sxMRE0zHkOI7pCK6cnBzTEVzp6emmIwTYs2eP6Qiu999/33QEV4sWLUxHCNC7d2/TESSF9zODPfMAAAAAAHgMzTwAAAAAAB5DMw8AAAAAgMfQzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAe43McxzEdojGVlJQoNTXVdAzUQ5s2bUxHcH333XemI8CDiouLlZKSYjpGvVAb4XVpaWmmI7h69OhhOkKAFStWmI4QgNqIcLvuuutMRwjw9NNPm44ADwpHbWTPPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAe44lmft68ecrNzVVCQoJ69eql1atXm44EAMZRGwEgGLURQHNhfTP//PPPa/z48Zo8ebI2bNigAQMGKC8vTzt27DAdDQCMoTYCQDBqI4DmxPpmfvbs2Ro5cqSuv/56HX/88ZozZ46ys7M1f/5809EAwBhqIwAEozYCaE5iTAc4nIqKCq1fv1533HFHwPKBAwdqzZo1td6nvLxc5eXl7uXi4uJGzYjw8fv9piMAR8VxnIg8D7URiNz7rT6qqqpMR7AatRHhVlFRYToCcNTCURutbua/++47VVdXq127dgHL27Vrp8LCwlrvM2vWLM2YMSMS8RBm33//vekIwFHZu3evUlNTG/15qI2AtG/fPtMRXPwm+/CojQi3Z5991nQE4KiFozZa3czX8Pl8AZcdxwlaVuPOO+/UxIkT3cs//PCDcnJytGPHjoh8kHhJSUmJsrOztXPnTqWkpJiOYw3GpW6MTe2Ki4vVoUMHpaWlRfR5qY3hxxyvG2NTO8albtTGpoN5XjfGpnaMS93CWRutbubbtGmj6OjooK2pRUVFQVtda8THxys+Pj5oeWpqKhOpDikpKYxNLRiXujE2tYuKisxhSKiNjY85XjfGpnaMS92ojU0H87xujE3tGJe6haM2Wn0AvLi4OPXq1Uv5+fkBy/Pz89W/f39DqQDALGojAASjNgJobqzeMy9JEydO1DXXXKPevXurX79+WrBggXbs2KFRo0aZjgYAxlAbASAYtRFAc2J9M3/55Zdr7969uueee7R7925169ZNb7zxhnJycup1//j4eE2bNq3Wr1A1d4xN7RiXujE2tTMxLtTGxsG41I2xqR3jUjdqY9PBuNSNsakd41K3cI6Nz7Hp3C4AAAAAAOCIrP7NPAAAAAAACEYzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4TJNu5ufNm6fc3FwlJCSoV69eWr16telIVvj222919dVXq3Xr1mrRooVOPvlkrV+/3nSsiFu1apUGDRqkrKws+Xw+vfzyy+51lZWVuv3229W9e3e1bNlSWVlZ+v3vf69du3aZCxwhhxsXSdq/f7/Gjh2r9u3bKzExUccff7zmz59vJmwEzZo1S3369FFycrLS09M1ZMgQbd68uc7b33jjjfL5fJozZ07kQtYTtbF21MafUBtrR22sHbWx6aM2UhcPh9pYu0jVxibbzD///PMaP368Jk+erA0bNmjAgAHKy8vTjh07TEczat++fTr99NMVGxurN998U59//rkeeughtWrVynS0iCstLVWPHj00d+7coOsOHDigjz/+WFOmTNHHH3+sl156SV999ZUuueQSA0kj63DjIkkTJkzQW2+9pWeffVZffPGFJkyYoHHjxumVV16JcNLIWrlypcaMGaO1a9cqPz9fVVVVGjhwoEpLS4Nu+/LLL+tf//qXsrKyDCQ9PGpj7aiN/0FtrB21sXbUxqaN2vgT6mLdqI21i1htdJqoU0891Rk1alTAsq5duzp33HGHoUR2uP32250zzjjDdAzrSHKWLVt22Nt8+OGHjiRn+/btkQllgdrG5cQTT3TuueeegGWnnHKKc/fdd0cwmXlFRUWOJGflypUBy//97387xxxzjPPZZ585OTk5zsMPP2wmYB2ojbWjNtaO2lg7amPdqI1NC7UxGHWxbtTGujVWbWySe+YrKiq0fv16DRw4MGD5wIEDtWbNGkOp7PDqq6+qd+/e+u1vf6v09HT17NlTTz75pOlYnlBcXCyfz9fstkYf6owzztCrr76qb7/9Vo7jaMWKFfrqq690wQUXmI4WUcXFxZKktLQ0d5nf79c111yjSZMm6cQTTzQVrU7UxrpRG0NHbfwJtfEn1MamhdoYGurif1Abf9JYtbFJNvPfffedqqur1a5du4Dl7dq1U2FhoaFUdvjmm280f/58derUSW+//bZGjRqlm266SYsXLzYdzWplZWW64447dOWVVyolJcV0HKMeffRRnXDCCWrfvr3i4uJ04YUXat68eTrjjDNMR4sYx3E0ceJEnXHGGerWrZu7/P7771dMTIxuuukmg+nqRm2sG7UxNNTG/6A2UhubImpjw1EXA1EbG7c2xoQjoK18Pl/AZcdxgpY1N36/X71799bMmTMlST179tSmTZs0f/58/f73vzeczk6VlZW64oor5Pf7NW/ePNNxjHv00Ue1du1avfrqq8rJydGqVas0evRoZWZm6rzzzjMdLyLGjh2rTz75RO+//767bP369XrkkUf08ccfW19nqI3BqI0NR20MRG2kNjZF1MaGoS4GozY2bm1sknvm27Rpo+jo6KCtqUVFRUFbXZubzMxMnXDCCQHLjj/++GZ/gJe6VFZWatiwYSooKFB+fn6z38J68OBB3XXXXZo9e7YGDRqkk046SWPHjtXll1+uBx980HS8iBg3bpxeffVVrVixQu3bt3eXr169WkVFRerQoYNiYmIUExOj7du365ZbbtGxxx5rLvDPUBvrRm1sGGpjIGojtbGpojbWH3UxGLWx8Wtjk2zm4+Li1KtXL+Xn5wcsz8/PV//+/Q2lssPpp58edFqEr776Sjk5OYYS2aumKG/ZskX/+Mc/1Lp1a9ORjKusrFRlZaWiogJLR3R0tPx+v6FUkeE4jsaOHauXXnpJy5cvV25ubsD111xzjT755BNt3LjR/cvKytKkSZP09ttvG0odiNpYN2pj/VEbg1EbqY1NFbWxfqiLtaM2Nn5tbLJfs584caKuueYa9e7dW/369dOCBQu0Y8cOjRo1ynQ0oyZMmKD+/ftr5syZGjZsmD788EMtWLBACxYsMB0t4vbv36+tW7e6lwsKCrRx40alpaUpKytLv/nNb/Txxx/rtddeU3V1tbvFPi0tTXFxcaZiN7rDjUuHDh101llnadKkSUpMTFROTo5WrlypxYsXa/bs2QZTN74xY8ZoyZIleuWVV5ScnOzOh9TUVCUmJqp169ZBH96xsbHKyMhQly5dTESuFbWxdtTG/6A21o7aWDtqY9NGbfwJdbFu1MbaRaw2HuVR9q322GOPOTk5OU5cXJxzyimnBJ0KoLn6+9//7nTr1s2Jj493unbt6ixYsMB0JCNWrFjhSAr6Gz58uFNQUFDrdZKcFStWmI7eqA43Lo7jOLt373ZGjBjhZGVlOQkJCU6XLl2chx56yPH7/WaDN7K65sPChQvrvI+Np19yHGpjXaiNP6E21o7aWDtqY9NHbaQuHg61sXaRqo2+//9kAAAAAADAI5rkb+YBAAAAAGjKaOYBAAAAAPAYmnkAAAAAADyGZh4AAAAAAI+hmQcAAAAAwGNo5gEAAAAA8BiaeQAAAAAAPIZmHgAAAAAAj6GZByRNnz5dJ5988lE9hs/n08svvxyWPABgA2ojAASiLsImNPOwwogRIzRkyBDTMWo1ffp0+Xy+w/5t27ZNu3fvVl5enum4AJoQaiMABKIuAv9BMw8cwa233qrdu3e7f+3bt9c999wTsCw7O1sZGRmKj483HRcAIoLaCACBqIuINJp5eMLKlSt16qmnKj4+XpmZmbrjjjtUVVXlXu/3+3X//ferY8eOio+PV4cOHfSnP/3Jvf72229X586d1aJFCx133HGaMmWKKisr6/XcSUlJysjIcP+io6OVnJwctOznX5natm2bfD6fXnjhBQ0YMECJiYnq06ePvvrqK3300Ufq3bu3kpKSdOGFF2rPnj0Bz7dw4UIdf/zxSkhIUNeuXTVv3ryjH0AATRK1EQACURfRnMSYDgAcybfffquLLrpII0aM0OLFi/Xll1/qhhtuUEJCgqZPny5JuvPOO/Xkk0/q4Ycf1hlnnKHdu3fryy+/dB8jOTlZixYtUlZWlj799FPdcMMNSk5O1m233dao2adNm6Y5c+aoQ4cOuu666/S73/1OKSkpeuSRR9SiRQsNGzZMU6dO1fz58yVJTz75pKZNm6a5c+eqZ8+e2rBhg2644Qa1bNlSw4cPb9SsALyF2khtBBCIukhdbHYcwALDhw93Bg8eXOt1d911l9OlSxfH7/e7yx577DEnKSnJqa6udkpKSpz4+HjnySefrPfzPfDAA06vXr3cy9OmTXN69OhRr/vm5OQ4Dz/8cNBySc6yZcscx3GcgoICR5Lz1FNPudc/99xzjiTn3XffdZfNmjXL6dKli3s5OzvbWbJkScDj/vGPf3T69etXr2wAmhZq40+ojQBqUBd/Ql2E4zgOe+ZhvS+++EL9+vWTz+dzl51++unav3+//v3vf6uwsFDl5eX61a9+VedjvPjii5ozZ462bt2q/fv3q6qqSikpKY2e/aSTTnL/3a5dO0lS9+7dA5YVFRVJkvbs2aOdO3dq5MiRuuGGG9zbVFVVKTU1tdGzAvAWaiO1EUAg6iJ1sbmhmYf1HMcJKMo1y6SfTu2RmJh42PuvXbtWV1xxhWbMmKELLrhAqampWrp0qR566KFGy1wjNjbW/XfN/8Ohy/x+vyS5/33yySfVt2/fgMeJjo5u7KgAPIbaSG0EEIi6SF1sbmjmYb0TTjhBf/vb3wIK9Jo1a5ScnKxjjjlGbdu2VWJiot59911df/31Qff/4IMPlJOTo8mTJ7vLtm/fHrH89dWuXTsdc8wx+uabb3TVVVeZjgPActRGAAhEXURzQzMPaxQXF2vjxo0By9LS0jR69GjNmTNH48aN09ixY7V582ZNmzZNEydOVFRUlBISEnT77bfrtttuU1xcnE4//XTt2bNHmzZt0siRI9WxY0ft2LFDS5cuVZ8+ffT6669r2bJlZv4nj2D69Om66aablJKSory8PJWXl2vdunXat2+fJk6caDoeAAOojdRGAIGoi9RF/IRmHtZ477331LNnz4Blw4cP16JFi/TGG29o0qRJ6tGjh9LS0jRy5Ejdfffd7u2mTJmimJgYTZ06Vbt27VJmZqZGjRolSRo8eLAmTJigsWPHqry8XL/+9a81ZcoU96imNrn++uvVokUL/dd//Zduu+02tWzZUt27d9f48eNNRwNgCLWR2gggEHWRuoif+JyaH5IAAAAAAABPiDIdAAAAAAAANAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMRzNHmhCysrKVFFRYToGAABAg8XFxSkhIcF0DMAzaOaBJqKsrEy5ubkqLCw0HQUAAKDBMjIyVFBQQEMP1BOnpgOaiJKSEqWmpmrnzp1KTk52l9e8xY/034bcNpT7HPrvxnjccNy2PvlM5j7SfU3mrS2rTWN3pDFsyOPaNP8iedsjPYbJ3PXJF6m8XnvvhjKGRzsept67jXlbm7KYHrtQ37uSVFxcrJSUlFqvAxCIPfNAE5OSkmJ9M9/UbmtTFhtua1MW/h8ZD8bD/G1tytLc/x9tZHM2wHYcAA8AAAAAAI+hmQcAAAAAwGNo5gEAAAAA8BiaeQAAAAAAPIZmHgAAAAAAj6GZBwAAAADAY2jmAQAAAADwGJp5AAAAAAA8hmYeAAAAAACPoZkHAAAAAMBjYkwHABBeJSUlchzHvVzz7yP9tyG3PdrHP9J9I5HhcLeJ9Hg0tbGzIcvhbnu45aZf88Pls3nsbHjN65MvUnm99t4NZQyPdjxMvXcb87Y2ZTE9dqG+dwE0DM080EQ4jqOkpCRlZ2ebjgIAANBgSUlJNPlAA9DMA02Ez+fT/v37tXPnTqWkpJiO02yUlJQoOzubcY8gxtwMxj3yGHMzGHczasbd5/OZjgJ4Bs080MSkpKSw8mEA4x55jLkZjHvkMeZmMO4AbMcB8AAAAAAA8BiaeQAAAAAAPIZmHmgi4uPjNW3aNMXHx5uO0qww7pHHmJvBuEceY24G424G4w40nM/hkJEAAAAAAHgKe+YBAAAAAPAYmnkAAAAAADyGZh4AAAAAAI+hmQeaiHnz5ik3N1cJCQnq1auXVq9ebTpSkzV9+nT5fL6Av4yMDNOxmpxVq1Zp0KBBysrKks/n08svvxxwveM4mj59urKyspSYmKizzz5bmzZtMhO2CTnSuI8YMSJo/p922mlmwjYRs2bNUp8+fZScnKz09HQNGTJEmzdvDrgN8z286jPmzPXwmz9/vk466SSlpKQoJSVF/fr105tvvulezzwHGoZmHmgCnn/+eY0fP16TJ0/Whg0bNGDAAOXl5WnHjh2mozVZJ554onbv3u3+ffrpp6YjNTmlpaXq0aOH5s6dW+v1DzzwgGbPnq25c+fqo48+UkZGhs4//3z9+OOPEU7atBxp3CXpwgsvDJj/b7zxRgQTNj0rV67UmDFjtHbtWuXn56uqqkoDBw5UaWmpexvme3jVZ8wl5nq4tW/fXvfdd5/WrVundevW6dxzz9XgwYPdhp15DjSQA8DzTj31VGfUqFEBy7p27erccccdhhI1bdOmTXN69OhhOkazIslZtmyZe9nv9zsZGRnOfffd5y4rKytzUlNTnccff9xAwqbp0HF3HMcZPny4M3jwYCN5mouioiJHkrNy5UrHcZjvkXDomDsOcz1SfvGLXzhPPfUU8xwIAXvmAY+rqKjQ+vXrNXDgwIDlAwcO1Jo1awylavq2bNmirKws5ebm6oorrtA333xjOlKzUlBQoMLCwoB5Hx8fr7POOot5HwHvvfee0tPT1blzZ91www0qKioyHalJKS4uliSlpaVJYr5HwqFjXoO53niqq6u1dOlSlZaWql+/fsxzIAQ084DHfffdd6qurla7du0Clrdr106FhYWGUjVtffv21eLFi/X222/rySefVGFhofr376+9e/eajtZs1Mxt5n3k5eXl6X/+53+0fPlyPfTQQ/roo4907rnnqry83HS0JsFxHE2cOFFnnHGGunXrJon53thqG3OJud5YPv30UyUlJSk+Pl6jRo3SsmXLdMIJJzDPgRDEmA4AIDx8Pl/AZcdxgpYhPPLy8tx/d+/eXf369dMvf/lLPfPMM5o4caLBZM0P8z7yLr/8cvff3bp1U+/evZWTk6PXX39dQ4cONZisaRg7dqw++eQTvf/++0HXMd8bR11jzlxvHF26dNHGjRv1ww8/6G9/+5uGDx+ulStXutczz4H6Y8884HFt2rRRdHR00FbroqKioK3baBwtW7ZU9+7dtWXLFtNRmo2aswcw783LzMxUTk4O8z8Mxo0bp1dffVUrVqxQ+/bt3eXM98ZT15jXhrkeHnFxcerYsaN69+6tWbNmqUePHnrkkUeY50AIaOYBj4uLi1OvXr2Un58fsDw/P1/9+/c3lKp5KS8v1xdffKHMzEzTUZqN3NxcZWRkBMz7iooKrVy5knkfYXv37tXOnTuZ/0fBcRyNHTtWL730kpYvX67c3NyA65nv4XekMa8Nc71xOI6j8vJy5jkQAr5mDzQBEydO1DXXXKPevXurX79+WrBggXbs2KFRo0aZjtYk3XrrrRo0aJA6dOigoqIi3XvvvSopKdHw4cNNR2tS9u/fr61bt7qXCwoKtHHjRqWlpalDhw4aP368Zs6cqU6dOqlTp06aOXOmWrRooSuvvNJgau873LinpaVp+vTpuuyyy5SZmalt27bprrvuUps2bXTppZcaTO1tY8aM0ZIlS/TKK68oOTnZ3TOZmpqqxMRE+Xw+5nuYHWnM9+/fz1xvBHfddZfy8vKUnZ2tH3/8UUuXLtV7772nt956i3kOhMLYcfQBhNVjjz3m5OTkOHFxcc4pp5wScHodhNfll1/uZGZmOrGxsU5WVpYzdOhQZ9OmTaZjNTkrVqxwJAX9DR8+3HGcn07XNW3aNCcjI8OJj493zjzzTOfTTz81G7oJONy4HzhwwBk4cKDTtm1bJzY21unQoYMzfPhwZ8eOHaZje1pt4y3JWbhwoXsb5nt4HWnMmeuN47rrrnPXVdq2bev86le/ct555x33euY50DA+x3GcSG48AAAAAAAAR4ffzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHkMzDwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAeQzMPAAAAAIDH0MwDAIAAO3fu1EUXXaSUlBQdf/zxeuWVV4Ju89e//lWtW7fWnj17DCQEAAA08wAAIMDw4cNVXFysF198UZdeeqmGDRumr7/+2r2+uLhYN998sx588EG1bdvWYFIAAJovn+M4jukQAADADgcOHFBSUpI++OAD9evXT5LUuXNnTZw4UaNGjZIkjRo1Sps3b9aKFStMRgUAoFmLMR0AAADYo6KiQo7jqGXLlu6ypKQklZWVSZLWrFmjxYsXa+PGjYYSAgAAia/ZAwCAn2nVqpW6du2qhx56SPv27dPLL7+s//3f/1X//v1VWVmpP/zhD7rzzjvVuXNn01EBAGjW2DMPAAAC/OUvf9Fll12mtLQ0RUVF6e6779app56qe++9V47j6PbbbzcdEQCAZo/fzAMAgCCVlZX65ptvlJGRodTUVG3ZskU9evTQO++8o/79+2vGjBlauHChKioq9Jvf/EYPPvigEhISTMcGAKDZoJkHAABH9Ktf/UodO3bUE088ob/85S+aMWOG3n33XSUlJenCCy/UkCFDNGPGDNMxAQBoNvjNPAAAOKxFixbp888/1/333y9JevPNN/Xb3/5WnTp1UmZmpkaOHKk33njDcEoAAJoXfjMPAADq9N133+nWW2/VvHnz1KpVK0mS4zgqLS11b7N//37xRT8AACKLZh4AANRp4sSJ6tu3r4YNG+Yuu+CCC3TLLbeoX79+SkpK0qOPPqrrr7/eYEoAAJoffjMPAABq9e6772rw4MHatGmTcnJy3OXV1dWaPHmynnnmGVVUVOiyyy7TI488osTERINpAQBoXmjmAQAAAADwGA6ABwAAAACAx9DMAwAAAADgMTTzAAAAAAB4DM08AAAAAAAeQzMPAAAAAIDH0MwDAAAAAOAxNPMAAAAAAHgMzTwAAAAAAB5DMw8AAAAAgMfQzAMAAAAA4DE08wAAAAAAeAzNPAAAAAAAHvP/ACq3dCSgDQmWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot four stations with different daily rainfall peak times\n",
    "\n",
    "cid = [1, 2, 3]\n",
    "\n",
    "nr = 1\n",
    "nc = 3\n",
    "row = 0\n",
    "nv = 4\n",
    "\n",
    "fig, axs = plt.subplots(nrows=nr, ncols=nc, figsize=(10, 4), layout=\"constrained\")\n",
    "axs_flat = axs.flatten()\n",
    "\n",
    "xbins = np.arange(0, 25, 1)     # hourly bins \n",
    "\n",
    "HC = np.zeros((24, 24))\n",
    "HCR = np.zeros((8, 24))\n",
    "\n",
    "vmin = 0\n",
    "vmax = 30\n",
    "\n",
    "ylimit1 = y1[nv]\n",
    "ylimit2 = y2[nv]\n",
    "dy = (ylimit2-ylimit1)/25\n",
    "ybins = np.arange(ylimit1, ylimit2, dy)\n",
    "\n",
    "# looo through clusters\n",
    "\n",
    "for cn in cid:\n",
    "    k = cn\n",
    "    ncl = k - 1\n",
    "    \n",
    "# loop through variables\n",
    "    \n",
    "    for v in vsh:\n",
    "\n",
    "        if v!= 'R3h':\n",
    "            continue\n",
    "\n",
    "        data_in_path = f\"/Users/cdz/Desktop/Khafid_etal/Manuscript/2Dpdf_clusters3/pdf2d_data_cluster{k}_{v}.npz\"\n",
    "        print(data_in_path)\n",
    "        \n",
    "        datain = np.load(data_in_path)\n",
    "        \n",
    "\n",
    "        \n",
    "        nlab = cn - 1\n",
    "    \n",
    "        axs[ncl].set_xlim(0, 24)\n",
    "        axs[ncl].set_xticks(np.arange(0, 25, 6))\n",
    "        \n",
    "        xedges = datain[\"xedges\"]\n",
    "        yedges = datain[\"yedges\"]\n",
    "\n",
    "        P = datain['H']\n",
    "        \n",
    "        xbins = np.arange(0, 25, 3)     # 3 hourly bins \n",
    "    \n",
    "        im = axs[ncl].pcolormesh(xbins, ybins, P.T, cmap=cmap, vmin = vmin, vmax = vmax)\n",
    "                 \n",
    "#    plt.colorbar(label=\"Counts\")\n",
    "        if k == 1:\n",
    "            axs[ncl].set_ylabel(y_labs[nv])\n",
    "        axs[ncl].set_ylim(ylimit1, ylimit2)\n",
    "        axs[ncl].set_yticks(np.arange(ylimit1, ylimit2, yd[nv]))\n",
    "        axs[ncl].set_title(labs[nlab], fontsize=10, loc='left')\n",
    "        \n",
    "        #    print(nv, mv, xlimit1, xlimit2, x_labs[nv])\n",
    "        \n",
    "        H = 0\n",
    "    \n",
    "\n",
    "    axs[ncl].set_xlabel(x_labs)\n",
    "        \n",
    "    ncl = ncl + 1  \n",
    "    datain = 0\n",
    "    \n",
    "cbar = fig.colorbar(im, ax=axs, orientation='horizontal', location='bottom',\n",
    "                    pad=0.01, shrink=0.6, aspect=40)\n",
    "cbar.set_label('%', fontsize=12)\n",
    "\n",
    "fig_path = f\"/Users/cdz/Desktop/Khafid_etal/Manuscript/2Dpdf_R3hcluster.png\"\n",
    "fig.savefig(fig_path)   \n",
    "plt.show()\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3b591d73-48c8-47c8-a240-f6b31fc37678",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5\n"
     ]
    }
   ],
   "source": [
    "print(nv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffe99be5-daad-4bb0-a337-fc7e447f5f86",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}