Last active
July 23, 2020 12:48
-
-
Save sfinkens/d76b153204236494fa0c2c3c8749a2f4 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"pygac azimuth angles\n", | |
"====================" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"%matplotlib inline" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import netCDF4\n", | |
"import h5py\n", | |
"import xarray as xr\n", | |
"import matplotlib.pyplot as plt\n", | |
"from pygac.utils import centered_modulus\n", | |
"import numpy as np" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"def comp(data):\n", | |
" return ((data - 180.0) / 0.01).astype('int16')\n", | |
"\n", | |
"\n", | |
"def decomp(data, offset):\n", | |
" data = data.astype(float)\n", | |
" data[data == -32001] = np.nan\n", | |
" return data * 0.01 + offset" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"There are two problems with the versions of pygac which compute azimuth angles in [0, 360]:\n", | |
"\n", | |
"1. The offset attribute in `/image5/what` (180) doesn't match the actual offset that was applied (0).\n", | |
"2. Azimuth values > 327.67 degrees are clipped because they don't fit into 16 bit integer." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 18, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAIICAYAAAB+XWiFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXyV5Zn/8e+dk30P2dhJBGSHAGERRRQX1GJbrHWp1qotLtXWWqul02XsdNO2/qatzviq+zajtqUuVWvVqoMLVVAQkEW2AGEN2ch+cs65f3882YCAICd5kpvPe17nlZPnbNfB8fHbm+u5bmOtFQAAAIDDi/G7AAAAAKA3IDgDAAAAR4DgDAAAABwBgjMAAABwBAjOAAAAwBEgOAMAAABHINbvAo5UTk6OLSgo8LsMAAAAOOyDDz7Ya63N7eyxXhOcCwoKtHTpUr/LAAAAgMOMMVsO9RitGgAAAMARIDgDAAAARyAqwdkY85AxZo8xZlWHY7cbY7YbY5a33M7r8NgPjDEbjDHrjDFzolEDAAAA0JWi1eP8iKR7JD12wPH/tNb+tuMBY8xoSZdIGiOpv6TXjDEnWmvDUaoFAACgR2hublZpaakaGxv9LgUHSExM1MCBAxUXF3fEr4lKcLbWLjLGFBzh078g6SlrbZOkzcaYDZKmSlocjVoAAAB6itLSUqWlpamgoEDGGL/LQQtrrcrLy1VaWqrCwsIjfl1X9zjfaIxZ0dLKkdVybICkbR2eU9py7CDGmGuMMUuNMUvLysq6uFQAAIDoamxsVHZ2NqG5hzHGKDs7+6j/JqArg/O9koZKKpK0U9JdR/sG1tr7rLXF1tri3NxOx+kBAAD0aITmnumz/HPpsuBsrd1trQ1bayOS7pfXjiFJ2yUN6vDUgS3HAAAAgB6ry4KzMaZfh1/nSWqduPG8pEuMMQnGmEJJwyW931V1AAAAANEQrXF0T8q7uG+EMabUGPN1Sb82xqw0xqyQdLqkmyXJWvuxpD9JWi3pZUk3MFEDAACga5SUlGjUqFGaP3++xowZo7PPPlsNDQ3asGGDzjzzTE2YMEGTJk3Sxo0b/S61x4vWVI1LOzn84GGe/wtJv4jGZwMAAPQGd75/p9ZWrI3qe47sM1Lfn/r9T33e+vXr9eSTT+r+++/XRRddpIULF+ruu+/WggULNG/ePDU2NioSiUS1NhdFa44zAAAAeqjCwkIVFRVJkiZPnqzNmzdr+/btmjdvniRvpjE+HcEZAACgGxzJynBXSUhIaLsfCARUVVXlWy29WVfPcQYAAEAPk5aWpoEDB+rZZ5+VJDU1Nam+vt7nqno+gjMAAMBx6PHHH9cf/vAHjR8/XjNmzNCuXbv8LqnHo1UDAADAYQUFBVq1alXb79/73vfa7r/++ut+lNRrseIMAAAAHAGCMwAAAHAECM4AAADAESA4A4DD/u2tf9PTa5/2uwwAcALBGQAc9t6u97S8bLnfZQCAEwjOAOC4pnCT3yUAgBMIzgDgMktwBo5nJSUlGjt2bKeP3X777RowYICKiopUVFSkl156qZur632Y4wwADrOyBGcAh3TzzTfvN9c5mkKhkGJjYw/5e2/EijMAOMzKKhgO+l0GAB+Fw2HNnz9fY8aM0dlnn62Ghoajev2dd96pcePGacKECVqwYIEkafny5Zo+fbrGjx+vefPmqbKyUpJ02mmn6Tvf+Y6Ki4v1+9///qDfO6qrq9PVV1+tqVOnauLEiXruueckSY888oguuOACnXPOORo+fLhuu+22tte8/PLLmjRpkiZMmKAzzjjjWP5YPpPeHfsBAIdlLSvOQI/y8OcOPjbmi9LU+VKwXvqfLx/8eNFXpImXSXXl0p+u2P+xq1781I9cv369nnzySd1///266KKLtHDhQl1++eWSpHvuuUePPfaYiouLdddddykrK2u/1/7973/Xc889p/fee0/JycmqqKiQJF1xxRW6++67NWvWLP3kJz/RT3/6U/3ud7+TJAWDQS1dulSS9Le//W2/3zv6xS9+odmzZ+uhhx5SVVWVpk6dqjPPPFOSF8yXLVumhIQEjRgxQt/61reUmJio+fPna9GiRSosLGyrpTux4gwADmPFGUBhYaGKiookSZMnT1ZJSYkk6frrr9fGjRu1fPly9evXT7fccstBr33ttdd01VVXKTk5WZLUp08fVVdXq6qqSrNmzZIkfe1rX9OiRYvaXnPxxRfv9x4H/t7qlVde0R133KGioiKddtppamxs1NatWyVJZ5xxhjIyMpSYmKjRo0dry5Yt+te//qVTTz1VhYWFbbV0N1acAcBxrDgDPcjhVojjkw//eEr2Ea0wHyghIaHtfiAQaGvVyM/Pbzs+f/58zZ0796jfuzMpKSmH/b2VtVYLFy7UiBEj9jv+3nvvHVRzKBSKSm3HihVnAHAcwRlAZ3bu3Nl2/5lnnul0+sZZZ52lhx9+WPX19ZKkiooKZWRkKCsrS2+99ZYk6fHHH29bfT4ac+bM0d133y1rrSRp2bJlh33+9OnTtWjRIm3evLmtlu7GijMAOMxaWjUAdO62227T8uXLZYxRQUGB/vjHPx70nHPOOUfLly9XcXGx4uPjdd555+mXv/ylHn30UV133XWqr6/XCSecoIcffvioP//HP/6xvvOd72j8+PGKRCIqLCzUCy+8cMjn5+bm6r777tMFF1ygSCSivLw8vfrqq0f9ucfCtKb8nq64uNh21lgOADi0mU/NVFO4Se9f9r7fpQDHpTVr1mjUqFF+l4FD6OyfjzHmA2ttcWfPp1UDABzGxYEAED0EZwBwmLVWYRtWKNIzLqwBgN6M4AwADrPy2vFYdQaAY0dwBgCXtVzGwmQNADh2BGcAOA4QnAHg2BGcAcBhra0aBGcAOHYEZwBwGMEZOL6VlJR0urGJJN1+++0aMGCAioqKVFRUpJdeeqmbq+t92AAFABzWOqufiwMBdObmm2/W9773vai8l7VW1lrFxLi7LuvuNwMAsOIMQOFwWPPnz9eYMWN09tlnq6Gh4Yhfe8MNN+j555+XJM2bN09XX321JOmhhx7SD3/4Q5WUlGjEiBG64oorNHbsWG3btk233nqrxo4dq3Hjxunpp5+WJL355ps67bTTdOGFF2rkyJG67LLL2v6H/UsvvaSRI0dq8uTJ+va3v625c+d2+h1uvfVWTZkyRePHj2/b5fBw77tkyRLNmDFDEyZM0NSpU1VTU/PZ/xBbsOIMAMcBgjPQQ5x22sHHLrpI+uY3pfp66bzzDn78yiu929690oUX7v/Ym29+6keuX79eTz75pO6//35ddNFFWrhwoS6//HJJ0j333KPHHntMxcXFuuuuu5SVlbXfa2fOnKm33npLn//857V9+3bt3LlTkvTWW2/pkksuaXv/Rx99VNOnT9fChQu1fPlyffTRR9q7d6+mTJmiU089VZK0bNkyffzxx+rfv79OPvlkvfPOOyouLta1116rRYsWqbCwUJdeemmn3+HBBx9URkaGlixZoqamJp188sk6++yzD/m+U6dO1cUXX6ynn35aU6ZM0b59+5SUlPSpf1afhhVnAHAYrRoACgsLVVRUJEmaPHmySkpKJEnXX3+9Nm7cqOXLl6tfv3665ZZbDnpta3BevXq1Ro8erfz8fO3cuVOLFy/WjBkzJElDhgzR9OnTJUlvv/22Lr30UgUCAeXn52vWrFlasmSJJGnq1KkaOHCgYmJiVFRUpJKSEq1du1YnnHCCCgsLJemQwfmVV17RY489pqKiIk2bNk3l5eVav379Id933bp16tevn6ZMmSJJSk9PV2zssa8Xs+IMAMcBVpyBHuJwK8TJyYd/PCfniFaYD5SQkNB2PxAItLVq5Ofntx2fP39+py0SAwYMUFVVlV5++WWdeuqpqqio0J/+9CelpqYqLS1N5eXlSklJ+Ux1hEJHvqOptVZ333235syZs9/xN99885je92ix4gwADmPnQACH0tp2IUnPPPPMIadvTJ8+Xb/73e906qmnaubMmfrtb3+rmTNndvrcmTNn6umnn1Y4HFZZWZkWLVqkqVOnHrKGESNGaNOmTW2r4K090QeaM2eO7r33XjU3N0uSPvnkE9XV1R32fXfu3Nm22l1TUxOVQM2KMwA4rLVVgxVnAAe67bbbtHz5chljVFBQ0HbB3YFmzpypV155RcOGDdOQIUNUUVFxyOA8b948LV68WBMmTJAxRr/+9a/Vt29frV27ttPnJyUl6b//+791zjnnKCUlpa214kDf+MY3VFJSokmTJslaq9zcXD377LOH/G7x8fF6+umn9a1vfUsNDQ1KSkrSa6+9ptTU1E/5Uzk803pS7emKi4vt0qVL/S4DAHqVSY9PUnOkWQumLtBloy7zuxzguLNmzRqNGjXK7zJ6tNraWqWmpspaqxtuuEHDhw/XzTff3C2f3dk/H2PMB9ba4s6eT6sGADiMVg0APd3999+voqIijRkzRtXV1br22mv9LumQaNUAAJe1/KUirRoAeqqbb76521aYjxUrzgDgMDZAAYDoITgDwHGA4Az4p7dcT3a8+Sz/XAjOAOAwepwBfyUmJqq8vJzw3MNYa1VeXq7ExMSjeh09zgDgMMbRAf4aOHCgSktLVVZW5ncpOEBiYqIGDhx4VK8hOAOAw+hxBvwVFxfXtp00ej9aNQDgOECrBgAcO4IzADiqY08lK84AcOwIzgDgqNY2DYkVZwCIBoIzABwHWHEGgGNHcAYAR3Vs1WDFGQCOHcEZABzVsVWDFWcAOHYEZwBwFMEZAKKL4AwAruqwURmtGgBw7AjOAOCojivOjeFGHysBADcQnAHAUYyjA4DoIjgDgONiY2LVFG7ab8oGAODoEZwBwFGtQTkxkChJao40+1kOAPR6BGcAcFRrq0Z8IF4SkzUA4FgRnAHAUQeuOBOcAeDYEJwBwHEp8SmSpLrmOp8rAYDejeAMAI5qbdXISsiSJFU1VflZDgD0egRnAHBUa6tGZkKmJKm6qdrPcgCg1yM4A4DjshK9FefKxkqfKwGA3o3gDACOamvVSKRVAwCigeAMAI5qDc6pcakKmADBGQCOUVSCszHmIWPMHmPMqg7H+hhjXjXGrG/5mdVy3Bhj/mCM2WCMWWGMmRSNGgAA+2vtcY4xMcpIyCA4A8AxitaK8yOSzjng2AJJ/7TWDpf0z5bfJelcScNbbtdIujdKNQAAOmFklJWQpapGgjMAHIuoBGdr7SJJFQcc/oKkR1vuPyrpix2OP2Y9/5KUaYzpF406AADtWlecjTGsOANAFHRlj3O+tXZny/1dkvJb7g+QtK3D80pbjgEAoqi1x1nyLhAkOAPoDXbV7dKvl/xaS3ct9buUg8R2x4dYa60xxn76M/dnjLlGXjuHBg8eHPW6AOB4YGSUmZCpj5o+8rsUADikdRXr9OTKh2VX/UUpkbBWJeWpuG+x32XtpyuD825jTD9r7c6WVow9Lce3SxrU4XkDW44dxFp7n6T7JKm4uPiogzcAHM9aV5yN8YJzVWOVrLUyxvhcGQB4rLV6d8e7evTjR7V452I9vrNMRY0Nauo/UQljvuZ3eQfpyuD8vKSvSbqj5edzHY7faIx5StI0SdUdWjoAAFHS2uMsebsHhmxIdc11So1P9bEqAPB2Mn1h49+0+sP7NX3nJyoddKJumnSThgdypKQ+SiiYKfXA/5EfleBsjHlS0mmScowxpZL+XV5g/pMx5uuStki6qOXpL0k6T9IGSfWSropGDQCA/bWtOMsoM9HbdruyqZLgDMAX1lp9VPaRXlr5qFJXP6fz91XrsuaQgvEpOnfGnYodPN3vEj9VVIKztfbSQzx0RifPtZJuiMbnAgA+XWuPs+St8gxKG/QprwCA6NnbsFcvb35ZC9cvVH3ZWr1YukOxkur7TZCmXaf40V+U4pP9LvOIdMvFgQCA7tdxHF1rcK5srPSzJADHiZpgjV4reVWfrHhchVuWSMYq/sRTdNmpP1Nkz1ZpxFwl557od5lHjeAMAI7qOI6uNTgzkg5AV2kMNWpR6SIt+fhJ5W94Q+fU7NO8UFjNgTjVj5qry+c+4j2x9+XlNgRnAHCcMUZZiVmSCM4AoqsmWKO3t7+tlWv/qucqVmhfqF4/qm7UlyurVDeoWHby1Yobdb4yEtL8LjUqCM4A4Ki2Vg0ZpcWnKcbEEJwBHLNddbv0xtbXtX7tM+q39X2dXlunc5ublT/5SxpRfI2mJA9SjIlRWrp7G0MTnAHAUR1bNWJMjDLiM1TRWOFjRQB6o2A4qA/3fKh3t7+rt3e8rfo9a/TArt26NBRWREZ1/ScoMv4SfW3shVJqrt/ldimCMwA4zsibhTowbaC27dvmczUAejprrbbs26J3dryjNRv/obQti3VS7T5lJiSqz6izdP7U7ypz1SvSiLmKGfk5paXk+F1ytyE4A4CjOk7VkKTCjEL9a8e//CwJQA9VG6zVe7ve0zvb39G7O97V57at1rl19bqsuVmS1JDWV9MnXq2rZ33fe8H4+T5W6x+CMwA4quMGKJIXnJ/f+Lxqg7VsggIc5xpDjVpRtkLLt7+r6vUvK2fXahUGm/TigCGa1n+6zqsLaEB6SBp5vnTiHCVlD+uRO/l1N4IzADiqY4+z5AVnSSrZV6KxOWP9KAmATxpDjfqo7CMt2bVES3cvVWbJYn21slJXNTUpTlLYBFTfd7TevuB5xSX38bvcHovgDACO69iqIUmbqzcTnAHHNYQatHzPcn1U+rZqNr6ujD1rNKGhXm9mZyuuX5Fm9z9Fw8MrFBl2ljT0DAUGT1dafIrfZfd4BGcAcFTHcXSSNChtkGJNrDZVb/KzLABdYG/DXn20e7lW7lmmD/auUM3OD3XH7t2aH2xWQFJERvV9CvTYab9R8vCz/C631yI4A4CjDmzViIuJ06D0QdpcvdmnigBEQ0OoQesq1umTnUtVtfl1xe1cqcE1ZRrXFNT69HRFhp+ss8ZcrvzIywoXnq5AwUzFDJyi1MR0v0vv9QjOAOCqltxsOlzQU5heSHAGepH65nqtrVirtXs+UuXWt7WzcqP+FqmUjYT17pZSpbb8zVJ1ap7CQyfo60VX6LrRn/dePP0HPlbuJoIzADjqwKkaktfnvGj7IjVHmhUXE+dXaQA6UROs0dqKtVpdvlqry1erz8b/U2HFNo1qCurCYFDxkrak9lHfM27R6OzRCu1cJ5s9QmbARGUkZfld/nGB4AwAjjqwVUOSRvQZoVAkpLXlazUud5wPVQGw1mpP/R6tr1qvkt3LVLt1seJ2r1b+vt3KCYf12375ykvO02/r6jW6IaSG3BFqHjRN8UNmasiAyboxc5D3RoNn+/tFjkMEZwBw1IEboEjStH7TJEmLdy4mOAPdoLqpWhurNmrr7mWq3r5E4T1r9FhsUOWhWt1YWaVrq/a1Pbc2MUNNuWP0xgWPKidtgNRUI8WnKoH5yT0GwRkAHNexVaNPYh+N6jNK7+54V9eMv8bHqgB31DfXa1vNNm2p2qiKXcu1trlK6+u2K2/nKn11z04VNjdrUiTS9vymGV9T5uCTNLGuVg015UoaNFXqO16pKTnab2uihLRu/y44PIIzADiqsx5nSTp5wMl6ZNUj7CAIHIXmcLNKa0u1dd9WlVRv1paaraotW63xW5crr75Khc0hndbcrHhJPxh8opIGTNDE/Cnq3/ChGnOGaV/fiUrrP1km90RdnzFIign4/ZXwGRCcAcBRnfU4S9KM/jP0wMoHtGTXEp0++PRurgrouSI2ot11u1Wyr0Rb923VtsoNitu+VPEVJcqqLVNBMKiRzSG9npWh13MG6pTYLF1avks1KTkK9h2imrxRSu03Sb868Rwpvb/fXwddgOAMAK7qZBydJBXlFikpNklvbHuD4IzjTlO4Sdtrt6u0plSlNaXaUbVZwb1rFFe+Wan7dmpQU4OWJiZqYXqq+pp4vbppg/e62ATVZw6UzRmuW4uu0E9Hfk6KhKUvhpUZG+/zt0J3ITgDgKMO1aoRF4jTeYXn6cVNL+qW4luUkZDhR3lAl7DWqrKpUttqtrXd9lZsVLB8vUz1NqXWlasyENALqSmStVq8dbtSO/Qf1yX3UdHQL+r6WQuUl5wnbXlHyh6mhNT8gy/SiwnQcnGcITgDgKNap2qokwvyLx15qRauX6i/rv+rrhp7VfcWBhyjcCSs3fW724LxrsoNqi9bq3BlieL37VBDJKin070L657YsUsTmoL7vX5Pv3G6+Nx7NDBtoFJWPiMlZUk5w6XsYUqJT1FKxycXnNJ9Xww9HsEZABx1qBVnyZvnPKXvFD219ildPvpyNkNBj1LfXK899Xu0p36Pdtfv1t6aUjXuXa9Q5SbFVG1TqKFCD2R4F7betbtMF9Y37Pf66rR8zZx9jwalDdKgdf+QQs1S1hApc7CUOUR5SVnKa109njq/u78eejGCMwA4rrPgLElXjblK3/znN/XE6idYdUa3sNZqX3CfdtXt0u763e3BuHanmitLFFO9TYk1e9SnqVb/nZkha4xuLa/UVftq9nufYCBe/efcocHpBTpx2zJFGmsVk1UgZRVImYOVkZKrWa3BeNr13f494S6CMwA4qm3F+RCbJ8wcOFOzB83WvR/dqzkFc9Q/lSkA+OzCkbD2NuxtC8Ott721O9VUtVVm3w4F6vYoO9io/FBY92emqzoQ0BXVNfpRRaU6dgpHTIwKz/yVsrJPVMHutWqu3Ka4PkOlTG/VOD6tn74cE+M9uWVTH6A7EJwBwFFtPc6H8YNpP9Dnn/28blt0mx6c86ASAgndUBl6m6Zwk/bUtQfi1nC8p36P9lWXKr2yRIl1FcoNNSs/HFJ+KKw/Z2Voa2KKLm+Uvrt9437vFw7Ea9o5P1PmkJOVs+cTBTa/1dJK4QXjmPQBmhNoiSj9pvrwjYHOEZwBwEXBoJLe+pekQ7dqSFLflL76+ck/1y3/d4t+/M6P9atTfqUAUwKOC82RZlU2VqqisULlDeX7/2wsV0XDXpU1lKuhZocmVu5Ufiis/HBY+aGQhobD+mNOX23KG6pZwYi+W7Ku7X1D8SkKpfbTlNPvVNrQs2QqNkklb0npA7zZxun9FUjM1MjWvwkp7CsVnurTnwJwdAjOAOCir39dg594QgN+c+Jhg7MknV1wtm6quUm///D3CkfC+tXMXyk+wFza3qYx1KjKxkpVNlW2BeKOv7fer26okK3bo+rmWlUEAkqNRHRV1T7lh8Ma0RqOwxH9ecCJWlpQrNFJA/XN9Q9LkpoTMxRJG6BAxiD958k3eYG3oVKavLwlGPdTbEKaYiUlthaWPdS7AQ4gOAOAi556qv3+4XOzJOkb476hWBOruz64Sztqd+iOU+/QkPQhXVcfPlVzpFkVDRXa27hX5Q3lbavBnYfhciU11apPOKI+4bCyIxHtDgS0JClRcYrRf5VVKCdilR0KKaO5UQFrtWzYTH1S/FXlxSbrtP+5QuGUXJk+/RWTMVAmrb+uHHW+riycKYWbpck3SWn9FBeXeHChSVnSUDbSwfGB4AwALho3TrUZSdqeu+9TV5xbXTn2Sg1IG6Db371dFz5/oa4ae5WuHHOlkuOSu7jY40coEmpridjbsFfljeUH3S9vKFd5Y7mqmqraXjeqKah+oZByw2H1ixiNUUCVyVlaVDhZBRkF+vFbjyg5WL/fZ9WOOFfhC/6o9Ph0mYfPleJTpJQ8Kb2flNZPEwdM1sQBk7wn/3ivYg/VohOIk/oUdtUfCdCrEJwBwDXBoLRqlRqvv0LSe0ccnCXprCFnaVzOON219C7d+9G9+p81/6Mvn/hlnXfCeRqeOfyQEzqOV9ZaNYQaDloB3tuwt9NAXNVUJdmIMiIRpUasSuO8/wx/sSGk4rBRXxuj3HBEfZqb1JiSrffPuFXZidk66bnvKam85QI7EyOl5Eq5U/WVM//LOxZfKMXEeME4JVdKyVVqap6UkO49fvXLh/8i9LUDR8QcyVXXPUFxcbFdunSp32UAQM8XiUgPPKDaxx/U3C/v0+3z7tVpg0476rdZUbZCD658UG+WvqmIjaggvUCzB89WUW6RxuWOU05STvRr70JN4ab9en8rGitU3VStYCSo5nCzmiMdbgf8HoqE2o4FI0HVBmvbwnJTuEmpkYhyQ2H1CYeVFYkoKxxWhgL6e//hyknK0eW7tmpixQ6lhhqV1FSvGBtWMCVXu697Q9mJ2Ur+09ekjW9IqXlSar6U1lfKHyPN/pFXfOlSL9ym9ZOSc6QA615AVzHGfGCtLe7sMf7NAwDXxMRIqalKfft9pc4d/pnfZnzueP1+9u+1t2GvXt/6ul4peUWPfvyoHrIPSZJyknLUP7W/+qf0V25yrtLi0pQSl6LU+FSlxKUoPiZesTGxbbe4mDgFTGC/YzHGm8VrWv6vdXHcyMgY07Za3vp7Rw2hBtU316uuuU51zXXaF9zXFogrGisOukiuPrR/K0Nn4mPilRgTq2xrlBMx2pmUpkBsvMbX1+mk2mplhkLKCDcrMxxSRiik/511vdKTsnXayhc1dMs/93svawK66YaXvH8eb/+ntO19KTnbWxFOzVN8Wj8NShvkPfnLj0pxSdKhVvQHdvrfcADdjOAMAK75y1+kP/9ZkhQTOfw4uiORk5Sji0ZcpItGXKTGUKPWVqzVirIVWl+1Xjtrd2p1+WqVlZapIdTw6W/WDeJi4pSVmKU+iX2UlZClQWmD1D8mUScEg8q2RlmRiDLCIaU0B6WTblBcxiDFr3legTd/LVO/V2ps7y3Wd1Z62zS/dZf0z//wLoRLzpEycqTkbN08/jopMV3KmiCN+0pLMM6RknNkkrO90CxJp9x8+KLj6SMHegOCMwC45je/kd5/v+UXG9W+5MTYRBXlFakor+igx0KRkOpD9aoL1qm2uVbBSFChSOjgm22/b2XbNmpp3enQWrvf8dbHrLUykbDimmqU0FSj9FBQkezhissYqMyaPeq36jklNtUq0FAls6dcqiuRLnpMKjhZWvkXaeHX9y84NlGadKUUl+z1BvcvkpL6tIfflBwvKEvS9BukGTcdukViYDGrwsBxgOAMAC555BEvNE+eLH3wgUw3XsYSGxOr9Ob83ioAACAASURBVPh0pcenH9kLmhska73V1mCdtHmR1FDlzQWu3yvV7ZXGfkk6YZa0a6X0yNz9V4Ml6YL7pb5TpOC/pE9ebVvtVc5wafBJ7cG3YKb01We8UNwajjuu8g6d7d0OpbMxbACOOwRnAHDFK69IN9wgnXaadMMNCl99law59laNIxKJSKXvS/UVXrhtqPJ+9p8kjThHaqqRnrhw/8dCjdLpP5Jm3eqF5ScvaX8/E+OF29ZV3JQ8adyFXihuXQ1OzpHyRnuPD54ufX/zoetLy/duAHAMCM4A0JsFg1J8vLRjhzR3rlRQIP3v/0r9+mnVrOHa/NLl0WvV+OQVqXyDVLFJqtwsVW3zVmnPvcN7/JHPSZHQ/q+ZfoMXnGOTvHnAOcOlxEwpKdP7WTDTe15qvjT/jfbjiZnt/cGSF3o/d1d0vgcAfEYEZwDoSRobpepqqarK+9nQIM2a5T22cKH04Yfe8V27pGXLpMxM6YMPpP79pZdflmbMkBK9toLWHuGjWnGOhKXdq7wJELtXeXOAz/6Z99hL35OqtnjH+hRKuSd6QVjyQu5Xn5USUtuDcUJ6+3zgQKx05QuH/txAnNS6GQcA9FAEZwDoLps2SStXeqvDO3ZIO3dKZWXSs896Y8jmz5ceeGD/16SlSfv2efdbp2VkZEjZ2dLEidL06e3Pnd2hR3fJEg35+Q/Vf2bwyOt7+QfSsiekppbPS8yUhp/V/vhlf/baI5L7dD42rXDmkX8WAPRCBGcAOFKRiBdiKyu9FeHWn2ed5QXcN9/0gm1lpVRR4a0K79ghrVnjBd0HH5R++UvvvWJipPx8qV8/b1U5OVn6/OelwkJvFTkjw/uZmdn++Y8+6rVhHEnrxY4dynr+VaVPHtr5ivO2JdLKP3sX5F27SIqNl9L7exfjDTlZGjxNyhi0/2fljjimPz4A6O0IzgCOX3V13opvdrYXfDdu9Nohysq825493s8HHpAmTPCC69VXH/w+H30kjR8vrV0rPf20lJXlBd7Bg6Vp07zALXkryhdc4LVV5OVJgQO2OT7/fO92KPHxR/7dWgJvjJX2y81ln3gtF5v/zxvHNuxMqaHC26luxreO/P0B4DhEcAbgjlBIKi09OPiefbYXfFeskL7xjfbj9S07yf31r9K8edL69dL3v+/1COfmere8vPZV12nTpP/3/7xQnJXVHpBPPNF7/LrrvNuhFBR4t+7QWrPt0ONcuUW6f7bXTzznl9Kkr3k9yQCAI0JwBtCzBIMHt0IMGiSNGSPV1Ej/8R/e8Y7P+da3pKuu8laMR448+D1TUrzgnJzsrS6PHOkF4tZwPKnlorTZs73PSEnpvB1i9Gjv1hu01G9shw1QMgdLM2+Wxl0kZQ7ysTgA6J0IzgC6RiQilZdLu3dLcXHSiJb+2F/9ylvtbe0DrqyUzjlH+rd/k5qbpYSEg9/rlluk3/7Wu/9f/7X/iu+AAVJqy6rpoEFeH3HHUJyX5wVhSRo2TPr73w9dc3z80bVD9GRJSQr2y1MoYLwV51CTFJsgzbzF78oAoNciOAM4Otu2eRe97d3rtTzs3i3l5LT3/s6Z47VElJVJ4bB3bN48rx1C8lodGhvbg29WVtv4NMXFecE6La29DSIrSxoyxHs8La29vaIzycmd9yAfj844Qx99+ILW/eNqJe7bJT0+Wpr3R2n4mX5XBgC9FsEZOB5Z610Yt3evd2tokGa2jBJ74AFp6VLveHm597NfP29XOkm68EJvS+eOTjmlPbCOGuVdFJef76325udLw4e3P3fnTin2MKeeBQui9z0hScrYtsTbwjp7qN+lAECvRnAGXBCJtO+ytm6d9PHH7RfAlZV5fbsPP+w9fv313v2mpvbX5+Z6z5e8TTQWLfJWkXNyvNDbevGbJP38595rs7O91+Xnt7dKSNLvfnf4Wg8XmhE9H36o4d+7SSec2ai4PlWSjDdeDgDwmfFfMKAnstbr/e04HeLcc6WkJOlvf/Nm+bYe37PHWxWurvZ6ee+7z2uHaJWR4a38hkJeaD3lFCk9vT0Y5+R4AbjVn/98+DnBZ5116MfQc5SXK/ONxUqfUai4hiopKcvbvQ8A8JlxFgW6g7XeqnAg4LUqLF588Mi0X/7S2/ziwQe9kWah0P7v8ckn3upvaanXSpGXJ51wgjciLS+vfVbwt74lffWr7RfHHXix22WXebdDOZLNNdDztc5xjkixjdVSSo7PBQFA70dwBj6LzlaEJ0/2LmL7+GOvnaFjKN67V3rxRW+e8OLF0pe+1P5erSvCFRVecB4/XrrttvaJEK0BePBg7/nXX+/dDqW75gSjZ2sdRyeptv8E9Rl8sr/1AIADCM5AK2u9sNHQ4G2dvHu3d9u1y/v5la9Ic+dKK1d6c38PXBF+8EHvArmmJumDD7ywe8IJ0vTp+wff00+Xli3zQnFOzsErwlOmeDfgWLTNcZb2jrtAg/Mn+VwQAPR+BGe4rbm5PfimpXmzhINBb0X3wGB8443ST3/qXUh33nnt75Ga6l0Ad2bLGK8BAzpfET7hBO/xSZO8topDaR3BBnSl1FTVDxuixnijmEiz39UAgBOMtdbvGo5IcXGxXbp0qd9loKeIRLz2h+3bpR07vFturvTFL3qPn3yytGGD1ybR+v/jV10lPfSQ93turhde8/Olvn29n+ee660oRyLSe+95x/Lz2zfPAHqZd7e/q+tfvUbLt2yXmbVAOu37fpcEAD2eMeYDa21xZ4+x4oyeaf16qaSkPRRv3+61Ndx+u/f42LHSmjX7v+ass9qD88iR3hbN/ft7M4j79m3fuc4YL3QfSkyMdNJJ0f5GQLezssqIRGRsRErM8LscAOj1CM7oPo2N7TvEvfGG1wfcGox37PC2Wn71Ve/xa6/1ntMqM7N9gw5J+s53vJaL/v291on+/b1w3OrBB7v++wA92YoVGvuN6zXlrDopTkzVAIAoIDgjOiIRr0+4b19vRfeFF6SXXpK2bPFupaVev3Ftrff4ww9Ljz/ubZHcGnxbt1WWvNFszc3tK8bJyft/3jXXdO/3A3qbmhplLFmpvif1l7JEcAaAKCA448g0N3vhd8AAbwrE669LTzzRHoy3bfNWgCsqvN7hxYulp5/2wvDw4dLs2V4IDoe9TTjuuku6+25vI47O5gZPn9793xFwScu/V+nhlvneyQRnADhWBGd46uqkrVulgQO96RPvvSf94Q/twXjHDm9VeflyacIEaeNGb2vmIUOk4mJvLvGQIe3bKf/sZ9IvfnHoz+u4Ux2A6GsJzntjAiorulS56f19LggAej+C8/HCWm/CREKCt+HGunVesN2wwQvBe/Z4z3vhBelzn5OqqrxV4yFDpDPO8H4OGeKtGkvS/Pne7VBiYrr+OwE4tJZ/B0vjYrX7lBuVm9zH54IAoPcjOLskEvHaJRITpfJy6de/9kJx662mRrr33vbtnN98Uxo6VDr/fG8G8ZAh0sSJ3nvNmSNt2uTr1wFwDFJTtW/CSNn4BgWC9X5XAwBOYI5zbxUOS3/+s7dyvHatd/vkE+nmm73tnquqvM05CgulYcO8gDx0qDeybfRov6sH0A3e3PamGp+6VKfH5SjhOyv9LgcAegXmOPdWZWXSqlXevOK1a72QPGGCt5IcE+NNlqit9cLxiBHeVs6zZnmvzcz0to4OBPz9DgB8Y61VVjisUHqmEvwuBgAc0OXB2RhTIqlGUlhSyFpbbIzpI+lpSQWSSiRdZK2t7OpaeqyaGunjj72QHAxK3/ymd/zUU73ALHnbPo8c2b6LnTHShx96F/O1zkY+EKEZOH6tXaviL81X9fQmhWdl+l0NADihu1acT7fWdtyqbYGkf1pr7zDGLGj53f29YCMRb2xb67ziH//YG+lWUtL+nGHD2oPzb37jXcw3erR3Ud6BY9uGDeuWsgH0Qg0NSlu9UaEJKYokpPldDQA4wa9WjS9IOq3l/qOS3pSLwbmkRFq0yFsZ/vBDb5RbQ4PXXtE63eKkk7zpFGPHSuPG7b8JyNy5vpUOoJdrmaphrKQYuvIAIBq642xqJb1ijLGS/mitvU9SvrV2Z8vjuyTld/ZCY8w1kq6RpMGDB3dDqZ9RMCitXNkekG+/XcrP9y7eu+02b9e7oiLpiiukSZO8lWdJ+t73fC0bgMNa/obqteQkTRs6WzRrAMCx647gfIq1drsxJk/Sq8aYtR0ftNballB9kJaQfZ/kTdXo+lKPUCTireYsWeJNsfjgA6mx0XssI8MLyPn50uWXe6vGJ55IvzGA7tUSnN9NStK4gpN8LgYA3NDlu1RYa7e3/Nwj6RlJUyXtNsb0k6SWn3u6uo7PLBiU3n9f+v3vpUsu8VopHn7Yeyw93RsLd/313vbSGzdKlZVe+4Uk9esnjRpFaAbQ/dLSVD6jSIEko0D98XvtNQBEU5euOBtjUiTFWGtrWu6fLek/JD0v6WuS7mj5+VxX1vGZ7dvnrRy3riYPHOiF4ta2kREjvN31AKCnKSjQB/9zp/7tsYtl3n9AGjjd74oAoNfr6laNfEnPGO+vDGMl/a+19mVjzBJJfzLGfF3SFkkXdXEdn016uvTDH3oB+aSTvOAMAL2ElVWMJHvgRB4AwGfSpcHZWrtJ0oROjpdLOqMrPztqfvQjvysAgKO3aZNmnn6ldHJQwfFd3pUHAMcFzqYA4KJQSMlbd8k0WclwqgeAaOBsCgAuamnPMFYEZwCIEs6mAOCiluD8YmqK6ob3js44AOjpCM4A4KKWnQM/SExQU/8in4sBADcQnAHARSkp2nXmNMWnGsXWlvldDQA4geAMAC7Kz9fSP/5IP0zfp/TlT/ldDQA4geAMAI6y1kqSDBcHAkBUcDYFABft2KGzp35VgRXNTNUAgCjhbAoALopElFBRrZhmK071ABAdnE0BwEUdt9mOYcttAIgGgjMAuKglOD+TmqKGYcxxBoBoIDgDgIta5jivjE9QKG+Uz8UAgBsIzgDgoqQkbT1/phIzpdiaXX5XAwBOIDgDgIsyMrTs1zfqtqxaJa950e9qAMAJBGcAcJS1kZZ7XBwIANFAcAYAF1VU6PyJV0hLgm39zgCAY8PZFAAcFWhqlsJWrDgDQHQQnAHARR1Xmdk5EACigrMpALiodY5zSoqCzHEGgKggOAOAi1pWnDfExSvSp8DfWgDAEQRnAHBRfLw2XnyGEvKMAtU7/K4GAJxAcAYAFyUkaOWCS/TtPnVK2PSG39UAgBMIzgDgImtlQiHJWi4OBIAo4WwKAC5qatIXplwvvR0U4+gAIDoIzgDgItMhLLPiDABRwdkUAFzUGpytZNg5EACigrMpALioJSw/k5qi5hNO97kYAHADwRkAXNSy4rwjNlZK6+dzMQDgBoIzALgoJkZrvjZHsQNiFFNd6nc1AOAEgjMAuMgYrfvmHF2bXa/Ybe/7XQ0AOIHgDACOittXJzVbKSbgdykA4ASCMwA46ryzfyK91bT/aDoAwGdGcAYAV5nWH5zqASAaOJsCgKOsMd4cZzZAAYCo4GwKAI6yxujZ1BSFC2f6XQoAOIHgDACuMkblMQGZ5Gy/KwEAJxCcAcBRK685V4ETAjLV2/wuBQCcQHAGAEeVXFysK7MbFLN7td+lAIATCM4A4KjEPVVSo2UcHQBECcEZABx15lfvkf6viakaABAlnE0BwFneSrMx7BwIANFAcAYAR1kjyYpWDQCIEoIzADgqEhPQ86kp0uCT/C4FAJxAcAYAVxmjGmOkhFS/KwEAJxCcAcBRy66bo5iRsTLV2/0uBQCcQHAGAEftnjNMl2Y3ylRs8rsUAHACwRkAHJWyba9UE5GYqgEAUUFwBgBHzbz5aW+OcwynegCIBs6mAOCo1nF0hlM9AEQFZ1MAcFXL+GYTQ6sGAEQDwRkAHBUOxOmllGSp/0S/SwEAJ8T6XQAAoGtYY9QoI8Um+F0KADiBFWcAcNSH182WJsRJ+3b4XQoAOIHgDACOqprWVxdkN0o1u/wuBQCcQHAGAEdlbNwjVUQkw6keAKKBsykAOGrKL/4h/V8TwRkAooSzKQA4yrbeITgDQFRwNgUAVxnjpWdj/K4EAJxAcAYARzXHxeu15BQpd6TfpQCAEwjOAOAoa6SQMRI7BwJAVBCcAcBRH157qiLT4qXaMr9LAQAnEJwBwFGNY9J0XmaD1FjtdykA4ASCMwA4KnPNbml3mIsDASBKfAvOxphzjDHrjDEbjDEL/KoDAFw1/r8XM8cZAKLIl7OpMSYg6b8knStptKRLjTGj/agFAJxHcAaAqPDrbDpV0gZr7SZrbVDSU5K+4FMtAOAka8QcZwCIIr+C8wBJ2zr8XtpybD/GmGuMMUuNMUvLyrgqHACORmN8khalpEoZg/wuBQCc0KP//s5ae5+1tthaW5ybm+t3OQDQq1hjvGVnVpwBICr8Cs7bJXVcAhnYcgwAECXLr5mu8KwExtEBQJT4FZyXSBpujCk0xsRLukTS8z7VAgBOskNidXp6vdTc6HcpAOCEWD8+1FobMsbcKOkfkgKSHrLWfuxHLQDgqj6rdkulIaZqAECU+BKcJcla+5Kkl/z6fABw3YgnVkgxQXqcASBKWIYAAEdZ0zKPjhVnAIgKzqYA4CojhZmqAQBRQ3AGAEfVJKbqw+RUKSnL71IAwAkEZwBwlZGM9bsIAHAHwRkAHLXq60UKnpPEODoAiBKCMwA4KjY3ohkpdZIN+10KADjBt3F0AICulb1ij7SdOc4AEC0EZwBwVOFzG6WGJoIzAEQJZ1MAcFXLGGeCMwBEB2dTAHBUJMaoWUZeggYAHCuCMwA4qioxXeuSUqQAXXkAEA0EZwBwlGWOMwBEFcsQAOCoT64YobztK/wuAwCcwYozADgqOS2oiYm1fpcBAM5gxRkAHNXnozJvjjMAICoIzgDgqAH/3C5T3uR3GQDgDFo1AMBVrXOcAQBRQXAGAEeFY2LUbJjhDADRQnAGAEftTc5UaXyS32UAgDMIzgDgKMtiMwBEFRcHAoCjtl46WH13NfpdBgA4gxVnAHBUWkKjhsXX+10GADiDFWcAcFT2snLFbmccHQBEC8EZAByV+68yxW+jVQMAooVWDQBwFXOcASCqCM4A4KjmmFjmOANAFBGcAcBRO9OyVB6X4HcZAOAMgjMAOMzQqgEAUcPFgQDgqIov5cjsbVB/vwsBAEew4gwAjko3DRoQaPC7DABwBivOAOCorGWVSiolOANAtBCcAcBRWSuqlLyOnQMBIFpo1QAAh3FtIABED8EZABzVEBevMKd5AIgazqgA4KjS1D6qDcT5XQYAOIPgDACOssbIWJo1ACBauDgQABwV/FyySmb1Ux+/CwEARxCcAcBRGZEGZdsmv8sAAGfQqgEAjsr8sFoZr9T4XQYAOIMVZwBwVNr6WqUur/O7DABwBivOAOAqY/yuAACcQnAGAEfVxCcoLMIzAEQLwRkAHFWSlq2gCfhdBgA4g+AMAI6KxBjZGFacASBauDgQAByVOlv65NTBmux3IQDgCIIzADgqrblJseGg32UAgDNo1QAAR6Uv26fcv1b4XQYAOIPgDACOStraqIwPav0uAwCcQXAGAEcZSbJ+VwEA7iA4A4CjyhOTyc0AEEUEZwBw1PqMHEXYAAUAoobgDACOaooPqD4l3u8yAMAZjKMDAEedOLFGG4oLmOMMAFFCcAYARyWGQ4rE8BeLABAtnFEBwFGpH9Vq0KM7/S4DAJxBcAYARyXsalbmcuY4A0C0EJwBwFWtAzUsQ+kAIBoIzgDgqB1Jqd4dgjMARAXBGQActTo9x7tDcAaAqCA4A4Cj6lNiVZGTQnAGgCghOAOAo04/Ybc2/aRAimXyKABEA8EZABwVYyMSW24DQNR0WXA2xtxujNlujFnecjuvw2M/MMZsMMasM8bM6aoaAOB4lrSiQcPuLZGCQb9LAQAndPXf3/2ntfa3HQ8YY0ZLukTSGEn9Jb1mjDnRWhvu4loA4LgSWxFS+ppaKczpFQCiwY9WjS9Iespa22St3Sxpg6SpPtQBAE4zrV0akYivdQCAK7o6ON9ojFlhjHnIGJPVcmyApG0dnlPacgwAEEUbkzO8O0zVAICoOKbgbIx5zRizqpPbFyTdK2mopCJJOyXd9Rne/xpjzFJjzNKysrJjKRUAjjsrM7K9OwRnAIiKY+pxttaeeSTPM8bcL+mFll+3SxrU4eGBLcc6e//7JN0nScXFxZz5AeAoNKTGatfgLPU1TNYAgGjoyqka/Tr8Ok/Sqpb7z0u6xBiTYIwplDRc0vtdVQcAHK8uyN+mzQsKpNRUv0sBACd05VSNXxtjiiRZSSWSrpUka+3Hxpg/SVotKSTpBiZqAED0GUmW1WYAiJouC87W2q8e5rFfSPpFV302AECKX9WosUvWSV/ZJ6Wn+10OAPR67MMKAI4K1ESUvLFOam72uxQAcAJbbgOAo4xarqlmjjMARAXBGQActSKtj3eHcXQAEBUEZwBw1DLmOANAVBGcAcBRwdQYbRuRJwUCfpcCAE4gOAOAo76eXqLNtwyVcnL8LgUAnEBwBgBHMccZAKKL4AwAjgqsadLkny6Xdu70uxQAcAJznAHAUTH1VqmlDVIw6HcpAOAEVpwBwFFtXRpM1QCAqCA4A4Cj3slsuSiQ4AwAUUFwBgBHLU3P8u4QnAEgKgjOAOAok2q0aeJAKSHB71IAwAkEZwBw1HfjN2vTTcOlAQP8LgUAnEBwBgAXtbVnMMcZAKKF4AwALrJWWtesGbe+K23c6Hc1AOAEgjMAuMhGpGYpuayROc4AECUEZwBwkY10uM9UDQCIBoIzALgoJqCX+uR69wnOABAVBGcAcFFMQB+kZXr3Cc4AEBUEZwBwUSSshJSIPpkxVEpN9bsaAHACwRkAXNRcr9titmrzN0dIBQV+VwMATiA4A4CLWtozrOE0DwDRwhkVAFxkI9KGkGZf94q0YoXf1QCAEwjOAOAiG5HCVvE1Qam52e9qAMAJBGcAcB1TNQAgKgjOAOCiuGQ9nZPv3Sc4A0BUEJwBwEVxifowpWUMHcEZAKKC4AwALgo1KSsppDVnjZb69PG7GgBwAsEZAFxUX6EFdoc2zR8tDRvmdzUA4ASCMwC4yEa8HzI+FwIA7iA4A4CLbETaHNK5ly2U3n3X72oAwAkEZwBwkpWsFGiOSOGw38UAgBMIzgDgIhtRW5cGUzUAICoIzgDgopRc/TEnz7sfifhbCwA4guAMAC6KT9FHKcnefVacASAqCM4A4KKmWg2Oa9bKL0yS+vf3uxoAcALBGQBctG+HFoR3a/OVY6QRI/yuBgCcQHAGABfZiNeiEbH0OAPoPVaulK65Rrr6ar8r6RTBGQBcZCNSaVif/9IT0quv+l0NAByatdIzz0inny6NHy89/riUmNgjr8+I9bsAAEBX6PAfnB74Hx8AUE2NlJYmGSPdeae0a5f38+tfl7Kz/a6uUwRnAHARc5wB9EShkPSPf0gPPOD9bdiWLV5IXrhQ6ttXCgT8rvCwaNUAABdlFejOvHzvPsEZgN927JAWLJAKCqS5c6V335VuuKH9GowBA3p8aJZYcQYANyWkaXVionefiwMB+GH7dqmuTjrxRKm+XrrrLunss6W77/bCc1yc3xUeNYIzALiovkKjYhu1/NIZKho61O9qABwv9u2T/vpX6YknpNdfly64QPrLX6Rhw6Q9e6SsLL8rPCa0agCAiyo2a0GwXJsvHy+NGuV3NQCOB9/+tpSXJ111lVRSIv3kJ9KvftX+eC8PzRIrzgDgJhuRIlaxtUEpGJTi4/2uCIBLduyQnn1Wevllb0U5Pl4aOlT65jelL39Zmj7dm5bhGIIzALjIRqQ9EZ1/8UPSX+dK8+b5XRGA3q6kRHrwQenFF6Vly7xjI0ZIW7d6rRg33eRred2BVg0AcJHtcEEgUzUAfBYVFdJTT0kffeT9vmOH13qRmur9/Phjae1aLzQfJ1hxBgAnWeY4Azg64bD0z39Kb7zhXdi3dKk3lefWW6UJE6Rp06SyMid6lT8rgjMAuKjvOH0/L093ajPBGUDnmpqkf/1Lqq6WPv95ryf58sulykqvR/lHP5LOPVeaMsV7fiBwXIdmieAMAG5KSNP6hJYZqcxxBtBq2TLp73/3VpXffltqbJSGD/eCc0yMt5vf0KFeOwYOQo8zALho3w6dpAZ9MP8MaexYv6sB4IeaGum116Rf/rL9f0Dfc4/0wx96M5WvvVZ67jnp/ffbXzNhAqH5MFhxBgAXla3TrcFqLbyoSJNHj/a7GgDdZflyb/LFO+94F/VFIl4LxsUXeyvJ//7v0h13SLm5flfaK7HiDAAushEpbJW0p9bb8haAWyorpVdekX7+c6/NonXVeMsW6eGHvV7kH/7Qm7NcWemFZkkaPJjQfAxYcQYAJ1mpKqLzLvuj9MRM6bLL/C4IwGfV2Cg1NHhheMMG6bzzpPXr2x8fNUoqL/fun3eeVFUlxRLxugJ/qgDgIBuJyLTu2sVUDaD3iES8i/aWLWu/rV4tffe70p13SgMGSOPGSVdfLU2dKk2eLGVktL8+Ls6/2o8DBGcAcJC14bYxzkzVAHoga6WdO72e5GXLvPB7441eP/K8ed7mI337ShMnSnPnemPhJCkpSVq40N/aj2MEZwBwkB0yQ/P75up+1bLiDPgpFJI2b5Z275ZOOcU7dvXV3jSLior2533u/7d35/FxV/X+x19nZjIzmexrl6Q7XaBUSimbLGKBUpFFQBZF0As/kU2Qq1euotf7uw8q/h7KovcqFxQEuaJAK4qAF0R2CoUCLaUt3Uu3kDZJm32bzPn9cb7JTNomHegkk07ez8fj+8jM9/udyZnDN+Hdk8/3nM/Hg/Pf/uZqkUeOTE+bpU8KziIiGcgGc9gc9P5kq+AsMvCamuLTuD30EPz5z2456rVrobMTSkvdqnsAEyfCRRe52uQjj3RTNsbHAAAAIABJREFUwOXnx9/rmGMGv/2SFAVnEZEMZOvWc2ashTe/dQ7HdK/6JSKpsXKlW5r6gw/iW1WVC8+RCKxa5c6ZNs3NeDFtmtusdSPKP/hBuj+BfEIKziIiGchUr+LGlgYWnHckx0yfnu7miBxc2tpgzRpYvbp3OH70UTet23PPwY03urrkQw+FuXNdMI5G3et//GO3ScZRcBYRyUTePM75m2pgxi43jZWIxFnrSicSg/FXv+rKJp5+Gi64wJ1nDIwb54Jxe7vbd9llbkGR8nJ3XIaNAwrOxpgLgX8HDgWOsdYuSTj2PeBKoAu4wVr7jLd/HvBzwA/8xlr7kwNpg4iI7M3aLmi2zP3qL+GeT8FVV6W7SSKDr7XVLQiycSNs2gSzZ8PRR7tSihNOcAuDdMvOhmOPdcH505+GP/zBjSZPnuzKLxLpH6LD1oGOOL8PnA/ck7jTGHMYcAkwHRgNPGeMmeId/iVwOrAVeMsY84S1duUBtkNERBLZGD3z0enmQMlUHR2webMLxRs3ujKKOXPcAiCHHgoffdT7/B/+0AXnUaPciHF37fG0aTBmDPi8BZVHjoRLLhn0jyND3wEFZ2vtKiA+yX7cucAfrbXtwEZjzDqg+xbRddbaDd7r/uidq+AsIpJCNtaV8ETBWQ5S0Shs2+ZC8caNUFwM557rrumpU2H9+t7zlH/tay44FxS4m/LGjIEJE2D8ePe1e3q3wkK4++50fCI5yA1UjXMF8EbC863ePoAte+w/tq83McZcBVwFMHbs2BQ3UUQkc3VNPYMvVoxkAesUnGXo6upyo8LdI8bGxJeHnzMHXnklfsMdwKmnuuBsjPsaibhA3L2NHu3OMwbuuWevbydyoPYbnI0xzwH7moH7FmvtX1LfpDhr7b3AvQCzZ8/Wb34RkSTZQDY7s7xf8Vo5UNKho8NN0bZ1a3xra4NbbnHHL77YrYDXlfDXkSlT4sH5tNPg+OPjo8UTJrgR5G4//emgfRSRbvsNztba0z7B+24DEq5uKr199LNfRERSxFSv4EttTbzxwy9z3Gc/m+7mSCaqqnIzUSQG4x073JRtxsDXvw6/+13v15SXx4Pzqae6G+8qK+PhOPGvy9///qB9FJFkDVSpxhPAw8aYO3A3B04G3sTdqjLZGDMBF5gvAb48QG0QERm2fNUruLq5gQWfnwWHHZbu5sjBorusxxi34t0rr7hAvG1bPBy/8IKrNf7lL2H+/Phri4pcCG5pgZwcN7XbySe7fd1b4up4mulFDkIHOh3decB/AmXAU8aYpdbaM6y1K4wxj+Ju+osC11lru7zXXA88g5uO7n5r7YoD+gQiIrIX01QNMUv+Bx/B2Co3i4AMb52dUF3tAm5ODrz/PvzP/8RD8bZtblu0yE3J9o9/wDXXuNeWl8dHhtva3L7LLnN1yJWVUFHh3jPRnDmD+vFEBoOxB8lNI7Nnz7ZLlizZ/4kiIkLnI5ex872/MvrH9XDnnfCtb6W7STIQ2tpcycT27S4U79gBn/mMm4pt2TK44Qa3b8cOqKtzr/nrX+Gss+Cpp+C889wNdRUVbqushG9+05VN1NVBQ4P7R1colN7PKTKIjDFvW2tn7+uYVg4UEclA/ur3WRMMMho0q8bBqqMD3njDBePucFxVBRde6KZaW7kS9rWc+n//twvOwaAruZgxw40Yl5fDiBHuOcC8eS54d89dvKfiYreJSA8FZxGRTNPZhmnawepQkFNAwXmoaGlxSzz7fG52iFgMbrvNTcdWXe22jz6Cr3zFLdTR0uJGj7sFg250+MQT3fMxY1yN8ahRbhs50gXjsjJ3/NBD4cUX+26P3z9gH1UkUyk4i4hkmqwwDTct5+HffppvgKajG0ibNrlR4J0741tFhav/BTel2rp1bn9Li9t36aWuttjng5/8xAXYESPcdsQRMHGiO6+gAJ591oXlUaNcbXLigmN5eZp5QmSQKTiLiGQin59Wv/cneI04719XF+zaBbW1rnzhiCPc/vvvh/fe6x2MJ0yAP/3JHT/rLFixxz3uc+fGg/OECS5Il5XFt8RZTmpr3UjyvhgDp5+e2s8pIgdEwVlEJNO89FNCbXV0ZPl45fbrOenMs9PdosFlbXxkdvVqt9XUxLfOTnfDJMDVV8Njj7nQ3P0PjHHj3EgyuDmJFy2Kh97KSpg2Lf69br/dhe7EYJw4u8Svf91/W/sKzSIyJCk4i4hkkq5OWPYw/sIxdPkNmz93XO+gdzCKxWD37t6jvmedBVlZbuR34cLex+rqoLHRlUDceWfvpZdDIVf6cMcdLlzPnOnOKymB0lL3NXHqvqee6r8W+IwzBu5zi8iQo+AsIpJJ/vEfULeB9pP/BfPujxm5aAWEN7qSgXSz1gXaXbvcDA/Z2bB+vVtQY9eu+FZX5wLv6NHwn/8JN93Ue1lmgC1b3Ojvhg1u5onu0eBZs9zjzk4XeL/9bbjyyvhocCTSu0746qv7b7NuoBORBArOIiKZINoOr94Ji34Bs6+kY9o8Akvmc+rXb4P5uQN3E1l9PSxZEp8yrbbWjQ5fey186lNuEY1vfMPt2707HoBffhlOOglef90tzQwQCLgb4IqKXIAePdoF4Ztv7l0KUVbmgjfAd77jtr5Mnjwwn1tEhiUFZxGRg1XTTmj6CEbOgI5mWHwPHH4BnPFjbLQZ2z2weiA3B3Z2wvLlbjGNNWvcCPH69XDLLXD++e7GuNNOi5+fleWC7znnuOBcUgLHHOP2FRbGg/Ehh7jzzzkHPvzQ7cvN7T0aDHDCCW4TERkCFJxFRIaKuo1QvwXaGqCt3m2xKJxwgzv+8s9gw4vQ3gBNO6CxCoonwQ3vQKQYrlsMuW4k1kabsN0hNNngHIu5IPzyyy7YnnGGW3TjqKPc8UDALbl8yCGuzALcYhovvOBGh0eOdFOkJYbfmTPh4Yf7/p75+W4TETkIKDiLiAyUjhY3Ilw43s3Zu+FFWP88NH7kQm9jNbTUwHfWueOv3gHv/K73e2RF4sG5o9kF6dyRUD4dRkyH0TPj53qhuVvPiHN/8zjHYvDzn7vw++qrrkQC4PrrXXAeO9bNOjFjBkya5MJzorw8OOWUj9kxIiIHJwVnEZH+RDvA+MAfcIG36j1o2+1Gg1t3u8cnfAtyy2D5Anjtrvix9gb3Ht9eA3kjYPMb8MbdLvjmjYSyKZB7MnR1gC8Mx14DMy6CcAGE893XUMJo7Gk/+lhN7xln3nPEubnZjSp/7nMusD/wgJu7+IIL4OSTXe3xuHHuXGPgi1/8JD0nIpJxFJxFZHjp6oTmGjfS2+xt446HgkrY9ja8coe3f6c7p60erngWxh4L61+AP+8xC0NWBGZ+2QXnrAjkV7rR4HCBC8t5oyAYceee9G34zM171/F2G3HYvvd/AtZarIGX7/1XTv7spfEDCxbANddAa6u7kS8Ugtdec/XFIiLSLwVnETm4RdvdiGpW2IXcNc8mhOKd0FILR18Jk+bA5sVw/9y93+PCB1xwjnZA3QaIlMCoIyCnFHLKIN+b1/eQ0+DK5yC70BsVLoRAwgIW0850W1/8WSn96P2xuEVAqj8zK34j3vz58IMfuJv17rrLhWZQaBYRSZKCs4ikl7WubtefBbEu2Px6vASiu+RhzDEw+XRo3QUPX+wdr3fnRNvg9P+AE26Eljr40/9x72t8LgDnlLmb7QCKxsMp34ccb3+k1IXjgjHu+Ljj4drX+25rbpnbDgLWK88Y/dxbEJ3uZsf44Q/hkktcaUZ3aBYRkaQpOItIanW2uhkfmnd6X3e4mt6p89zxP14KTdUu/LbucuH3yMvg7Lvc8Qc+v8cbGheKJ58OgTAEQq42OFzojfwWwtjj3akFlXDdmy4Uhwtd/S4uRLZFW2nNyqJ19ldo7WylNdpKW1cbre3VtGzfSDQWJWZjdMW6iNoosViMLttFwBcg6A8S9AUJ+UME/UHCgTCRQIRwIEx2IJtwIEzYHyYcCGOtJWZj7j1sjGgsSpftoivWRZft6nke9AUpzi4m5B+YAGu9Cufjb7gdbs5yq+GVlMCvfqXQLCLyCSk4i0jfrHUzOXQ0u3pdgNX/C7s2utHd5p1uKxoPZ8x3x//raDelWqIp8+LBubMVgrku5HaH38pj3DGfH776JF3BHJoDWTT5AzQaaIq20LTlJRo7G2k+9ivua2czrdFWWqONtG5YQOuah7znrbREW3oet0ZbaYu29QTJoSYnK4exeWOZUTqDGWUzOHrk0VTkVqTuGxjj/jtefz1cfrmmfhMROQAKziLDTdNOaNzuan9b6tzXWBccf607/swtbtq0llq3dXVA+WHxEoZX74AtiwHjSiFyy3tPg3bK94jGOmkJ5dCUFaExGKYu4Kdh0zM0djTSdOQ5NHQ00NjRSGNnI00dNTRufISG1b92xzubaO5s3u/HCJgA2YFst2Vl9zzODeZSFimLH/O27tHhSCCyz2NZviz8Pj8BE8BnfPh9fvzGTzQWpb2rnY5YB51dnbRGW2nvau8J5D0j19FW2qPtYIi/h/HjMz4CvkDP5jd+/D4/7dF26trqqG2rZf3u9Ty98WkeXfMoANNLpnP+5PM5a+JZRLIin+g/c3ephu0OzqDQLCJygIw9kBWlBtHs2bPtkiVL0t0MkaGlo9kb9a2FilludHHt32HjS/EZI1pq3HnXv+Ves/DrsPzR3u+TXQQ3b3KPn5+PrX6frnABrcEILcFsGsJ5bBhzJPXt9UTrt1Df1UaN7aSxs4WGTi8EdzTS1NFEY0cjbV1t/TbbZ3zkBfPIy8ojL5hHfjCf3GCu2+ftzw3mkpuV2/M1L5jX63nIH8L0NTvFQShmY2ys38jLW1/mqQ1PsXrXagpDhVxx+BVceuilBP3B/b9JgqqmKuYunMvSq1bj7+h0czCvXdv3jB4iIgKAMeZta+3sfR3TiLPIUNNS51aQ654OrXt2iM98183k8NZv4NW73P5oa/x1N3/oyh42vQJv/hpyyohFiolmF9JeUMHG6ndoiLbgHzODroISan2GGhOj2kb5KNbG7r9dTn17vdui9UQbor3bteH3PQ8DvgD5wXwXeL1QOzIyMh58E7bEc7qfZweyMyr0poLP+JhUOIlJhZP42vSvsXTnUu557x7uePsOFq5dyI+O/xFHjzw66ffrKU3p7uemJoVmEZEDpBFnkYHSXR/cUuNGhEsmupHd6hWw9GFXBpE4ZdqXH3Erwb11Hzz1z73eKuYPsf6SB6nJKSS07nlKNrxCY1aY+kCAOr+fGh+8k1vArmgLze311Hc00dDZSHtXe79NzA5kUxAqoDBUSEGwgPxQPgWhAgqCBT3780P5Pc8LQgXkBfMI+8MKvoNk0fZFzH9jPlsat3D1EVdz9RFX4zO+/b5uW9M25i2cxy/DX+bk2x6Bri5YvnwQWiwicnDTiLNIKnW2QvXKvUeED7/AlUtseRMWXLHXiHDnJb+ndsxsolsWMerNe+kI59MSjNCUFWJ3Th5Pv30n6wMQaqimZNw0ttt2tsXaqfX7aDUGFn033gYDRMHf5Sc/mO9GczsMecE8yiPlPaPBe478Jn4tCBV87D//y+D79OhP89jZjzF/8XzuXnY3m+o3ceuJt+73v133oEjdUYdCTg6Ew4PRXBGRjKbgLMObtW4+4OadbqaH/FFumrTX/2vvKdVO/GeY/U+ujOI3c3q9TZc/yDumnWV170DtBmZFItTljmcHMaqIsi3WztuLvkud34+xFjtmRM+fzQ0dFIYiFLRXU0ABwcIxdIyYzqRgAUeG8ntCcH4wn/xQfq8ArJKH4SGSFeHWE25lUuEk7nz7TlqiLdz12bsI+Pr+Fd5dqlH55KuwaBFcfPFgNVdEJGMpOEtm6myDxipo/Mh9bd4JxZNg8mlupbn758VDcXc5w4k3ETv132ho303Byz+jM7uAtlAejcEIu3PyeWPLs7xe+xoNzdWMrxzH1q426vx+6rpHhKuegapnCPvDFJUUURgqpDhcTGG4kBGhIi4NF1EULqI45PYVhYsoChWRH8zH7/Ont79kyDPGcMXhV5ATyOHWxbcyf/F8/u24f+v7H05eFd4RP/iVezBnzr7PExGRpCk4y8Gpahns3uKFYy8glxwCJ3m1wT+bDO0NvV7SfOjZrCkoYWfLTg6jk8a8YuoKy6n2WbYRZWnVs7zz0F+Ixjrxja8k1hNIWgn7LaWx3ZR2BagoOoT80cdycnYppQlbcbiYonAR2YHswe0LGVYunnYxVc1V3Pf+fcwqn8XZk87e53m9bg688Ua46qpBbKWISGZScJahIRZzQTe70D1f9Vd3E13iqHHBGLjEm9nhT9+AnasAsMZPR6SIbc3bee69XGpaa5gy9nB2d7WyyXawoauZrXSyu3Up9m+Xu9eH3CwGJeESSrPLKMkuYXR2KZ/yQnBJdgml4VLKImWUZpcSCURUEiFDxjeP/CZLqpfwkzd/wvGjj6c0u7T/F8Rig9MwEZEMp1k1ZOC1N0L9Nmitg3Gfdvveug82vOCFYm/LH4W98T0aOhrIevgSIh8uoj2cT3M4n93BbD6M5POH0YdQ01bDiLrN1HU0Uu33s8vvSxgdhrxgXnwkOOyF4D1Gh0uzSykMFapEQg5aG+o3cOETFzJ3/FxuO+m2vY5/2PAhZz1+Fsu/9r7b8cQTcPa+R6dFRCROs2rIwLHWTaW2e7NbZnna2eDzweJ74d3fuXKKtt3uVF+A965+jpq2OsasfYrSquXsDobZGcrio+yxbPJ18eD/HEVnrJN8G6Nl/BiiXiDO8sUoy86itLOZytxKSstmMiNxdDi7lDJv5DjkD6WzR0QGxcSCiXxp2pd4aNVDXDvzWsbkjel1fK9BkYKCQWydiEhmUnCW/sW6XJnE7i0uGE+ZB+F8WPYIvPIztz9hyrWnz7+LLbaDig//zviOOrbl5bExP5sNtoMdfj/vPH2ZWwIYMKXZFIWLeo0Cf8Urkeg1OhwpJS8rT6USInu4fPrlPPzBwzy44kF+cNwPeh3rrnFed82FHHL3Y1BWlo4miohkFAXn4a4r6gLxrk3u66RToaAC1jwDT38HGrZDLL6C3JOnfptlWX4KNy/myM46NuWG2ewLURUIsD0QYN07txM1hqJQEeXjD6M8Us6InBFMiJRzXGQEV+xxM11/02mJSP/KI+WcM+kcHl/7ONfNvI6icFHPse7g3FHi7VNwFhE5YEotw0FHiwvGuzZC+aFQPBG2v+sW6di9uVcwXnTSNSwuHIHZvpTZvk42FBWxiSjbvWC8Zf2jZAfzqcyrZN3McxiZM5KRkRF8ygvI5ZFyyiPlKpcQGSQXT72YhWsX8uKWFzlv8nnxA16lxsTfLHQPiosHv3EiIhlGwTlTdLRA7ToI5bpg3PgRPPZPLiw3VvWctmL2Zbw29lM07ljJSaaD9SXlfGDb2RzwUxUIUL3lSdieRUVuBasO/QyVuZWMyavkuNwKKvMqqcitoCCkWkmRoWJa8TRG5Yzi+c3P9wrO3SPO4Y9q3Q7f/pfpFhGR/ik4H0ysdcs9ByOuxOLZW6BmDdSsg/rNAGw5/FxeOmwu1fWbOHf3Oj7M8rG2pJQNPsuWrAAbd75Ac+1LlGaXsnTSLCryKqjMreRILxSPyRtDWXaZZpsQOUgYY5gzdg4L1iygpbOFSFak1/G28mLC516QptaJiGQWBeehbO3foWop7FwDtWuhZh3RQz7L2lO/z8b6jZyw/FEafD42ZPlZUVTE+iw/K3e/xda33iUSiPDaITOozHWB+Ii8Ss7KG0NFbgWjc0drkQ6RDDJnzBx+v+r3vL79dU4ddyoQn1Uj0NgMnZ3pbJ6ISMZQcE4na6Gp2i30sWOV23x+7Nk/Z1f7LkLPfI+cmrXUh/PZFo6wJjeHV3e+xjNPXgSAb0QelXljmFgwkQkFEzixYAKXFUxgbP5YikJFmoVCZJiYNWIW+cF8Xtz6Yjw4e6UagdZ2eOAB+O1v09hCEZHMoOA8WLqiULMaatbC9C8AYBdeiXl/Yc8pjcEIH+QUcNMjJ1HfXk9lOErduEpsMIfx+eOZUDCBQwomMNcLyuPyxxH0B9P1iURkiAj4Aswom8HK2pU9+7qD8/obv8KkE85JV9NERDKKgvNA2vQavL8AqpZhq1dgom0A3L7rJpY3bWJk3RIKiotYG8xiXTALk1PGxIKJzC2Y0DOKPLFgIiNyRuAzurFHRPo2rWgai6sW09nVSZY/q6dUY8NNlzNp3Olpbp2ISGZQcD5Q0Q6oXg5b34aqZVC1jN1f+AXLu5ox797HUSueYm04wrs5QVYFc1gZDFK96SkmF08lZ8bFTCyawrziKUwsmKjZKkTkE5taPJVoLMr6+vVMK57Ws9+gki0RkVRRcP64Gj8CfxAixbDhJezDF/WMJDdmhfkgFOK2J7/E2mCQgLWMmDaLqcXTmFI8hdOKpnJd0VQq8io0giwiKTW1eCoAq+tWM614Wk+phoKziEjqKDj3pyvqRpG3vglb3sRufRNTv5X3Z1/KguJyNm17gzmRLJaFc1kWChIoGMv00sM5p3QG00unM614GnnBvHR/ChEZBsbljSPsD/NB3Qecy7nxA8rNIiIpo+Dcj1i0BXPfaRgbY3coh6WhIG8WF/JK1fPUNhYzs2wmzVPO5Qslh3NL6XSKw1qZS0TSw+/zM7loMmt2rQHi09FpxFlEJHUUnPtxxQs3ECovYX0wC5s3mlkjZjFrxCx+Vj6LyUWTVW4hIkPK1OKpPLvpWay1KtUQERkACs79uGjKRUQnn8+R5UdSmVupeZFFZEibWjSVBWsWUN1SHQ/O+r0lIpIyCs79OHPimelugohI0kbnjgaguqWagNGvdxGRVFOtgYhIhijNLgWgprVGpRoiIgNAwVlEJEOUZZcBUNtaG785UKUaIiIpo+AsIpIhisJFGAw7W3emuykiIhlJwVlEJEMEfAGKwkUq1RARGSAKziIiGaQ0u7R3cFaphohIyig4i4hkkLLsMmpaarQAiojIAFBwFhHJICXZJdS01aS7GSIiGUnBWUQkg3SXasRsDNCIs4hIKik4i4hkkNLsUqKxKPXt9W6HcrOISMooOIuIZJDuuZy7p6TTiLOISOooOIuIZJCS7BLArR4ImlVDRCSVFJxFRDJI97LbGnEWEUk9BWcRkQzSE5xbFJxFRFJNwVlEJIPkZuUS9oe17LaIyABQcBYRySDGGIrDxdS21vY8FxGR1FBwFhHJMJGsCC3RlnQ3Q0Qk4yg4i4hkmLA/TGu0FVCNs4hIKik4i4hkmHAgTDQWBVSqISKSSgcUnI0xFxpjVhhjYsaY2Qn7xxtjWo0xS73tvxOOHWWMWW6MWWeM+YXRb3URkZTKDmT3PNaIs4hI6hzoiPP7wPnAy/s4tt5aO9Pbrk7YfzfwdWCyt807wDaIiEiCcCDc81hjEyIiqXNAwdlau8pauzrZ840xo4B8a+0b1loL/A74woG0QUREeksccRYRkdQZyBrnCcaYd40xLxljTvL2VQBbE87Z6u0TEZEUUamGiMjACOzvBGPMc8DIfRy6xVr7lz5eVgWMtdbWGmOOAv5sjJn+cRtnjLkKuApg7NixH/flIiLDUtgf3v9JIiLyse03OFtrT/u4b2qtbQfavcdvG2PWA1OAbUBlwqmV3r6+3ude4F6A2bNn24/bDhGR4Ug1ziIiA2NASjWMMWXGGL/3eCLuJsAN1toqoMEYc5w3m8blQF+j1iIi8gmoVENEZGAc6HR05xljtgLHA08ZY57xDp0MvGeMWQosAK621tZ5x64FfgOsA9YDfzuQNoiISG+9RpwVnEVEUma/pRr9sdY+Djy+j/0LgYV9vGYJcPiBfF8REelbrxFnlWqIiKSMVg4UEckwujlQRGRgKDiLiGQY1TiLiAwMBWcRkQyTWOOs3CwikjoKziIiGUYjziIiA0PBWUQkwyg4i4gMDAVnEZEMowVQREQGhoKziEiGSZxVQyPOIiKpo+AsIpJhet0cKCIiKaPgLCKSYSKBSLqbICKSkRScRUQyTMAXwG/8gGqcRURSScFZRCTDGGN6yjVU4ywikjoKziIiGah7SjoFZxGR1FFwFhHJQN0za6hUQ0QkdRScRUQykEo1RERST8FZRCQDJa4eKCIiqaHgLCKSgXqCswacRURSRsFZRCQDqVRDRCT1FJxFRDJQz82BCs4iIimj4CwikoF6pqPTrBoiIimj4CwikoFUqiEiknoKziIiGUgLoIiIpJ6Cs4hIBtJ0dCIiqafgLCKSgbpLNTTgLCKSOgrOIiIZSLNqiIiknoKziEgGKs4uxmd8KtkQEUkhBWcRkQx02tjTeOzsxyjJLkl3U0REMoaCs4hIBgr4AkwpmpLuZoiIZBQFZxERERGRJCg4i4iIiIgkQcFZRERERCQJCs4iIiIiIklQcBYRERERSYKCs4iIiIhIEhScRURERESSoOAsIiIiIpIEBWcRERERkSQoOIuIiIiIJEHBWUREREQkCQrOIiIiIiJJUHAWEREREUmCgrOIiIiISBIUnEVEREREkqDgLCIiIiKSBAVnEREREZEkKDiLiIiIiCRBwVlEREREJAkKziIiIiIiSTDW2nS3ISnGmJ3Ah2n41qVATRq+78FIfZUc9VPy1FfJUT8lT32VPPVVctRPyTtY+mqctbZsXwcOmuCcLsaYJdba2elux8FAfZUc9VPy1FfJUT8lT32VPPVVctRPycuEvlKphoiIiIhIEhScRURERESSoOC8f/emuwEHEfVVctRPyVNfJUf9lDz1VfLUV8lRPyXvoO8r1TiLiIiIiCRBI84iIiIiIklQcO6HMWaeMWa1MWadMeZf092eocQYs8kYs9wYs9QYs8TbV2yM+bsxZq33tSjd7UwHY8z9xpgdxpj3E/bts2+M8wvvGnvPGDMrfS0fXH30078bY7Z519VSY8yZCce+5/XTamPMGenkJrCCAAAE+ElEQVRpdXoYY8YYY14wxqw0xqwwxtzo7dd1laCfftJ1tQdjTNgY86YxZpnXV//X2z/BGLPY65NHjDFBb3/Ie77OOz4+ne0fTP301QPGmI0J19VMb/+w/PnrZozxG2PeNcY86T3PqGtKwbkPxhg/8Evgc8BhwJeMMYelt1VDzmettTMTppb5V+Af1trJwD+858PRA8C8Pfb11TefAyZ721XA3YPUxqHgAfbuJ4A7vetqprX2aQDvZ+8SYLr3ml95P6PDRRT4trX2MOA44DqvT3Rd9dZXP4Guqz21A3OstUcAM4F5xpjjgP+H66tDgF3Ald75VwK7vP13eucNF331FcC/JFxXS719w/Xnr9uNwKqE5xl1TSk49+0YYJ21doO1tgP4I3Bumts01J0LPOg9fhD4QhrbkjbW2peBuj1299U35wK/s84bQKExZtTgtDS9+uinvpwL/NFa226t3Qisw/2MDgvW2ipr7Tve40bc/5Qq0HXVSz/91Jdhe11510aT9zTL2ywwB1jg7d/zmuq+1hYApxpjzCA1N6366au+DMufPwBjTCXweeA33nNDhl1TCs59qwC2JDzfSv+/gIcbCzxrjHnbGHOVt2+EtbbKe/wRMCI9TRuS+uobXWd7u9778+b9Jl7uo37yeH/OPBJYjK6rPu3RT6Drai/en9SXAjuAvwPrgd3W2qh3SmJ/9PSVd7weKBncFqfPnn1lre2+ruZ719WdxpiQt284X1d3Ad8FYt7zEjLsmlJwlk/qRGvtLNyfpK4zxpyceNC66Vo0Zcs+qG/6dTcwCffn0Crg9vQ2Z2gxxuQCC4FvWWsbEo/puorbRz/putoHa22XtXYmUIkbaZ+W5iYNWXv2lTHmcOB7uD47GigGbk5jE9POGHMWsMNa+3a62zKQFJz7tg0Yk/C80tsngLV2m/d1B/A47pdudfefo7yvO9LXwiGnr77RdZbAWlvt/Q8qBvya+J/Nh30/GWOycGHw99baP3m7dV3tYV/9pOuqf9ba3cALwPG4soKAdyixP3r6yjteANQOclPTLqGv5nmlQdZa2w78Fl1XJwDnGGM24cpb5wA/J8OuKQXnvr0FTPbuBg3ibiB5Is1tGhKMMTnGmLzux8Bc4H1c/3zVO+2rwF/S08Ihqa++eQK43LsL+zigPuFP78POHnWA5+GuK3D9dIl3F/YE3E03bw52+9LFq/u7D1hlrb0j4ZCuqwR99ZOuq70ZY8qMMYXe42zgdFxN+AvAF73T9rymuq+1LwLP22GyEEQfffVBwj9aDa5uN/G6GnY/f9ba71lrK62143GZ6Xlr7aVk2DUV2P8pw5O1NmqMuR54BvAD91trV6S5WUPFCOBxr4Y/ADxsrf1fY8xbwKPGmCuBD4GL0tjGtDHG/AE4BSg1xmwFfgT8hH33zdPAmbibklqAfxr0BqdJH/10ijelkwU2Ad8AsNauMMY8CqzEzZxwnbW2Kx3tTpMTgMuA5V6dJcD30XW1p7766Uu6rvYyCnjQm0XEBzxqrX3SGLMS+KMx5lbgXdw/RPC+PmSMWYe7qfeSdDQ6Tfrqq+eNMWWAAZYCV3vnD9efv77cTAZdU1o5UEREREQkCSrVEBERERFJgoKziIiIiEgSFJxFRERERJKg4CwiIiIikgQFZxERERGRJCg4i4iIiIgkQcFZRERERCQJCs4iIiIiIkn4/71SFPpKWekmAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 864x648 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"# Read h5 file, generated with a version that computes azimuth angles in [0, 360] (2d2a16293c7ae54d6de2d4163f3d23718fc9646f)\n", | |
"h5 = h5py.File('/data/avhrr_gac/pygac_h5/ECC_GAC_sunsatangles_tirosn_99999_19781125T1642148Z_19781125T1721013Z.h5', mode='r')\n", | |
"azi_h5 = h5.get('/image5/data')[:]\n", | |
"azi_h5_corr = centered_modulus(decomp(azi_h5, 0), 360) # correct encoding\n", | |
"azi_h5_wrong = centered_modulus(decomp(azi_h5, 180), 360) # incorrect encoding\n", | |
"\n", | |
"# Read nc file generated with current version (no h5 files involved)\n", | |
"nc = xr.open_dataset('/data/avhrr_gac/output/AVHRR-GAC_FDR_1C_TSN_19781125T164207Z_19781125T172101Z_R_O_20200101T000000Z_0100.nc')\n", | |
"azi_nc = nc['sensor_azimuth_angle']\n", | |
"\n", | |
"# Plot azimuth \n", | |
"plt.figure(figsize=(12,9))\n", | |
"plt.plot(azi_nc.isel(y=0), label='nc', color='C2')\n", | |
"plt.plot(azi_h5_corr[0, :], label='h5 corr enc', linestyle='--', color='C1')\n", | |
"plt.plot(azi_h5_bad[0, :], label='h5 wrong enc', linestyle='--', color='red')\n", | |
"plt.legend()\n", | |
"plt.show()" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The red line is with the wrong offset applied. The orange line is with the correct offset applied. But you can clearly see the clipping for large values. I'd say the current version (green line) is correct, you could say even better than before :)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Luckily, both problems have been fixed by https://github.com/pytroll/pygac/pull/35 and https://github.com/pytroll/pygac/pull/73" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "pygac-fdr", | |
"language": "python", | |
"name": "pygac-fdr" | |
}, | |
"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.2" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 4 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment