George J. Bendo
\n",
"02 January 2022
This script takes three FITS images covering the same part of the sky, colourizes each one using a different colour (red, green, or blue), and then combines the three images together to produce a colour image that is exported as a png file. This type of image is called a representative colour image because each colour represents a specific wavelength of light. However, the colours may not be the same as the observed colours, especially when the colours correspond to wavelengths of light outside the visible part of the electromagnetic spectrum (as is the case for the images used in this demonstration scripts, which are all infrared images).
\n", "\n", "The example used in this script is based on the following files:
\n", "These files should be downloaded to the directory containing the Juptier Notebook before attempting to execute the script.
\n", "\n", "This script could be adapted to work with other FITS files, although most of the lines would need to be modified. Modifications would therefore only be suggested for people who are familiar with Python programming or who have studied this script carefully.
\n", "\n", "The following packages need to be installed to use this script: astropy, matplotlib, numpy, and PIL. The matplotlib and numpy packages are standard python utilities. The astropy package contains the tools needed to import FITS files. FITS_tools is used to shift all of the images to the same coordinate system. PIL is used to export images as png files.
\n" ] }, { "cell_type": "markdown", "id": "4f60d5ee", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "id": "b1393d31", "metadata": {}, "source": [ "Perform a series of prepratory steps first." ] }, { "cell_type": "code", "execution_count": 1, "id": "6aec7ed9", "metadata": {}, "outputs": [], "source": [ "# Import packages.\n", "import numpy\n", "import matplotlib.pyplot as pp\n", "from astropy.io import fits\n", "from FITS_tools import hcongrid\n", "from PIL import Image" ] }, { "cell_type": "markdown", "id": "4e7aa4df", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "id": "20bc1f21", "metadata": {}, "source": [ "Read in and prepare the images." ] }, { "cell_type": "code", "execution_count": 2, "id": "48438234", "metadata": {}, "outputs": [], "source": [ "# Read the image files. Each fits image is assigned to a different colour \n", "# (imgr, imgg, and imgb).\n", "file=fits.open('ngc3031_spitzer3.6micron.fits')\n", "imgb=file[0].data\n", "headb=file[0].header\n", "file.close()\n", "\n", "file=fits.open('NGC_3031_I_MIPS24_bgm2012.fits')\n", "imgg=file[0].data\n", "headg=file[0].header\n", "file.close()\n", "\n", "file=fits.open('ngc3031_herschel250micron.fits')\n", "imgr=file[0].data\n", "headr=file[0].header\n", "file.close()" ] }, { "cell_type": "code", "execution_count": 3, "id": "f0c33bba", "metadata": {}, "outputs": [], "source": [ "# Use congrid to match the coordinate systems of all three images to \n", "# the coordinate system of the blue image (which has the smallest \n", "# pixels in this case).\n", "imggshift=hcongrid.hastrom(imgg,headg,headb,preserve_bad_pixels=True)\n", "imgrshift=hcongrid.hastrom(imgr,headr,headb,preserve_bad_pixels=True)" ] }, { "cell_type": "code", "execution_count": 4, "id": "803c0fba", "metadata": {}, "outputs": [], "source": [ "# Crop the images.\n", "ix1=747\n", "ix2=2187\n", "iy1=837\n", "iy2=2997\n", "imgbcrop=imgb[iy1:iy2,ix1:ix2]\n", "imggcrop=imggshift[iy1:iy2,ix1:ix2]\n", "imgrcrop=imgrshift[iy1:iy2,ix1:ix2]" ] }, { "cell_type": "markdown", "id": "c009d013", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "id": "874e8d4d", "metadata": {}, "source": [ "Rescale the images." ] }, { "cell_type": "code", "execution_count": 5, "id": "54df8e0b", "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\georg\\AppData\\Local\\Temp/ipykernel_15108/846428621.py:13: RuntimeWarning: invalid value encountered in log10\n", " imgglog=numpy.log10(imggcrop+0.13)\n", "C:\\Users\\georg\\AppData\\Local\\Temp/ipykernel_15108/846428621.py:14: RuntimeWarning: invalid value encountered in log10\n", " imgrlog=numpy.log10(imgrcrop+3.0)\n" ] } ], "source": [ "# Convert the image to log values, which will reveal more detail when the \n", "# images are converted to png. \n", "\n", "# These specific image are either background-subtracted or otherwise have \n", "# pixel values equal to or less than 0. To make the colour scaling look \n", "# more natural, it is appropriate to add an offset to the data, which is \n", "# based on 3 or 5 times the background noise for these examples (or, in\n", "# statistical terms, 3 or 5 times the standard deviation of the background\n", "# pixel values.)\n", "#\n", "# [NOTE: Ignore the RuntimeWarning.]\n", "imgblog=numpy.log10(imgbcrop+0.90)\n", "imgglog=numpy.log10(imggcrop+0.13)\n", "imgrlog=numpy.log10(imgrcrop+3.0)" ] }, { "cell_type": "code", "execution_count": 6, "id": "875a0c81", "metadata": {}, "outputs": [], "source": [ "# Rescale the images so that the pixel values run from 0 to 255. The numbers\n", "# in this step could be adjusted to either brighten or dim the red, green, or\n", "# blue parts of the image.\n", "imgbscl=(imgblog+0.0)*2.0*255\n", "imggscl=(imgglog+0.5)*1.0*255\n", "imgrscl=(imgrlog-0.55)*1.1*255" ] }, { "cell_type": "code", "execution_count": 7, "id": "590317f9", "metadata": {}, "outputs": [], "source": [ "# Set out-of-range values in the rescaled images to fall between 0 \n", "# and 255.\n", "imgbscl[numpy.isnan(imgbscl)]=0\n", "imggscl[numpy.isnan(imggscl)]=0\n", "imgrscl[numpy.isnan(imgrscl)]=0\n", "imgbscl[imgbscl<0]=0\n", "imggscl[imggscl<0]=0\n", "imgrscl[imgrscl<0]=0\n", "imgbscl[imgbscl>255]=255\n", "imggscl[imggscl>255]=255\n", "imgrscl[imgrscl>255]=255" ] }, { "cell_type": "code", "execution_count": 8, "id": "6b72b94f", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAJQCAYAAABhOHhJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApyklEQVR4nO3deZhkd1no8e+bmUyWySTDzIQISTAQAhpAAowRL3hlEwGRiLgEFxAiAcUlCjyA+3MF4SKLgj5oWEzwQlgEJcomIosiCJMIhATQCMEkBKpnQpLJZJvlvX+cUz1F01Nd3XWWOqe+n+eZp6tOLf2maPrbp6p+pyIzkSQJ4LC2B5AkzQ6jIElaZBQkSYuMgiRpkVGQJC1a3/YA09i2bVuecsopbY8hSZ1yySWX7MzM45e7rNNROOWUU9ixY0fbY0hSp0TEVw91mU8fSZIWGQVJ0iKjIElaZBQkSYuMgiRpkVGQJC0yCpKkRUZBkrTIKEiSFhkFSdIioyBJWmQUJEmLjIIkaZFRkCQtMgqSpEVGQZK0yChIkhYZBUnSIqMgSVpkFCRJi4yCJGmRUZAkLVrf9gCq19698Nznwq5dbU8iqUqPfzycfXb192sUeu7yy+HVr4bv+A7YuLHtaSRV5f73r+d+jULPDQbF13e8Ax760HZnkTT7fE2h54ZRuPOd251DUjcYhZ4zCpJWwyj03GAAhx8Oxx3X9iSSusAo9NxgUOwlRLQ9iaQuMAo9N4yCJE2itihExMkR8eGIuCIiLo+IXy+3/0FEXBsRnyn/PW7kNi+MiCsj4ksR8cN1zTZPjIKk1ajzLan7gOdk5qURsQm4JCI+WF72qsx8+eiVI+J04GzgPsBdgX+KiHtl5v4aZ+y9wQC+67vankJSV9S2p5CZ12XmpeXp3cAXgBPH3OQs4K2ZeXtmfgW4EjizrvnmhXsKklajkdcUIuIU4AHAv5ebfiUiPhcRb4yIO5XbTgSuHrnZNSwTkYg4NyJ2RMSOhYWFOsfuvD174NZbjYKkydUehYg4BngncF5m3gS8FjgVOAO4DnjFau4vM8/PzO2Zuf3444+vetxeGa5R8GGSNKlaoxARh1ME4c2Z+S6AzPxGZu7PzAPA6zj4FNG1wMkjNz+p3KY1cuGapNWq891HAbwB+EJmvnJk+11GrvZE4PPl6YuBsyPiiIi4O3Aa8Km65psHRkHSatX57qOHAD8PXBYRnym3/Rbw5Ig4A0jgKuCZAJl5eUS8HbiC4p1Lz/adR9MxCpJWq7YoZOa/Asuto33vmNu8GHhxXTPNG19TkLRarmjuscEAjjkGjj667UkkdYVR6DHXKEhaLaPQY0ZB0moZhR5bWDAKklbHKPTYYOCLzJJWxyj01IED7ilIWj2j0FM33AD79hkFSatjFHrKhWuS1sIo9JRRkLQWRqGnjIKktTAKPWUUJK2FUeipYRS2bWt3DkndYhR6ajCArVthfZ3HwZXUO0ahp1yjIGktjEJPuZpZ0loYhZ7yYHiS1sIo9JRRkLQWRqGH9u6F6683CpJWzyj00M6dxVejIGm1jEIPuXBN0loZhR4yCpLWyij0kFGQtFZGoYcWFoqvRkHSahmFHhoMisNbbN7c9iSSusYo9NBwNXNE25NI6hqj0EMuXJO0Vkahh4yCpLUyCj1kFCStlVHoIaMgaa2MQs/s2VP8MwqS1sIo9IxrFCRNwyj0jKuZJU3DKPSMewqSpmEUema4p+BHcUpaC6PQMz59JGkaRqFnBgM4+mjYuLHtSSR1kVHoGdcoSJqGUegZoyBpGkahZ4yCpGkYhZ4xCpKmYRR6JLNYp2AUJK2VUeiRG2+EvXuNgqS1Mwo94hoFSdMyCj3iamZJ0zIKPeKegqRpGYUeMQqSpmUUemQYhW3b2p1DUncZhR4ZDOBOd4ING9qeRFJXGYUeceGapGkZhR4xCpKmZRR6xChImpZR6BEPcSFpWkahJ/btg127XLgmaTpGoSd27SoOiOeegqRpGIWecOGapCoYhZ4wCpKqYBR6wihIqoJR6AmjIKkKRqEnBgNYt644zIUkrZVR6InBoHg76mH+LyppCv4K6QkXrkmqglHoCQ9xIakKRqEnhk8fSdI0jEJPuKcgqQpGoQduvRV27zYKkqZnFHpgYaH4ahQkTcso9IAL1yRVxSj0gFGQVBWj0ANGQVJVjEIPGAVJVTEKPbCwAEcdBRs3tj2JpK4zCj0wXLgW0fYkkrrOKPSAC9ckVcUo9IBRkFQVo9ADRkFSVYxCx2UaBUnVMQodd9NNcMcdRkFSNYxCx7lGQVKVjELHGQVJVTIKHecRUiVVySh0nHsKkqpkFDpuGIVt29qdQ1I/GIWOGwzguOPgiCPankRSHxiFjhse90iSqmAUOm4wgBNOaHsKSX1hFDrO1cySqmQUOs4oSKpSbVGIiJMj4sMRcUVEXB4Rv15u3xIRH4yI/yq/3qncHhHx6oi4MiI+FxEPrGu2vti/H3buNAqSqlPnnsI+4DmZeTrwYODZEXE68ALgQ5l5GvCh8jzAY4HTyn/nAq+tcbZe2LWrOCCeUZBUldqikJnXZeal5endwBeAE4GzgAvLq10I/Fh5+izgTVn4JLA5Iu5S13x94GpmSVVr5DWFiDgFeADw78AJmXldedHXgeF7Z04Erh652TXlNh2Cq5klVa32KETEMcA7gfMy86bRyzIzgVzl/Z0bETsiYsfC8E/lOTWMgusUJFWl1ihExOEUQXhzZr6r3PyN4dNC5dfyVxvXAieP3Pykctu3yMzzM3N7Zm4/fs5/G7qnIKlqdb77KIA3AF/IzFeOXHQx8NTy9FOBd49sf0r5LqQHAzeOPM2kZQwGcNhhsGVL25NI6ov1Nd73Q4CfBy6LiM+U234LeCnw9og4B/gq8FPlZe8FHgdcCdwCPK3G2XphMCgOhLduXduTSOqL2qKQmf8KxCEufuQy10/g2XXN00cuXJNUNVc0d5hRkFQ1o9BhRkFS1YxChy0sGAVJ1TIKHXX77XDjjUZBUrWMQkd5iAtJdTAKHeXCNUl1MAod5SEuJNXBKHSUewqS6mAUOsooSKqDUeiowQCOOAI2bWp7Ekl9YhQ6arhwLQ51IBFJWgOj0FGuZpZUB6PQUa5mllQHo9BR7ilIqoNR6KBMoyCpHkahg26+GW67zShIqp5R6CBXM0uqi1HoIBeuSaqLUeggoyCpLkahg4yCpLoYhQ7yNQVJdTEKHbSwAMceC0ce2fYkkvrGKHSQaxQk1cUodJBRkFQXo9BBRkFSXYxCBxkFSXUxCh1z4EDxQrPvPJJUB6PQMddfX4TBPQVJdTAKHePCNUl1MgodYxQk1ckodIxRkFQno9AxCwvFV6MgqQ5GoWMGA4iArVvbnkRSHxmFjhkMYNs2WLeu7Ukk9ZFR6BgXrkmqk1HoGKMgqU5GoWOMgqQ6GYWOGQw8xIWk+hiFDrnjDrjhBvcUJNXHKHSIaxQk1c0odIhRkFQ3o9AhHuJCUt2MQocYBUl1MwodYhQk1c0odMhgABs2wLHHtj2JpL4yCh0yXLgW0fYkkvrKKHSIq5kl1c0odIirmSXVzSh0iHsKkupmFDrEKEiqm1HoiD174NZbjYKkehmFjnCNgqQmGIWOMAqSmmAUOsIoSGqCUegIoyCpCUahI4ZRcJ2CpDoZhY4YDGDTJjjqqLYnkdRnRqEjXKMgqQlGoSM8xIWkJhiFjlhYcE9BUv2MQkf49JGkJhiFDjhwwD0FSc0wCh1www2wb59RkFQ/o9ABLlyT1BSj0AFGQVJTjEIHGAVJTTEKHWAUJDXFKHTAYAARsHVr25NI6juj0AGDAWzZAuvXtz2JpL4zCh3gGgVJTTEKHeBqZklNMQodYBQkNcUodIBRkNQUozDj9u6F6683CpKaYRRm3M6dxVejIKkJRmHGuXBNUpOMwowzCpKaZBRmnFGQ1CSjMOMWFoqvRkFSEyaKQkQ8NCKeVp4+PiLuXu9YGhoM4PDD4bjj2p5E0jxYMQoR8fvA84EXlpsOB/5fnUPpoMEAjj++OCCeJNVtkj2FJwJPAPYAZObXgE11DqWDXLgmqUmTROGOzEwgASJiY70jaZRRkNSkSaLw9oj4S2BzRDwD+CfgdfWOpSGjIKlJKx6hPzNfHhE/BNwE3Bv4vcz8YO2TCTAKkpo10ce2lBEwBA275RbYs6d4oVmSmrBiFCJiN+XrCcAGincf7cnMY+scTAcXrp1wQrtzSJofkzx9tPhOo4gI4CzgwXUOpYKrmSU1bVUrmrPwd8APr3TdiHhjRAwi4vMj2/4gIq6NiM+U/x43ctkLI+LKiPhSRKx4//PA1cySmjbJ00c/PnL2MGA7cNsE930B8GfAm5Zsf1VmvnzJ9zgdOBu4D3BX4J8i4l6ZuX+C79Nbwz0FX1OQ1JRJXmj+0ZHT+4CrKJ5CGiszPxYRp0w4x1nAWzPzduArEXElcCbwiQlv30tGQVLTJnlN4WkVf89fiYinADuA52TmN4ETgU+OXOeacttcGwxg48binyQ14ZBRiIjXcPBdR98mM39tDd/vtcAflvf7h8ArgKev5g4i4lzgXIC73e1uaxihO1yjIKlp4/YUdlT9zTLzG8PTEfE64B/Ks9cCJ49c9aRy23L3cT5wPsD27dsPGa0+MAqSmnbIKGTmhVV/s4i4S2ZeV559IjB8Z9LFwFsi4pUULzSfBnyq6u/fNYMBnHzyyteTpKpM8u6j4ykOnX06cORwe2Y+YoXbXQQ8DNgWEdcAvw88LCLOoHj66CrgmeV9XR4RbweuoHgx+9nz/s4jKKLwoAe1PYWkeTLJu4/eDLwN+BHgWcBTgYWVbpSZT15m8xvGXP/FwIsnmGcuZPr0kaTmTbJ4bWtmvgHYm5kfzcynA2P3EjS9G26AffuMgqRmTbKnsLf8el1E/AjwNWBLfSMJXKMgqR3j3pJ6eGbuBV4UEccBzwFeAxwL/EZD880tD3EhqQ3j9hSujYiLgYuAmzLz88DDmxlLHgxPUhvGvabw3cCngd8Bro6IP40Ij47aEKMgqQ2HjEJm7srMv8zMh1Mch+jLwKsi4r8jwncJ1WwYhW3b2p1D0nyZ6NDZmfk1ireTvhbYDfxinUOpiMKWLXD44W1PImmejI1CRBwZET8ZEe8CrqR4K+oLKFYdq0auUZDUhnHvPnoL8CjgoxQL2H4mMyf5HAVVwChIasO4dx+9H3hmZu5uahgdNBjAfe7T9hSS5s24F5rfZBDa456CpDas6jOa1Yx9+2DXLlczS2qeUZhBu3YVX91TkNS0cS80//i4G2bmu6ofR+DCNUntGfdC84+OuSwBo1AToyCpLeM+ee1pTQ6ig4yCpLas+JpCRJwQEW+IiPeV50+PiHPqH21+GQVJbZnkheYLgA9wcBXzfwLn1TSPKKKwfj1s3tz2JJLmzSRR2JaZbwcOAGTmPmDuPz+5ToNB8XbUw3xvmKSGTfJrZ09EbKV4cZny8Nk31jrVnHPhmqS2TPJxnM8BLgZOjYiPA8cDP1nrVHNuuKcgSU1bMQqZeUlE/CBwbyCAL5VfVZPBAO5xj7ankDSPJnn30UeAkzLz8vIjOc+g+EQ21cSnjyS1ZZKnj14CvD8iXg2cCDwOcA1DTW69FW6+2ShIasckTx99ICKeBXwQ2Ak8IDO/Xvtkc2phofhqFCS1YZKnj34XeA3wv4E/AD4SET9S81xzy4Vrkto0ydNHW4EzM/NW4BMR8X7g9cB7ap1sThkFSW2a5Omj85ac/yrwQ3UNNO+MgqQ2jTt09p9k5nkR8feUC9dGZeYTap1sThkFSW0at6fw1+XXlzcxiAqDARx9NGzc2PYkkubRuChcHhHnAfcELgPeUB73SDVyNbOkNo1799GFwHaKIDwWeEUjE805F65JatO4PYXTM/N+ABHxBuBTzYw03wYDuOtdV76eJNVh3J7C3uEJnzZqzsKCewqS2jNuT+H+EXFTeTqAo8rzAWRmHlv7dHMm06ePJLVr3Gc0r2tyEMFNN8EddxgFSe3xs71miGsUJLXNKMwQoyCpbUZhhhgFSW0zCjNkGAUXr0lqyyGjEBHnRMTzRs5fGxE3RcTu8vMVVDGjIKlt4/YUngW8ceT8oHwb6vHAk2udak4NBrB5M2zY0PYkkubVuChEZu4aOf8OgMy8DTiq1qnmlGsUJLVtXBQ2j57JzD8CiIjDgG01zjS3jIKkto2Lwj9GxIuW2f5/gH+saZ655iEuJLVtXBSeB5waEVdGxDvLf1dSHEr7uc2MN1/cU5DUtnGHudgDPDki7gHcp9x8RWb+dyOTzZn9+2HnTqMgqV3jPo7zbuXJfcBnl27PzP+pd7T5smtXcUA8oyCpTeOOkvoeis9mjpFtSfGW1DsDHjCvQq5mljQLxj19dL/R8xFxCvB84FHAH9U71vxx4ZqkWbDiYS4i4rSIuAB4H3AJxSeyvabuweaNewqSZsG41xTuC/w2xYvMLwPOycz9TQ02b4yCpFkw7jWFzwJXU7y2cCZwZsTBlxcy89fqHW2+DAZw2GGwZUvbk0iaZ+OicA7FC8tqwGBQvJ5wmMetldSicS80XxARxwPfCVyZmTc0NtUccjWzpFkw7tDZvwhcDrwG+GJEPKGxqeaQq5klzYJxT1acB9wnM78f+F/ACxuZaE4ZBUmzYFwU7sjMBYDM/DJwRDMjzSejIGkWjHuh+aSIePWhzvvuo+rcdhvcdJML1yS1b1wUnrfk/CV1DjLPFhaKr+4pSGrbuHcfXdjkIPPMhWuSZoXvip8BRkHSrDAKM8AoSJoVkxwQ7yGTbNPaGQVJs2KSPYXljojqUVIrNBjAkUfCMce0PYmkeTfuKKnDRWvHR8Rvjlx0LH7ATqWGh7iIWPm6klSncW9J3QAcU15n08j2m4CfqHOoeePCNUmzYtxbUj8KfDQiLsjMrzY409wZDOCEE9qeQpLG7ykMXRAR33YI7cx8RA3zzKXBAO5737ankKTJovDckdNHAk8C9tUzzvzJ9OkjSbNjxShk5tLDW3w8Ij5V0zxzZ/duuP12oyBpNqwYhYgY/YDIw4AHAcfVNtGccY2CpFkyydNHl1B8LGdQPG30FYqP6lQFjIKkWTLJ00d3b2KQeWUUJM2SSZ4+OhL4ZeChFHsM/wL8RWbeVvNsc8EoSJolkzx99CZgNwcPbfEzwF8DP1nXUPNk+FkKfsCOpFkwSRTum5mnj5z/cERcUddA82YwgOOOgyP8sFNJM2CSA+JdGhEPHp6JiO8DdtQ30nwZDNxLkDQ7JtlTeBDwbxHxP+X5uwFfiojLgMzM76ltujngwjVJs2SSKDym9inm2GAA97xn21NIUmGSp49elJlfHf03uq3uAfvOPQVJs2SSKNxn9ExErKd4SklT2r8fdu40CpJmxyGjEBEvjIjdwPdExE0Rsbs8/w3g3Y1N2GPXXw8HDhgFSbPjkFHIzJdk5ibgjzPz2MzcVP7bmpkvbHDG3nLhmqRZM8kLze+LiP+9dGNmfqyGeeaKUZA0ayaJwvNGTh8JnElxkDw/ZGdKRkHSrJnkgHg/Ono+Ik4G/qSugebJ8BAXRkHSrJjk3UdLXQN8d9WDzKPBACJgy5aVrytJTZjkKKmvoTg6KhQROQO4dILbvRF4PDDIzPuW27YAbwNOAa4CfiozvxkRAfwp8DjgFuAXMnPF79F1gwFs2wbr1rU9iSQVJtlT2EHxGsIlwCeA52fmz01wuwv49tXQLwA+lJmnAR8qzwM8Fjit/Hcu8NoJ7r/zXLgmadZM8kLz24DhgRiunPRzFDLzYxFxypLNZwEPK09fCHwEeH65/U2ZmcAnI2JzRNwlM6+b5Ht1lVGQNGvGLV5bHxEvo3gN4UKKz1W4OiJeFhGHr/H7nTDyi/7rwAnl6ROBq0eud025bbm5zo2IHRGxY2H4Sm1HGQVJs2bc00d/DGwB7p6ZD8rMBwKnApuBl0/7jcu9glzxit9+u/Mzc3tmbj++48ecNgqSZs24KDweeEZm7h5uyMybgF+ieEF4Lb4REXcBKL+W79TnWuDkkeudVG7rrdtvhxtvNAqSZsu4KGT51/zSjftZw1/4pYuBp5ann8rBYyhdDDwlCg8Gbuz76wmuUZA0i8ZF4YqIeMrSjRHxc8AXV7rjiLiI4t1K946IayLiHOClwA9FxH8BjyrPA7wX+DJwJfA64JdX9V/RQa5mljSLxr376NnAuyLi6RRvRwXYDhwFPHGlO87MJx/iokcuc90sv9/ccE9B0iw6ZBQy81rg+yLiERz8TIX3ZuaHGpms54Z7Ch1/rVxSz0xy7KN/Bv65gVnmik8fSZpFazn2kSowGMCGDXDssW1PIkkHGYWWDNcoRLQ9iSQdZBRa4sI1SbPIKLTEKEiaRUahJUZB0iwyCi3INAqSZpNRaMGePXDbbUZB0uwxCi1wjYKkWWUUWmAUJM0qo9ACD3EhaVYZhRa4pyBpVhmFFrinIGlWGYUWDAawaRMcdVTbk0jStzIKLXCNgqRZZRRaYBQkzSqj0AKjIGlWGYUWLCwYBUmzySg07MABoyBpdhmFhn3zm7B/v1GQNJuMQsNcoyBplhmFhrmaWdIsMwoNMwqSZplRaJhRkDTLjELDBgOIgK1b255Ekr6dUWjYYFAEYf36tieRpG9nFBrmGgVJs8woNMxDXEiaZUahYUZB0iwzCg0zCpJmmVFo0B13FIe5MAqSZpVRaNDOncVXD3EhaVYZhQa5cE3SrDMKDTIKkmadUWiQUZA064xCg4yCpFlnFBq0sACHHw7HHdf2JJK0PKPQoOEahYi2J5Gk5RmFBrlwTdKsMwoNMgqSZp1RaJBRkDTrjEKDBgNXM0uabUahIXv2wC23uKcgabYZhYa4RkFSFxiFhhgFSV1gFBqysFB8NQqSZplRaIh7CpK6wCg0ZBgF330kaZYZhYYMBnDMMXD00W1PIkmHZhQa4sI1SV1gFBpiFCR1gVFoiFGQ1AVGoSEe4kJSFxiFBhw4UKxTcE9B0qwzCg244QbYt88oSJp9RqEBrmaW1BVGoQGuZpbUFUahAUZBUlcYhQYYBUldYRQaMIzCtm3tziFJKzEKDRgMYOtWWL++7UkkaTyj0ABXM0vqCqPQAFczS+oKo9AA9xQkdYVRaICHuJDUFUahZvv2wa5dRkFSNxiFmu3cWXw1CpK6wCjUzIVrkrrEKNTMKEjqEqNQM6MgqUuMQs2MgqQuMQo1GwyKw1ts3tz2JJK0MqNQs+HCtYi2J5GklRmFmnmIC0ldYhRq5mpmSV1iFGrmcY8kdYlRqJlRkNQlRqFGt9wCN99sFCR1h1Go0cJC8dUoSOoKo1AjF65J6hqjUCOjIKlrjEKNjIKkrjEKNRpGwcVrkrrCKNRoYQGOPho2bmx7EkmajFGokWsUJHWNUaiRUZDUNUahRkZBUte0EoWIuCoiLouIz0TEjnLbloj4YET8V/n1Tm3MViWjIKlr2txTeHhmnpGZ28vzLwA+lJmnAR8qz3dWpofNltQ9s/T00VnAheXpC4Efa2+U6d14I+zdCyec0PYkkjS5tqKQwD9GxCURcW657YTMvK48/XVg2V+nEXFuROyIiB0Lw4MLzSAXrknqovUtfd+HZua1EXFn4IMR8cXRCzMzIyKXu2Fmng+cD7B9+/ZlrzMLjIKkLmplTyEzry2/DoC/Bc4EvhERdwEovw7amK0qRkFSFzUehYjYGBGbhqeBRwOfBy4Gnlpe7anAu5uerUrDZ7Z8oVlSl7Tx9NEJwN9GxPD7vyUz3x8RnwbeHhHnAF8FfqqF2Soz3FPYtq3dOSRpNRqPQmZ+Gbj/Mtt3AY9sep66DAZwpzvBhg1tTyJJk5ult6T2igvXJHWRUaiJUZDURUahJkZBUhcZhZoYBUldZBRqsG8f7NplFCR1j1Gowa5dxQHxjIKkrjEKNRguXDMKkrrGKNRguHDN1cySusYo1MDjHknqKqNQA6MgqauMQg0GA1i3rjjMhSR1iVGowfBjOA/z0ZXUMf7aqoEL1yR1lVGogVGQ1FVGoQZGQVJXGYUaGAVJXWUUKnbbbbB7twvXJHWTUaiYh7iQ1GVGoWIuXJPUZUahYkZBUpcZhYoZBUldZhQqZhQkdZlRqNhgAEcdBRs3tj2JJK2eUajYcI1CRNuTSNLqGYWKuXBNUpcZhYotLBgFSd1lFCo2PGy2JHWRUahQpk8fSeo2o1Ch3bvh9tuNgqTuMgoVco2CpK4zChUyCpK6zihUyChI6jqjUCGjIKnrjEKFhlHwLamSusooVGgwgM2bYcOGtieRpLUxChVaWHAvQVK3GYUKuXBNUtcZhQoZBUldZxQqZBQkdZ1RqMj+/bBzp1GQ1G1GoSLXXw8HDhgFSd1mFCriwjVJfWAUKmIUJPWBUaiIUZDUB0ahIkZBUh8YhYosLMBhh8GWLW1PIklrZxQqMhjAtm1FGCSpq/wVVhEXrknqA6NQEaMgqQ+MQkWMgqQ+MAoVMQqS+sAoVOD22+HGG42CpO4zChVYWCi+GgVJXWcUKuDCNUl9YRQq4J6CpL4wChUY7in4+cySus4oVMCnjyT1hVGowGAARxwBmza1PYkkTccoVGC4RiGi7UkkaTpGoQIuXJPUF0ahAkZBUl8YhQoYBUl9YRSmlGkUJPWHUZjSzTfDbbcZBUn9YBSmNFzN7MI1SX1gFKbkwjVJfWIUpmQUJPWJUZiSUZDUJ0ZhSh4MT1KfGIUpDQZw7LFw5JFtTyJJ0zMKU3KNgqQ+MQpTMgqS+sQoTMkoSOoTozClhQWjIKk/jMIUDhwoouA7jyT1hVGYwje/Cfv3u6cgqT+MwhRcuCapb4zCFIyCpL4xClMwCpL6xihMwShI6hujMIXBACJg69a2J5GkahiFKQwGsG0brFvX9iSSVA2jMAVXM0vqG6MwBReuSeobozAF9xQk9Y1RmIJRkNQ3RmGN7rijOMyFUZDUJ0ZhjXbuLL4aBUl9MnNRiIjHRMSXIuLKiHhB2/McigvXJPXRTEUhItYBfw48FjgdeHJEnN7uVMszCpL6aH3bAyxxJnBlZn4ZICLeCpwFXFHlN/niF+Hii6e7j89+tvhqFCT1yaxF4UTg6pHz1wDfN3qFiDgXOBfgbne725q+yWWXwfOfv8YJR2zbBieeOP39SNKsmLUorCgzzwfOB9i+fXuu5T6e+ETYs2f6WTZsgPWdewQl6dBm7VfatcDJI+dPKrdVav16f5lL0nJm6oVm4NPAaRFx94jYAJwNTPnsvyRpUjP193Jm7ouIXwE+AKwD3piZl7c8liTNjZmKAkBmvhd4b9tzSNI8mrWnjyRJLTIKkqRFRkGStMgoSJIWGQVJ0iKjIElaZBQkSYuMgiRpkVGQJC0yCpKkRUZBkrTIKEiSFhkFSdIioyBJWmQUJEmLjIIkaZFRkCQtMgqSpEVGQZK0KDKz7RnWLCIWgK8u2XwccOMyV19u+9JtS89vA3ZOOeZKDjVv1bdd6brjLp/ksZtkWxOP56HmqPq2k1xvNT+Lh9ruz+hkl1f1Mwr9eUzHXe87M/P4ZS/JzF79A86fdPvSbcuc39HWvFXfdqXrjrt8ksdukm1NPJ5NPaaTXG81P4uTPqb+jE5+2Vq39eUxXev36OPTR3+/iu1Ltx3qtnWa5nuu5rYrXXfc5ZM8dqvZVrcmHtNJrrean8VDbfdndLLL/Rmt6Ht0+umjukXEjszc3vYcfeHjWT0f0+rN+2Paxz2FKp3f9gA94+NZPR/T6s31Y+qegiRpkXsKkqRFRkGStMgoSJIWGQVJ0iKjsEYRcY+IeENE/E3bs3RVRGyMiAsj4nUR8bNtz9MH/lxWKyJ+rPz5fFtEPLrteZowl1GIiDdGxCAiPr9k+2Mi4ksRcWVEvGDcfWTmlzPznHon7Z5VPrY/DvxNZj4DeELjw3bEah5Tfy5XtsrH8+/Kn89nAT/dxrxNm8soABcAjxndEBHrgD8HHgucDjw5Ik6PiPtFxD8s+Xfn5kfujAuY8LEFTgKuLq+2v8EZu+YCJn9MtbILWP3j+Tvl5b23vu0B2pCZH4uIU5ZsPhO4MjO/DBARbwXOysyXAI9veMTOWs1jC1xDEYbPML9/oKxolY/pFQ2P1zmreTwj4gvAS4H3ZealzU7aDv+PeNCJHPyrFYpfWCce6soRsTUi/gJ4QES8sO7hOu5Qj+27gCdFxGtp5/gzXbbsY+rP5Zod6mf0V4FHAT8REc9qY7CmzeWeQhUycxfF84xao8zcAzyt7Tn6xJ/LamXmq4FXtz1Hk9xTOOha4OSR8yeV2zQ9H9vq+ZhWy8ezZBQO+jRwWkTcPSI2AGcDF7c8U1/42FbPx7RaPp6luYxCRFwEfAK4d0RcExHnZOY+4FeADwBfAN6emZe3OWcX+dhWz8e0Wj6e43mUVEnSorncU5AkLc8oSJIWGQVJ0iKjIElaZBQkSYuMgiRpkVHQmkTEzTXc5ykR8TNjLrs1Ij4z8m9DRPxCRPxZRPz2yPb9I6d/LSLuHREfKc9/ISLOX+H+r4iIv4iIwyLiCSsdRn3Mf89HImL7km2/HxEvWbLtjPLAa4e6nz+IiOeuZYZl7uuoiPhoRKwr/5szIl40cvm2iNgbEX+2mu8dEWdHxG9POMOGiPhYRHiYnRlkFDRLTgGWjULpvzPzjJF/dwwvyMwXD7cDt45cZ3jsmleV578beM24+we+h+LwyT+WmRdn5ksr+G8buohvPy7/2eX2JjwdeFdmDg9V/hXgR0Yu/0lgLYu2Hgu8f6UrRcT68n+3DzEnn0/QNUZBU4mIh5V/Ef9NRHwxIt4cEVFedlVEvCwiLouIT0XEPcvtF0TET4zcx3Cv46XAD5R/rf9GhWPeheKolwBk5mXjrlyubv034J7DPZFyzndHxFPK08+MiDeXpx8dEZ+IiEsj4h0RccyY+/5P4JsR8X0jm38KuCginhERn46Iz0bEOyPi6KW3H937KP+qv6o8vS4i/ri8/eci4pmHGOFngXePnL8F+MLIHs1PA29f5vueGhGXjpw/bXi+/N/7DODSiPjBkb20/4iITeXPyL9ExMUcPLT335WzaMYYBVXhAcB5FH9d3wN4yMhlN2bm/YA/A/5khft5AfAv5V/0r1rm8lNHfuGs5gNPXgX8c0S8LyJ+IyI2j7ty+cv4kcDSeJwL/F5E/ADwHOBXI2IbxQewPCozHwjsAH5zhXkuotg7ICIeDFyfmf9F8Rf892bm/SkOtbCaT1A7h+Kx/l7ge4FnRMTdl/x3bQDukZlXLbntW4GzI+Jkig87+trSO8/M/wZujIgzyk1PA/6qPP0A4LNZHB7hucCzyz2uHwBuLa/zQODXM/Ne5fnPl3NqxhgFVeFTmXlNZh6g+MCcU0Yuu2jk6/dP+X1Gnz569qQ3ysy/Ar4beAfwMOCTEXHEMlc9NSI+A3wceE9mvm/J/XwD+D3gw8BzMvN64MEUMfx4edunAt+5wkhvozg+/2F861NH9y3/or6M4q/o+0z63wg8GnhKOcO/A1uB05ZcZxtwwzK3fT/wQ+UsbxvzPV4PPC2KTyn7aeAt5fbHAMPH6uPAKyPi14DN5V4XFD8jXxneUfn01R0RsWmi/zo1xhd6VIXbR07v51t/rnKZ0/so/yApfzFuqHU6IDO/BrwReGMUn817X+CSJVcbvqYwzv2AXcBdy/MBfDAzn7yKWa6OiK8APwg8iYOxvIDidYzPRsQvUARsqcXHDjhyZHsAv5qZHxjzrW9dcpvhPHdExCUUez+nc+jPy34n8PvAPwOXlJ/dAEWQnlTe10sj4j3A4yhC+cPldfYsc39HALeNmVctcE9Bdfvpka+fKE9fBTyoPP0E4PDy9G6g8r8co/hA9sPL099B8Vf0qo+VHxFnUryg+gDgueXTM58EHjLyesnGiLjXmLsZuojiaa0vZ+bw9Y5NwHXlrId6vv0qDj52PzGy/QPAL438d94rIjaO3jAzvwmsi4hvCwPwCuD55d7PsjLztvL7vJbyqaOIOA5YPwxERJyamZdl5v+lOBz1dy13XxGxFdiZmXsP9f3UDqOgut0pIj4H/DowfPH4dcAPRsRnKf5KHv4V+Tlgf/lCa5UvND8a+Hz5/T4APC8zv76aOyifbnod8PRyr+M5FHseO4FfoHih+HMU4Vv2F+ES76B4emj0XUe/S/HUz8eBLx7idi+n+OX/HxRPBw29nuJF3EvLPaG/ZPlnAv4ReOjSjZl5eWZeuMz11/Ote4JvBg6U9wPF007/NHL5eRHx+fKx2MvBp5WWejjwnkNcphZ56GzVpnxnzPbM3Nn2LCpExAOB38jMn5/w+n8LvC4z31uefy5wXGb+bnn+9cDrM/OTq5zjXcALyndjaYb4moI0RzLz0oj4cESsG1mrsKzyBe//pNwrKANxKvCIkfv7xdXOUL4L6u8MwmxyT0GStMjXFCRJi4yCJGmRUZAkLTIKkqRFRkGStOj/A4yHR5uTZYb3AAAAAElFTkSuQmCC\n", "text/plain": [ "