read data_plot mass density.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  view the infomation of fits file we use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filename: /root/cosmos/galaxies-mzr/data/manga_firefly-v2_4_3-STELLARPOP.fits\n",
      "No.    Name      Ver    Type      Cards   Dimensions   Format\n",
      "  0  PRIMARY       1 PrimaryHDU       4   ()      \n",
      "  1  GALAXY_INFO    1 BinTableHDU     31   4675R x 11C   [20A, 20A, I, 20A, E, E, E, E, 20A, 20A, 20A]   \n",
      "  2  GLOBAL_PARAMETERS    1 BinTableHDU     41   4675R x 16C   [E, E, E, E, E, E, E, E, E, E, E, E, E, E, E, E]   \n",
      "  3  GRADIENT_PARAMETERS    1 BinTableHDU     41   4675R x 16C   [E, E, E, E, E, E, E, E, E, E, E, E, E, E, E, E]   \n",
      "  4  SPATIAL_INFO    1 ImageHDU         9   (4, 2800, 4675)   float32   \n",
      "  5  SPATIAL_BINID    1 ImageHDU         9   (76, 76, 4675)   int16   \n",
      "  6  LW_AGE_VORONOI    1 ImageHDU         9   (2, 2800, 4675)   float32   \n",
      "  7  MW_AGE_VORONOI    1 ImageHDU         9   (2, 2800, 4675)   float32   \n",
      "  8  LW_Z_VORONOI    1 ImageHDU         9   (2, 2800, 4675)   float32   \n",
      "  9  MW_Z_VORONOI    1 ImageHDU         9   (2, 2800, 4675)   float32   \n",
      " 10  E(B_V)_VORONOI    1 ImageHDU         8   (2800, 4675)   float32   \n",
      " 11  STELLAR_MASS_VORONOI    1 ImageHDU         9   (2, 2800, 4675)   float32   \n",
      " 12  STELLAR_MASS_REMNANT    1 ImageHDU         9   (4, 2800, 4675)   float32   \n",
      " 13  SURFACE_MASS_DENSITY_VORONOI    1 ImageHDU         9   (2, 2800, 4675)   float32   \n",
      " 14  STAR_FORMATION_HISTORY    1 ImageHDU        10   (3, 8, 2800, 4675)   float32   \n",
      " 15  SIGNAL_NOISE_VORONOI    1 ImageHDU         8   (2800, 4675)   float32   \n"
     ]
    }
   ],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "vac0 = fits.open('/root/cosmos/galaxies-mzr/data/manga_firefly-v2_4_3-STELLARPOP.fits') \n",
    "vac0.info()    # 查看fits文件信息"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ##  1. read data from the file    \n",
    "> ##  2. calculate the metallicity using the formula from the paper   \n",
    "> ##  3. use Marvin to calculate the mask   \n",
    "> ##  4. plot the unmasked surface mass density   \n",
    "> ##  5. plot the masked surface mass density   \n",
    "> ## PS. only need to change the parameter of the plate_ifu and prop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Could not memory map array with mode='readonly', falling back to mode='denywrite', which means that the array will be read-only [astropy.io.fits.file]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: No release version set. Setting default to DR15\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/root/anaconda3/lib/python3.7/site-packages/marvin/tools/quantities/map.py:400: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  return 1. / ((1. / ivar1 + 1. / ivar2))\n",
      "/root/anaconda3/lib/python3.7/site-packages/marvin/tools/quantities/map.py:424: RuntimeWarning: divide by zero encountered in power\n",
      "  return np.log10(np.e) * ivar**-0.5 / value\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 850x1000 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "import copy\n",
    "\n",
    "#read in the vac file\n",
    "vac0 = fits.open('/root/cosmos/galaxies-mzr/data/manga_firefly-v2_4_3-STELLARPOP.fits') \n",
    "vac = copy.copy(vac0)    # 去除只读属性\n",
    "\n",
    "#plate and ifu ID of the galaxy\n",
    "# such as 9883-3701\n",
    "plate = 9883\n",
    "ifu = 3701\n",
    "\n",
    "#choose the quantity to make the map, could be changed to other properties for which the extension names are listed in the data model  \n",
    "prop = 'SURFACE_MASS_DENSITY_VORONOI'    # property 特性，根据上面的info来填写即可\n",
    "#SURFACE_MASS_DENSITY_VORONOI\n",
    "\n",
    "binid = vac[5].data\n",
    "basic = vac[1].data\n",
    "galid = (basic['plate']==plate)&(basic['ifudsgn']==str(ifu))    # 星系编号 galaxy ID\n",
    "mass = vac[prop].data[galid,:,0][0]    # 获取质量密度数据\n",
    "bin1d = vac[4].data[galid,:,0][0]    # 获取一维的bin数据\n",
    "\n",
    "# 导入Marvin的mask\n",
    "import numpy as np\n",
    "from matplotlib.colors import ListedColormap\n",
    "from marvin.tools import Maps\n",
    "import marvin.utils.plot.map as mapplot\n",
    "\n",
    "# 获取数据\n",
    "new_maps = Maps(''.join([str(plate),'-',str(ifu)]))    # 通过Marvin导入星系数据\n",
    "ha = new_maps.emline_gflux_ha_6564    # 获取Hα\n",
    "nii = new_maps.emline_gflux_nii_6585    # 获取Nii\n",
    "oii7 = new_maps.emline_gflux_oii_3727    # 获取O\n",
    "oii9 = new_maps.emline_gflux_oii_3729    # 获取O\n",
    "\n",
    "# 计算金属丰度\n",
    "RR = (oii7+oii9) / nii    # Oii/Nii\n",
    "# 论文中的关系式【LJ Kewley, MA Dopita 2002】【https://arxiv.org/pdf/astro-ph/0206495】\n",
    "logRR = np.log10(1.26602*RR+0.167977*RR*RR+1.54020)\n",
    "oh =  logRR + 8.93\n",
    "# 计算mask\n",
    "masks_bpt, __, __ = new_maps.get_bpt()    # 使用bpt图筛选star-forming区域\n",
    "masks_bpt['sf']['global']\n",
    "#RR.pixmask.schema    # 查看mask的种类列表\n",
    "mask_non_sf = ~masks_bpt['sf']['global'] * RR.pixmask.labels_to_value('DONOTUSE')\n",
    "mask_bad_data = RR.pixmask.get_mask(['NOCOV', 'UNRELIABLE', 'DONOTUSE'])\n",
    "min_snr = 3.\n",
    "mask_nii_low_snr = (np.abs(nii.value * np.sqrt(nii.ivar)) < min_snr)\n",
    "mask_oii7_low_snr = (np.abs(ha.value * np.sqrt(oii7.ivar)) < min_snr)\n",
    "mask_oii9_low_snr = (np.abs(ha.value * np.sqrt(oii9.ivar)) < min_snr)\n",
    "new_mask = mask_non_sf | mask_bad_data | mask_nii_low_snr | mask_oii7_low_snr | mask_oii9_low_snr\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAERCAYAAADYEnSgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXmcHFW5v5/vLEkIBEJIgJAQFlllhwHZDZsCIrggi8h+zQ+uIHivInBVFFHRqyJe1BhklU1AgYiILIKgAhog7FvYAwESIAuEJLO8vz9ONXT61EyqZrp7umbeh099pvutU6dODZO3T7+rzAzHcRynPjT19wIcx3EGE650Hcdx6ogrXcdxnDriStdxHKeOuNJ1HMepI650Hcdx6ogrXcdxnDriStdxHKeOuNJ1HKcwSPqKpMckPSrpSknDKs4fJWm2pOnJ8R/9tdbucKXrOE4hkDQO+DLQZmabAs3AISlDf2dmWybHb3pxn0l9XGqPuNJ1HKdItADLSWoBhgOv1uAekyXdJGmNGsxNSy0mrQejR4+2tddeu7+X4TgDlvvvv3+OmY3pyxwf3215e/Otzmz3e3jxY8CiMtEUM5tSemNmr0j6MfAS8B5wi5ndkjLVZyXtCjwNfMXMXs657O8ApwOPSDrRzK7IeX2PqKgFb9ra2mzatGn9vQzHGbBIut/M2voyR9sWw+xff5mQaWzz2Gd6vJ+klYHfAwcDc4FrgGvN7LKyMasA75jZYknHAQeZ2e551y1pC+ASYDPgOuA4M5uTd5403LzgOE7NMKAr438Z2BN43sxmm1k78Adgx6XuZ/ammS1O3p4PbNOrdZs9BLQB3wf2Bx6V9BlJa6QdeeYurHnBcZzGxzDaLZt5IQMvAdtLGk4wL+wBLPV1V9JYM5uVvN0feKK3NzOzDuCbkq4DbibsrFOHkkOXutJ1HKemZNzFLhMzu0/StcADQAfwIDBF0pnANDObCnxZ0v7J+beAo/pyz8RccSowmqDAH+jLfOBK13GcGmIYnVX0G5nZGcAZFeJvlZ0/DTitGveSdAAwGRgD/Aj4lpkt6eu8rnQdx6kpXRTLWS9pJeD/gMOAZ4HPmNk91Zrfla7jODXDgM6CKV3gMWAs8EvgFDN7r5qTu9J1HKemFG2nC3QCe5rZHbWY3JWu4zg1w4D24uUCPFvLyT1O13GcmmEYnRmPBmIiwbxQE3yn6zhO7TDobCh92v+40nUcp2aEjDSnHFe6juPUENGJ+nsRvaFm+3NXuo7j1IzgSCuk0r1U0iUZx5qZDc06cb8oXUnNhJzpV8xsP0nrAFcBowhpdodXI/PDcZz+JcTpFlLpPgvMrsXE/bXTPYmQx7xi8v6HwDlmdpWkycCxwK/6aW2O41SRrmLudM+sdh3dEnUPGZM0HvgE8JvkvYDdgWuTIZcAn6r3uhzHqT6lnW6WY7DQHzvdnwGnACOS96sAc5MyagAzgXFpFya9iyYBTJiQrTCy4zj9hyE6PR1gKer625C0H/CGmd1fLk4Zmuo5NLMpZtZmZm1jxvSpi4jjOHWiy5TpGCzUe6e7E7C/pH2BYQSb7s+AkZJakt3ueGrTbM5xnDpjiCXW3N/LyMsXgXtrNXldd7pmdpqZjTeztQmtk/9qZocBdwAHJsOOBG6o57ocx6kNITmiKdPRKJjZBWb2XK3mb5Q43a8DV0k6i1AN/oJ+Xo/jOFWiaE4ySXm7Q5iZZe7F1m9K18zuBO5MXj8HbNdfa3EcpzaYiU5rnF1sRlbFM9IcxykqXQXb6ZrZ+FrO70rX6ZarZrRFskPWm5YysrG58OmdItkxG/yjH1Yy+AiONFcz5fhvw3GcmlFypA0UkhIG2xByCR43s6fyzjFwfhuO4zQknaZMR6MgaVdJP5W0eoV8LeDfwD2EDNrHJU3JO78rXcdxakYpIy3L0UAcDRxoZq9VyC8GtgTuI3QLfgo4VtIReSZ384LjODWlq3jRC9sBU8sFkjYEPgr8HfiomZmk4YQQ1yOBS7NO7kp3gHPmo5+MZBOGzIlkR21wTyRbvWVepns8+fIakWxUyr+zVcd5ouFgIxS8KZzSXR14pkI2kfA455uFTptmtlDSlcCX8kzuStdxnJphiPbipQEPAxZWyLZNfv6tQv4isFKeyV3pOo5TM8woYnLEy8CHK2Q7E4p1vVQhXw7I9pUwoXC/DcdxioToyng0EHcDR0jaGEDS/sAGwF9Sxm4GvJJncle6juPUDCPsdLMcWZD0FUmPSXpU0pWShlWcHyrpd5JmSLpP0tq9WPbZBBPDI5JeBf4ALAF+UnGvJmB/IFemjZsXBhCH3DMpkn14RDyuWXFa+RXPxKUv1h8St6m7ZkZc12PP4fE/mJkdkYjWV+PsynldnZHskSWjI1lzSip8q+KbjGhaFMnWbB0WydKy1B55N15fmj3yvK1r0sVlwFItR5qkccCXgQ+b2XuSriZUK7y4bNixwNtmtp6kQwitwA7Ocx8ze1bS7sC3gfWA24DvmtkjFUP3AOaTsyqiK13HcWqGUfUC5S3AcpLageHEtbcPIChLCAkM50lSKeIgK2Z2H7DPMsbcCmycZ15wpes4Tg0JLdgzq5nRksqLe0wxs/czvszsFUk/Bl4C3gNuMbNbKuYYR3CEYWYdkuYRWoLFcZL9hCtdx3FqSK6mk3PMLK6yVJpJWpmwk10HmAtcI+kLZnbZUjeMybXLlbR5nvEAZvZw1rF1VbqJ0fsuYGhy72vN7AxJFxOyPUqhF0eZ2fR6rs1xnOpjVDUjbU/geTObDSDpD8COQLnSnQmsCcyU1EKIoX0r532mk09RGzl0ab13uouB3c3sHUmtwN8l/Tk59zUzu7aHa51lMHroO7GsZUGma9dofTuSvdAeO7RGNr8bydqtK5Kt09Iayd7qih1zC7qyBc5vMeTNSPbQklUyXTu7Y8VIluZwG9oUO+YWd8TP8eUHD13q/c+3ujLTOgYrVewc8RKwfZJ++x7BkVVZa3QqIS33HkILsL/mtecC3yeb0v0YHyRNZKauSjd5+JJmaE2OmlVodxynfzFT1Xa6ZnafpGuBB4AOQt2DKZLOBKaZ2VRCq6/fSppB2OEe0ov7fKOn85K2JkRFbAu8Dfwgz/x1t+km9SjvJ4Ri/CL5RR4PfE/St4DbgVPNbHHKtZOASQATJkyo46odx+kNwZFWvTRgMzsDOKNC/K2y84uAz1XthmVIWhf4XjL/EuB/gR+Y2dw889Q9OcLMOs1sS0Kr9e0kbQqcBmxE+OQYRWhUmXbtFDNrM7O2MWPG1G3NjuP0FlU1OaI/kDRa0s+BxwkK91JgAzP7el6FC/2YkZYs9k5gbzObZYHFwEV4k0rHGRAER5oyHY2GpOGSvgk8C5wA3AJsbmbHmNnM3s5b7+iFMUC7mc2VtBzBG/lDSWPNbJYkAZ8CHq3nuorI0f86OpItn/J/c0FXnI01rKk9kg0hzgzbcVhcirFZ8T+O9hQ/xTuk3CPl2g1ah0Sy4YpDKuPVwSKLnVw/n7VXJNtz1OMp82Xbb8xtXy6SXbH9+ZmudQJFK+2YmEAnEcwWqwH3AqeY2d+rMX+9bbpjgUuSh2oCrjazGyX9NVHIIoRrHFfndTmOUwNqkJFWDx4n+JyeAo43s+urOXm9oxceBrZKke9ez3U4jlM/CtiYcn2CZWQN4EJJFy5jvJlZtvhFPCPNcZwaYgbtXYVTuv+khqGsrnQdx6kZwbxQLKVrZjvXcn5XugVl9uIVItkGy78eyTYcOiuSLeiKnUNrtLwXyca2xPd4eEmcybVuyl/RA0tiB96KikKvoSWeb1iKCTDNgbd+6+xI9v3xN0ayZ9rjbipPLo77um0yPK5FPbo1zvJz8lHFjLQBQbE+ghzHKRRFDBmTtJ2kURnHrpO3BbsrXcdxakgwL2Q5Goh7gL1LbyStLOktSbukjN2RkFuQGTcvOI5TUxqs/1kWKhfcBIwk1IrpM650HcepGSF6oXAt2GuKK92CsuGI2Gk2a0nsMNp4WOwcWqU5dg4tymhT23xI7CBrtzhfbNd4GPO64hKQ7SmBOW+mhBiNbIqvfTelI0FnV/wcazTH5S3bh7wRyV5I6c02rzN2Oh75r2OWen/JdssK4xy8FDQ5oqa40nUcp6YU0LxQU1zpOo5TM0rRCwOEqiRMuNJ1HKemNFhkQlZOlnRg8noIQeGeKakyOHx83old6RaUue3DI9m2Kz6f6drdlottugtikylPt8eteYYp/rDvTPn8X74p3t2s2rx8JHtsSZyUsUZzPGF7ylfUtVvi9j8rNcVVy17siG3OwxRXQRvZvDCSzWofGckmLLd0y6397z4hGjN1l/Mi2WDETHQUU+m2JUc5O3UzNtcO2JWu4zg1pYDmhaqEhnWHK13HcWpGEW26ZinhOFXEla7jODWlaEo3DUlDgRWB+Wn9G/NQSGOL4zjFoBSnW6TaCyUkrSXpl5JeBBYCrwELJb0g6TxJveqOW+92PcOAu4Chyb2vNbMzJK0DXEVoSvkAcLiZxV4S531WTKkKlvaHm5YI8WJH/KtdrTn+/N2gNU4MmNMZO9faU/wIaU6zWR3xWkalfOynOc3S5uu02PvXrHjCVuKNyZjm+DmeXDw25R7L3pc8/3am2iiDliLG6Ur6KHADYXfbDjwJzE/erwf8J/B5SZ80s3/kmbveO93FwO5mtgWwJbC3pO0JPeTPMbP1CX3kj63zuhzHqQFm0NHVlOloFCStCFxJqMFwPLCSmW1iZjuY2SbASgSl2wxcJWlEnvnr+qRJx9/Sdqc1OQzYHbg2kV9CaE7pOM4AoIDmhcOB1YEDzOzXZrZU0WczW2Rmkwl6ahzwhTyT1/3jRVKzpOnAG8CthPbGc82sIxkyk/AgaddOkjRN0rTZs+MC1o7jNBYFtenuA9xpZnf2NMjM7gDuBD6RZ/K6K10z6zSzLQmZHNsBG6cN6+baKWbWZmZtY8aMqeUyHcepEmbKdDQQmwF3ZBx7RzI+M/0WMmZmcyXdCWwPjJTUkux2xwOv9te6+pu1L/tBJDtks2mRbGjKx+VbHXF7nVc6Vo5ka6RkXs1oj8vvjWyKHV8TWmLnWquyle5La//TF9KcZmms0TI0ks3ujJ1rWw57KZKlZalV8kc2zbSOwUoBHWmjyK6DXgEydwKGOu90JY2RNDJ5vRywJ/AE4dOilOd8JMFr6DhOwTGrnk1X0oaSppcd8yWdXDFmoqR5ZWO+1YtlLw/E4UHpLAbinUgP1HunOxa4RFIzQeFfbWY3Snqc4AU8C3gQuKDO63IcpyaIzipFJpjZU4SoJxId8gpwXcrQu81sv77ero/Xd0tdla6ZPQxslSJ/jmDfdRxngFEje+0ewLNm9mItJge+L+lrGcbF9rtl4GnAjuPUjJy1F0ZLKndgTDGzKd2MPYQQS5vGDpIeIthlv2pmj2VdQMKrBN2Y1VufywflSrfBWH9c3Ebmnc7YEZTGqkPmR7LmlG9J87riIkqtirO7PtQaO77eSMlIS8sWy0paq5+sjrmsDFVa0ajYkTZMHZHs5fZlZ5utMKxPqfgDGwt23YzMMbPKcooRkoYA+wOnpZx+AFjLzN6RtC9wPbB+5hUAZpa7Rm4eGicNxHGcAUkXynTkYB/gATOLGgWa2fxSApaZ3QS0Soqb39UIScMlnd5TXQZXuo7j1AxLHGlZjhwcSjemBUmrS1LyejuCjnuzzw+SnRHAdwn1GVJx84LjODUlh3lhmUgaDuwF/L8y2XHhPjaZEHp6vKQOQtjXIWbVXEG2ZfZ00pWu4zg1pZrRC2a2kIpkhETZll6fBzR0ryRXug3GrbudE8l++sRema5Nc5qllXbsTPkgHt0U9wxLoy9OszSq7TTLyviU7LhXO+PfQatiR9/bHUv/DlYauiga4wTMahYyVlhc6TqOU1MarJhNv+NK13GcmlJ3i2qD40rXcZyaYYiuBipQ3gi40nUcp6b4RndpXOkWgOfey5aNuOnyr0Sy4U1xttSoprhHWv+4sxqLhV1x5t87ncP6YSUDiMHnSOskFOLp1rvqStdxnNoywLa6klYldLeZYWYLys+Z2RxgzZ6ud2OL4zg1pWidIyRtJunLlenDkkZJ+iMwC5gGzJZ0et75e9zpSvprjrnMzPbIuwDHcQYuBnR1NY5Czch/Ap8FflEh/w2hH9oLwMPAjsB3JT1iZn/MOvmyzAtNZP9yULjfrOM4NcaABtrFZmQHYKrZByXwkgI2BwCPAB8xs0WSVgPuJ6QkV0fpmtnE3qy4OyStCVxKaG/cRaiXea6kbwNfBEotfk9PKgQ5QJOyfe7N64y7hrzWsVIk23C5uP5Ha4or7dn2OJstrdzjQOHDQxZEsnfthUj2l7lL9yFcb4R3pu6JAsbpjgN+WyHbg7Cx/EWpJbuZvS7pt4QWY5mptyOtA/hvM3tA0gjgfkm3JufOMbMf13k9juPUmuIp3RWAuRWybQlPUtkleAa1bEwpaZykn0qaJul5SZsm8pMlfWRZ15vZLDN7IHm9gNCUclyeNTiOUySyOdEayZFGCPn6UIVsR2CumT1TIW8B4sr+PZBZ6UrahGDPOJzQnmICMCQ5vRZwUp4bS1qb0C/tvkR0gqSHJV0oKbXvkKRJicKfNnu2f6VznEJgGY/GYRpwRGKzJdlQbg6kBRZsTM52PXl2uj8h7EzXAT7D0o6zfwLbZ51I0grA74GTzWw+8CvCJ8uWhHCMn6RdZ2ZTzKzNzNrGjMnavshxnH7DwLqU6WggfkjwOz0p6Z/A7QQf1LkpY/fjg41jJvLYdHcGDk16D1V6XV5PFrlMJLUSFO7lZvYHCAbpsvPnAzfmWNeAp5m4f1kaKzW/F8nSSjs+uDguz7hpihNpIDvN0ni3K95uLeoaEsmGNsW91JyeaCiFukzM7EFJBwI/ANqA54AzzOzv5eMk7U0wj96cZ/48Srenf/mjCVXaeyRpo3EB8ISZ/bRMPtbMZiVvPw08mmNdjuM0Mo1lOsiEmV1PaGrZ05ibgThkaBnkUbr/Ao4mPR7tIOAfGebYiWATfkTS9ER2OnCopC0J/3teoKwVh+M4BaeASreW5FG63wVuk3QLcAXhV7mnpJMIu9NdlzVBsj1P+67hMbmOMxApYHKEpBZgOPCOmXVVnPsscCzBrPAYcJaZPZ5n/syONDP7G/ApgiPtQoLyPBvYBfiUmeUyJjuOMzgILXuWfTQQ3wDmEMym75NsMK8G9gY2Aw4B7k4isTKTKznCzP4E/EnSesCqwJtm9lSeOZz8vL54xUzjtl3h+Uj2bkq5wtVbYqfZwgb7q68m7Rb3OUvrzTavqzWSpTkiV2heulzmdza7oQ+rGwQ0VmRCFnYBbjezN0oCScsB3yHE5B4E3Ev4hj8Z+BrwpayTZ1a6kppKW20zm0HIxCg/v7KZvZ11PsdxBgcZs9gbiQ0IDv9ydgdWBM5OHGgAF0naE9gzz+R54nTvToo+REiaCDyU58aO4wwCsiZGNJZiHg28VCH7CGGVf6qQ3wuMzzN5HqU7GnhI0sElgaRmSd8DbiMkTjiO45Sh4EjLcjQO84CRFbKPELpC3F8hj+1PyyCP0t0KuAG4MknV3ZyQifZV4FQz+3jemzuOMwgo3k53BrB/6Y2klQi1Fx4oVRgrYzzwBjnIbNM1s4XAUZJuBi4ilDN7FtihVMTGqQ1XbH9+JPvv6QdFsjSn2ZudcVbZyOaFkayVgZtltdjaM8k6M+5BWhU75pweyJZQ2UicT7DXXkWot3AQIYTsspSxu5LzW36u6AVJw4GPA0OBNwklzdYGXOk6jhNTwDhdgnLdn6BsS7ubPxMiFd5H0jrAbsApeSbPU2VsG+BBQqzuIYR43T8B10j6dRJS4TiOsxSybMcy55E2lDS97Jgv6eSKMZL0c0kzkqqFW+ddr5l1mtlnCUW8Dgd2NLNPmFnl18FW4DDSd8Ddkmen+09CybOPmdmLieyIxNzwS0Js24fz3NxxnEFAley1SU7AlhCc+IS6t9dVDNsHWD85PkKoYLjMWt/d3O9fhPIH3Z1/GnhaUqZiXyXyONJ+BOxapnBLN76C4GSrrLTuOI5TK/YAnq3UR4Q+Zpda4F5gpKSx1b55Erl1gKSpQOUaeiSPI+2bPZx7XtIueW7s9I00Z86CrmGRbERTpbMVltfAcJot7FoSyYY3xaUYF6VkpLWnZOB1ks0R+eS7q2VdokOu5IjRkqaVvZ9iZlO6GXsIcGWKfBzwctn7mYlsVsrY3EjaCDiGYHZYleAm/FueOarWI628c6bjOA6Q9GDP7EibY2ZtyxokaQjB0XVa2uluVtFrJC1PUPLH8EGzhnuA/wFuMLO402sP5I1e+DhwHLAhEG2rzGzdPPM5jjMIqH4M7j6EmNnXU87NBNYsez+enO10SkjaiVBR7EBCs8pHgJ8C/wX8tNSEIS95ohf2JZRgHA5sBDxJSJVbk15ssR3HGRxUK3qhjENJNy0ATCU4+CVpe2BeWYOEbOuVTpH0JHA38ElCXkKbmW1BCBvrUwxcHkfaN4FfAPsm779hZhOBTYBmQhxbj0haU9Idkp6Q9FhSKg1JoyTdKumZ5GdqY0rHcQpIFTPSklyBvYA/lMmOk3Rc8vYmQnudGYQkh//sxYrPJujGzwBrmNlJ1UwAy2Ne2Aj4FmFXa6VrzexpSd8mKOWrlzFHB/DfZvaApBHA/ZJuBY4ilFI7W9KpwKnA1/M8SFHZ5qb/Wer9/ft+L9N1aX26nl4YR658bGTc+eiFjvgzranlrUjWF5fvnM64K/Xo5rg3W19Ic5qlMSyljOOCrvj3NyzFOdmZsqkZP8wDdXJRRfNCkhm7SoVsctlrI0eZxW54A1iP0CNtQ0mXmdkrfZzzffLsdLuAjuShZhNasJd4lbhPfISZzSp9YpjZAkL63DhCmMclybBLCAkYjuMUnKymhQYr/ziOsMudAZwFvCDpFkmfJ5hX+0Sene5ThJRfCEkSJ0v6B8nuldDbLDNJtfWtCO2LVyvZXcxslqRVu7lmEjAJYMKE1CqTjuM0GgUrYp5EYl0PXJ8kPhxN+DZ+GbCYsHdfR5KSTWgu8ux0Lwc2Tl6fQbDlzgReIxT4/VbWiSStQGjDfrKZzc96nZlNMbM2M2sbM2ZM5oU7jtN/FHCn+z5m9pqZ/cDMNgQmEkyo7xGSxV6TNFnSx/LMmSc54hdlr++XtBkhdGM54LaszdkktRIU7uVlIRevl9qwJ9kjuUqlOY7TwDSoQs2Lmd0F3CXpBEIExbGEb95fJAQTZKLXyRFmNpPgHcyMJBHaYDxhZj8tOzWVUCry7ORnoZpOVTrD8pDVcVZJWl+uHzy2byR7syPOqGpNyUhbuyXutPTEkrgE5HqtcdZWWr+xvjjNnm+P60Kv0xo/RxrvdMUZeCs0xZl6K6R8x0sr9/hCeywb3hxnwjnd0MC72N6S+KOmAFMkbUpImshMbqUraTdgB4Kx+RXgn2Z2Z8bLdyKkzz0iaXoiO52gbK+WdCwh9vdzedflOE6DMsCUbjlm9ighWSIzeRpTjgKuIdg1DHgbWDmc0p3A58wsjjtaeoF/p/vA4j2yrsVxnOKgghUxT6IU8mBm1l2yRkSene7PgW0JO9VrzKw9sc8eRCjteG5yznEcp8hcxgf7855CLyw5b3SfIReRR+l+EjgtKeUY7mjWDlye7ILPyjGX4ziDhWKaFxYRst7+TGhIWTXyKN1O4Jluzj1FlRdWJHrrDKs2TSnf45ZvWhzJhjXFzqGXOyqbn8IwxePWbY3n60xpx9KsPNGIS5PVaZZGmtMsK+0phfKaU36n8zq8SUpmiulI+zohLvcwQjjspcBFSdHyPpPnX8YNwMHdnDuEEEzsOI6zNAXrBmxm/2tmmxAc/38CjgeekPR3ScckeQa9Jo/S/SOwp6Q/STpK0j7Jz5sITrCpknYvHX1ZlOM4A4iCKd0SZnaPmU0ilCE5GmgnhMnOknShpF61J8tjXrg2+bkmISmikt8nP0uG5czBwo7jDExE8aIXKjGz9wgmhkuTDsCTCfkELwBn5p0vj9LdLe/kjuMMcopp041IajAcTtjxbkRo/zO9x4u6IY/SvRtoKm9DnHSS2BT4q5k92JsFONWjy2Jr0WsdK0Wy5pTvcusPfS2SrdUSl8WYl7F4yZKUOiDjW/pkCqs5TSnWto1a50Syl4ctnaX+xWlHZJp/9qL4+a/f+ZcZV1dgCqp0JbUQoraOAT5OeJIbga8BfzazXu3h8yjdKwkVdo5IFnQcIT4XoF3SJ8zstt4swnGcAUzBlK6kzQk72sOA0YQ2PacAl5lZ/CmckzxKd3uWLiz+NeA3hLKOUwhN2lzpOo6zFAU0L0wnVBK7jtCq59+lE5JWTLsgT7XEPEp3VUKtBSStB6wDnGdmCyRdBFzR08WO4wxSiqd0IVRPPDQ5lsX7nXSykEfpzueDNhkTCe2SH07ed5LSHdhxnEGOFTJ64fJaTp5H6f4TOFVSB3AyoQFcifUIBc2dOrHf3SdGsrHLjY5kW63wUiRbc8ibkSyt3ONbXXEPsoUpWVvjU/6KhqtY3QIAmlPW/K7FD7d667ylx3TE5S7nLKpuP7hCU7CdrpnVtIZMnuSIU4BRhNq3w4Bvl507GLinestyHGegUOTOEXmRNDLpp7ZFd2PydI54BthA0ipmVrlVOonQtsdxHGdpBohCzchQYE8qOhaXk7uIeYrCxcweyTuP4ziDgAZN8e1Pel8Kqhck+cpvSHq0TPZtSa9Imp4ccc8Zx3EKiRhc5oUs9LpHWi+5GDiPkMdczjlm9uM6r6XQ3LjL/0Wysx7dL5KNbI77nC2vuMfX8JQSkJ0p9ZvHNMeu6PaUfzALUzLSWlNKRbaklOjoS1nIvrCwK15fa8q4Uc1L93BbbWgcotmVUu7yqh2m9HptRWYwKdQs1PWvO+mm2WNLH8dxBhgFrTJWK/pnSxFzgqSHE/PDyt0NkjRJ0jRJ02bPnl3P9TmO01tc6S5FIyjdXwEfArYkVO75SXcDzWyKmbWZWduYMWPqtT7HcXpLRntuVhNEEpJ1raQnJT0haYeK8xMlzSvzEX2rFo/VF+pt040ws9dLryWdT6ji4zjOQKG6u9hzgZsUxrR7AAAU8UlEQVTN7EBJQ4DhKWPuNrPYwdEg9LvSlTTWzGYlbz8NPNrTeKd7vrFp/Hl1/lO7RLJ3Lc40a0qpUrdK03uR7IWO+Nr1W2IH1KKUf2nzumIH3qrNcebWvK74vis11b4v2aKU38Hsrrgc48vtS4dgfnj4q9GY+d5H7X2qlQacFJvZldC/DDNbAsR/VP3Le4Q04lndDah3yNiVhMy1DSXNlHQs8CNJj0h6mFAo/Sv1XJPjOLUlh3lhdMlnkxyTKqZaF5gNXCTpQUm/kZSWb72DpIck/VnSJrnXK22XdDjPMnYdSe8XVDaz+WZ2uJk90d01dd3pmllaxZ4L6rkGx3HqSD4n2Rwza+vhfAuwNXCimd0n6VzgVOCbZWMeANYys3eSmP/rgfVzrvoeQpeIKwAS5/6zwAFmdnfF2B0J5R8rw2C7pREcaY7jDGSqF70wE5hpZvcl768lKOEPbhV2mu8kr28CWiXFlaB6pjLIugkYSXrYdm5c6TqOUzOqmZFmZq8BL0vaMBHtATy+1P2k1aVQLk7SdgQdF5fV60f63ZHmVI/JT300kk1ofTuSjUhxkC2y+EN8kcXZYu0pstkpDrLXOmMH1Lot70SytzvjjLmhqu6fZZpjrj3FadaaUtrx1fY4bPzdrqVLOc5pH5FpHcf++6hIdsG2F2e6tsioq6rhCycClyeRC88BRyetwzCzycCBwPFJCdr3gEPMUtIj+xFXuo7j1I4qJz6Y2XSg0u47uez8eYRSAw2LK13HcWrKAKq9UJUncaXrOE5tKabSPVnSgcnrIYSnOFNSZf2B8XkndqXrOE5NKehOt43YjLFTN2NzPaEr3YKSVsaxNSVrq7M1dg6lOc1GpjjX0ko7jkwpAdlucRDMmKbYQZZWArIz47/IhSnOuuFNcXZcGotTnGZp3LsoruexJMVx+NLibpsCvM9gcJBlpnhKtyqhYd3hStdxnNpRwG7AZindV6uIK13HcWpGKU63aEjaE1ickoFWPmYXYIiZ3Z5nbk+OcByntphlOxoESROBvwBrL2PoWsAtkrqz9abiStdxnJpSwB5pxwCPmdlvexpkZpcBjwBfzDO5mxcKSloZx7SMtOYUL0ZniuNrkcV/Cp0pn8mp86U43NZojh1uaaa92Z3xtas0p5SKTNkJvdvREcmWb4rv8nJHXHJ1p2HZ9hvNKdpgfsewpd6ft/UVmeYalBSzK8TOwGUZx04FDsszuStdx3FqStEcacAawPMZxz6fjM+MK13HcWpKAZXuEiBrFfrlgPirWQ+4TddxnNphFM6RBrwAbJdx7HbJ+MzUu3PEhZLekPRomWyUpFslPZP87LYbsOM4xaOAjrRbgIMkrdfToOT8QcDNeSavt3nhYkIFoPIq66cCt5vZ2ZJOTd5/vc7rGhAct+HfItmFT8fRLCOaFkWyrpTP31bFjqo3u+LuKCOb4+yzOZ3ZknqWpNx3QXzb1JKS7SnOv9au+OIFXfE3xclz47rW8zpjh9t9b68dya7f+ZfxAp3uaSyFmoVzgEnArZKOMbM7KgckYWUXAh2EZpmZqXe7nrskrV0hPgCYmLy+BLgTV7qOMyAoYnKEmb0i6QvAVcBtkp4HpgPzgRHAVsA6wGLgIDN7Jc/8jeBIW63UDdjMZklatbuBSaO6SQATJkyo0/Icx+k1ZtUuYl4XzGyqpI8AZwF7E5pilmgHbgS+aWYP5527EZRuZsxsCjAFoK2trXj/Jx1nMFLQf6lm9ghwgKTlgA2BFQm73afMLK4QlZFGULqvSxqb7HLHAm/094IGEq+3rxTJ7nnvQ5Fs8xVmRrLVW+ZFsrREiGkL141kmy73ciRrTokdSrPLvtGxYnxtSmpFWrW0VsW1StKuXdA1LJK9snhkJNt6ZPwcZz76yaXef2vTP0ZjnA8omnmhkkTBTk87J6mV0CX42qzzNULI2FTgyOT1kcAN/bgWx3GqiQFdlu0oEJI2k3QO8CrwuzzX1nWnK+lKgtNstKSZwBnA2cDVko4FXgI+V881OY5TY4qlT7tF0orA5wm1GbYh+AkfosGjFw7t5tQe9VyH4zj1o+jmhSQ87FjgM4QMNAN+A/zQzJ7LO18j2HQdxxnAFDF6QdI44CjgaEJ42FvABcDdBHPCX3qjcMGV7oDi9Ic/kyKNHUYrtsTJEW93xEkPWRnRHM/3cnvc0iYtwWHN1rciWZfFzro5HbFD8M32eM2Lu7IlZaS1GBraFCdWuJOsjxSwypikm4C9CCv/M/A14EYza5cUe6Fz4krXcZyaEZIjCqZ1Q1zuDOBgM3uw2pM3QvSC4zgDma6MR+NwPaErxL8k/VnSIZLir4y9xHe6juPUlKLtdM3sM5JGE0JYjwauAOZLuhq4q6/z+07XcZzaYTmOBsLM5pjZT8xsU2An4BrgYEKxLgP2k7RJb+b2ne4A4vub/6HX1x53/+GRrG14XDx/33UfjWRfefDgSHbOVnG8+NceikOw0yqU/edGd0ayUx/6bCR7p3NoJDu/7dJI9sVpR/R6nNNXill7oRwzuwe4R9JJBMV7DCGy4UhJM4Dfm9npWefzna7jOLWlikXMJY2UdK2kJyU9IWmHivOS9HNJMyQ9LGnr6j2GLTSzi8xsF2Aj4MeEegy5qiK60nUcp3ZYaNeT5cjIucDNZrYRsAXwRMX5fYD1k2MS8KsqPclSmNnTZvZ1YDzw6TzXutJ1HKe2VGmnm6Th7kpIUsDMlpjZ3IphBwCXWuBeYGRSSCszkg6StGa2R7NOM5uaZ35Xuo7j1JbsjrTRkqaVHZMqZloXmA1cJOlBSb+RVJkhMw4oLw03M5Hl4Upgl9IbSStK+qekbXLOk4o70hwAJm/z215fm+Y0S+N/t7im1/c4e4vf9/rarKQ515y+o67MtoM5ZtbWw/kWYGvgRDO7T9K5hPZe3yy/Xcp1eT15lXO0AtsDcVpkL/CdruM4tcOoZnLETGCmmd2XvL+WoIQrx5SbBsYTyi82DK50HcepGcKQZTuWhZm9BrwsacNEtAfweMWwqcARSRTD9sC8UjuwRsHNC47j1JbqZqSdCFwuaQjwHHC0pOPCbWwycBOwL6F2wkJCRllD4UrXcZzaUkWla2bTgUq77+Sy8wZ8qQq3OiLZKUMo1WfACZI+lb4sOynrxA2jdCW9ACwAOoGOZRjUHScz7iDrR0o23eLxseQoJ03hQnjK4indhN3MbE5/L8JxnOqRI3qhUVinlpM3mtJ1HGdAkT3Ft1EwsxdrOX8jRS8YcIuk+1OCogGQNKkUOD179uw6L89xnNwYVa290GgkURJxtageaCSlu5OZbU3Inf6SpF0rB5jZFDNrM7O2MWPG1H+FjuPkp3hFzJdJomw/T6j9cHGeaxtG6ZrZq8nPN4DrgO36d0WO41SDasXp1hNJu0i6QdLjkv4u6f+Vnfs48CjwW2As8MM8czeETTfJn24yswXJ648BZ/bzshzHqQYNplCXhaSdgNsI6b8ldkh00zDgLGAu8F3gZylFd3qkIZQusBpwnSQIa7rCzG7u3yU5jtNnzKCzYLaDUB93MXAgcDuwHqFjxDeAEcCvgdPyKtsSDaF0k/7xW/T3OhzHqQEF2+kCHwF+bWZ/TN4/LOmrwC3AJWZ2fF8mbwil6zjOAKZ4SncV4LEKWen9DX2d3JWu4zi1w4Di9UhrApZUyErv5/d1cle6juPUEAMrnE0XYHlJo8rel16PqJADYGZvZZ3Yla7jOLXDKKIjDUIRnckp8rSW20YOXepK13Gc2lI8m+4ltZzcla7jOLWlYErXzGpag9eVruM4NaS4dRVqhStdx3FqhwHFK+1YU1zpOo5TWwq205XURc8dhN8jtAq6HviRmb2TZ35Xuo7j1JBCpgHfRc9Kd3lgfUJa8Kcl7WhmC7JO7krXcZzaYWAFi9M1s4nLGiOpGfgP4FeEWg3fyDp/w5R2dBxngNJl2Y4CYWadZvZr4Grg03mudaXrOE5tGcCdI4B/AmvnucDNC47j1A4zj16owJWu4zi1pbi72CzsALyQ54KGMS9I2lvSU5JmSDq1v9fjOE41MKyzM9NRJJIeaV8EDiKEjmWmIXa6iSfwF8BewEzg35Kmmtnj/bsyx3H6RAFLO0r66zKGDCeEjI0EngTOzjN/QyhdQhPKGUkHCSRdBRwAuNJ1nKJTsJAxYOIyzi8CnidUIfthnhhdaBylOw54uez9TELLjKWQNAmYBDBhwoT6rMxxnF5jgFVxpyvpBWAB0Al0mFlbxfmJhO4OzyeiP5hZria3ZlZTs2ujKF2lyKL/U2Y2BZgC0NbWVqzvLI4zGLGaFDHfzczm9HD+bjPbr9o3BZDUQvhmPg543Mwq2/osk0ZxpM0E1ix7Px54tZ/W4jhOFSmaI03SREk/l7R6hXwd4H7gbuAqQsPKC3PPbw0QzpF8ejwN7AG8Avwb+HxPnyKSZgMvAqOBnj71isJAeI6B8AwwMJ6jGs+wlpmN6csEkm5O1pKFYQR7aYkpybfb8vmeB94mfBP+dcr5icDvCRu5V4Gv5t2NSrqYsJteq0J+F7Az8A/gPuDjwIeBY8wsc+HzhlC6AJL2BX4GNAMXmtn3Ml43rdKuU0QGwnMMhGeAgfEcA+EZ0pC0hpm9KmlV4FbgRDO7q+z8ikCXmb2T6JRzzWz9nPd4DLjDzE4ok21EcOzfVarNIGk54EHgFTPbI+v8jWJewMxuMrMNzOxDWRWu4ziDCzN7Nfn5BnAdwb5afn5+qdSimd0EtErKutMuMZbwzbuciYTd9W/K7vUecAWweZ7JG0bpOo7j9ISk5SWNKL0GPgY8WjFmdUlKXm9H0HFv5rzVUELN3HK2TX7+rUL+MrBSnskbJXqhL0xZ9pBCMBCeYyA8AwyM5xgIz1DJasB1iU5tAa4ws5slHQdgZpOBA4HjJXUQFOchlt+G+hKwSYVsZ+ANM3u5Qj4cmJtn8oax6TqO4zQCkqYQlPdHzewRSZ8mOOcuNrNjKsZOBj5iZltlnt+VruM4zgckoWGPEqIp3gRWAdqBbcojIZLyBS8DvzezE7PO7zZdx3GcMszseeCjwE0EpftnYGJK6Nluyfkb8sxfaKVb1Mpkki6U9IakR8tkoyTdKumZ5OfK/bnGZSFpTUl3SHpC0mOSTkrkhXkOScMk/UvSQ8kzfCeRryPpvuQZfidpSH+vdVlIapb0oKQbk/eFe4ZGwsymmdknzWxjM9vPzO5NGXObmW1mZrflmbuwSresMtk+hADlQyV9uH9XlZmLgb0rZKcCtycxhbcn7xuZDuC/zWxjYHvgS8nvv0jPsRjY3cy2ALYE9pa0PfBD4JzkGd4Gju3HNWblJOCJsvdFfIZBQWGVLmWVycxsCSEt74B+XlMmkmDutyrEBwClrJZLgE/VdVE5MbNZZvZA8noB4R/8OAr0HBYotc9uTQ4DdgeuTeQN/QwAksYDnyCJIU1Cpgr1DIOJIivdtMpk4/ppLdVgNTObBUGhAav283oyI2ltYCtCamShniP5Wj4deIOQ4fQsMNfMOpIhRfi7+hlwClCqLLMKxXuGQUORlW6mymRObZG0AiGc5mQzm9/f68lL0tV1S0KRpe2AjdOG1XdV2ZG0HyF+9P5yccrQhn2GwUaRkyMGWmWy1yWNNbNZksYSdl4NjaRWgsK93Mz+kIgL9xwAZjZX0p0E+/RISS3JTrHR/652AvZP6gwMA1Yk7HyL9AyDiiLvdP8NrJ94aYcAhwBT+3lNfWEqcGTy+khyhqHUm8RueAHwhJn9tOxUYZ5D0hhJI5PXywF7EmzTdxCC46HBn8HMTjOz8Wa2NuHfwF/N7DAK9AyDjUInR/S2Mll/I+lKQgGN0cDrwBmE5nZXAxMIaYifM7NKZ1vDIGlnQl3RR/jAlng6wa5biOeQtDnBydRM2IBcbWZnSlqX4JgdRagi9QUzW9x/K81GUtbwq2a2X1GfYTBQaKXrOI5TNIpsXnAcxykcrnQdx3HqiCtdx3GcOuJK13Ecp4640nUcx6kjrnSdwiPpziSxwXEaHle6juM4dcSVruM4Th1xpesAIGkDSdclxdUXSXpJ0jWSWiRNlGSSPivpYklvS5ov6XJJq1TMc4KkeyS9JWmupHslfaJizFmSlkjatky2fFKQ/h5JLWXyj0q6XdICSe9K+oukTWv/G3Gc2uBK1ylxI6H83/HAxwnFxxez9N/IzwjVqg4F/gfYnw9qtpZYm1DX9XPAwcA04EZJ+5SN+XYivyKpUgahIP3qwOdLJQkTZX078A7wBeDzwAjgbknlxY4cpzB4GrCDpNHAbOAAM4uKBiU5/XcAfzGzvcvkhwGXAXua2e0p1zURlPZNwHtmdkDZubWB6YRCLDcDVwCHmdkVZWNmAC+a2R5lshWB54DLzOzkRHYngJlN7M3zO0498Z2uA6G53nPA2ZK+KGn9bsZdXfH+GkKxmx1KAknbSLpR0uuElj7twF7AhuUXmtkLwHHAEcBFwKUVCnd94EPA5YmJoyUxOywE7gF27e3DOk5/4krXwcLXnb0IX/l/ADwt6TlJx1cMfb3iuiWE/lvjIDSrJJgDRgEnAjsC2xJ2ssNSbv0ngsIfCpxTca7UceICguIuP/YjdEdwnMJR5CLmThUxs+eAI5I6uVsAJwC/lPQC8F4ybLXya5I6xisDrySivYGVgIPMbGbZuOHd3PYXhLKKzwJTJO1kZu3JuTeTn6cBad1Wl2R/OsdpHHyn6yxF0qxxOvBfiag8UuCgiuGfI/wN3ZO8LynXkuJE0gaE7gZLIenzwOHAJILDbUvgzLIhTwEvAJsk7bArj4d783yO09/4TtcpFfM+F/gdMIOw+zyKYJP9KyFiAGATSRcRimNvAHwP+FuZE+225JpLJf0EGAt8h1DM/P0PeEnrAL8CLjCzaxLZ/xBsyreY2R1mZpK+BNyQ7KivBuYQdts7Ai9VdKxwnELgO10H4DWCYvwvQrudK4E1gP0qGh6eRGh6+Dvg+4Qws1JLGMzsMeAwYK1knlMIoWd3lcYkzrArknueVDb3jwn24N+WYn/N7CaCw2x5QhjaX4AfEULL7sFxCoiHjDnLpCxkbC8zS7OvOo6TEd/pOo7j1BFXuo7jOHXEzQuO4zh1xHe6juM4dcSVruM4Th1xpes4jlNHXOk6juPUEVe6juM4deT/AxWf4rel8+NyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 画图\n",
    "image_sz = np.shape(new_mask)[0]    # 图像边长\n",
    "maps = np.zeros((image_sz,image_sz))-9999    # 初始化画图数组\n",
    "for i in range(image_sz):\n",
    "    for j in range(image_sz):\n",
    "        idbin = (bin1d==binid[galid,i,j])    # 搜索二维bin的值在一维bin数组中的 位置\n",
    "        if len(bin1d[idbin])==1:    # 挑选出 一维2800的bin和二维76*76的bin相同的数值\n",
    "            maps[i,j] = mass[idbin]\n",
    "\n",
    "# 加mask\n",
    "masked_array = np.ma.array(maps,mask=new_mask)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "p = ax.imshow(masked_array , origin='lower',vmin=5.5,vmax=8.5)\n",
    "ax.set_xlabel('spaxel',fontsize=16)\n",
    "ax.set_ylabel('spaxel',fontsize=16)\n",
    "cb = fig.colorbar(p)\n",
    "cb.set_label(prop[0:len(prop)-8],fontsize=18,labelpad=10)\n",
    "#plt.xlim((0,53))    # x轴坐标范围\n",
    "#plt.ylim((0,53))    # y轴坐标范围"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAERCAYAAADYEnSgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXmYHVW1t99fDxlJCCEBQkIIyKAyYzMLRiZRERwQUWT+zAdXELxXEbwKXpzA6/w5YJRRJgEFIiIiCIIa0QBhngIEEgiQABkgUw/r+2NXw+mzq7urus85fap7vXnq6T6rdu3a1eleZ581ysxwHMdxakPDQC/AcRxnKOFK13Ecp4a40nUcx6khrnQdx3FqiCtdx3GcGuJK13Ecp4a40nUcx6khrnQdx3FqiCtdx3EKg6TPS3pY0kOSrpQ0ouz8sZIWS5qbHP9noNbaHa50HccpBJImA58DWsxsW6AROCJl6G/MbMfk+FUf7jOjn0vtEVe6juMUiSZgpKQmYBTwQhXucb6kmyRtXIW5aarGpLVgwoQJNm3atIFehuMMWu65554lZjaxP3O8772j7ZVX27Pd74E1DwOrS0QzzWxm5wsze17Sd4HngFXALWZ2S8pUH5O0D/AE8HkzW5Bz2f8DfBl4UNIpZnZFzut7REUteNPS0mJz5swZ6GU4zqBF0j1m1tKfOVp2GGH/+tPUTGMbJz3Z4/0krQf8FvgEsBS4BrjWzC4rGbM+8LqZrZF0InC4me2bd92SdgAuAbYDrgNONLMleedJw80LjuNUDQM6Mv7LwP7AM2a22Mxagd8Be3a5n9krZrYmeflL4F19WrfZ/UAL8C3gEOAhSR+VtHHakWfuwpoXHMepfwyj1bKZFzLwHLC7pFEE88J+QJePu5Immdmi5OUhwKN9vZmZtQFflXQdcDNhZ506lBy61JWu4zhVJeMutlfM7G5J1wL3Am3AfcBMSecAc8xsFvA5SYck518Fju3PPRNzxRnABIICv7c/84ErXcdxqohhtFfQb2RmZwNnl4nPKjl/JnBmJe4l6VDgfGAi8B3gLDNb2995Xek6jlNVOiiWs17SusD/A44EngI+amazKzW/K13HcaqGAe0FU7rAw8Ak4GfA6Wa2qpKTu9J1HKeqFG2nC7QD+5vZ7dWY3JWu4zhVw4DW4uUCPFXNyT1O13GcqmEY7RmPOmI6wbxQFXyn6zhO9TBoryt9OvC40nUcp2qEjDSnFFe6juNUEdGOBnoRfaFq+3NXuo7jVI3gSCuk0r1U0iUZx5qZDc868YAoXUmNhJzp583sYEmbAVcB4wlpdkdVIvPDcZyBJcTpFlLpPgUsrsbEA7XTPZWQxzw2eX0e8AMzu0rS+cAJwM8HaG2O41SQjmLudM+pdB3dTmoeMiZpCvBB4FfJawH7AtcmQy4BPlzrdTmOU3k6d7pZjqHCQOx0fwicDoxJXq8PLE3KqAEsBCanXZj0LpoBMHVqtsLIjuMMHIZo93SALtT0pyHpYOBlM7unVJwyNNVzaGYzzazFzFomTuxXFxHHcWpEhynTMVSo9U53L+AQSR8ARhBsuj8ExklqSna7U6hOsznHcWqMIdZa40AvIy+fAf5ZrclrutM1szPNbIqZTSO0Tv6LmR0J3A4clgw7BrihlutyHKc6hOSIhkxHvWBmF5jZ09Wav17idL8EXCXpG4Rq8BcM8Hocx6kQRXOSScrbHcLMLHMvtgFTumZ2B3BH8v3TwK4DtRbHcaqDmWi3+tnFZmQDPCPNcZyi0lGwna6ZTanm/K50nW65al5LJDtiizkpI/vOkwviCnpbbrIoZWTfufCJvSLZ8Vv9vaL3cNIJjjRXM6X4T8NxnKrR6UgbLCQlDN5FyCV4xMwezzvH4PlpOI5Tl7SbMh31gqR9JH1f0kZl8k2BfwOzCRm0j0iamXd+V7qO41SNzoy0LEcdcRxwmJm9WCa/GNgRuJvQLfhx4ARJR+eZ3M0LjuNUlY7iRS/sCswqFUjaGngP8DfgPWZmkkYRQlyPAS7NOrkr3UHOOQ99KJJNHbYkkh271exItlHTskz3eGzBxpFsfMrf2QaT40TDwv05OrkIBW8K97+8EfBkmWw64XF+aRY6bZrZSklXAp/NM7krXcdxqoYhWouXBjwCWFkm2yX5+tcy+bPAunkmd6XrOE7VMKOIyRELgHeWyd5NKNb1XJl8JJDtI2FC4X4ajuMUCdGR8agj7gKOlvQOAEmHAFsBf0oZux3wfJ7JXek6jlM1jLDTzXJkQdLnJT0s6SFJV0oaUXZ+uKTfSJon6W5J0/qw7HMJJoYHJb0A/A5YC3yv7F4NwCFArkwbNy8UlM/f94lItmh1bFp655hIRKPitPIrnoxLX2w5LG5Td828uK7H/qPiP5iFbZGI5hfi7MrGlA3O75/eLh6XkgrfrPgmYxpWR7JNmkdEsrQstQffiNeXZo/8yc5V6eIyaKmUI03SZOBzwDvNbJWkqwnVCi8uGXYC8JqZbSHpCEIrsPiPpQfM7ClJ+wJfA7YAbgW+bmYPlg3dD1hOzqqIrnQdx6kaRsULlDcBIyW1AqOIa28fSlCWEBIYfiJJnREHWTGzu4H39zLmz8A78swLrnQdx6kioQV7ZjUzQVJpcY+ZZvZmxpeZPS/pu8BzwCrgFjO7pWyOyQRHGGbWJmkZoSVYHCc5QLjSdRyniuRqOrnEzOIqS50zSesRdrKbAUuBayR92swu63LDmFy7XEnb5xkPYGYPZB1bU6WbGL3vBIYn977WzM6WdDEh26Mz9OJYM5tby7U5jlN5jIpmpO0PPGNmiwEk/Q7YEyhVuguBTYCFkpoIMbSv5rzPXPIpaiOHLq31TncNsK+ZvS6pGfibpD8m575oZtf2cK3TCxOGvx7LmlZkunbj5tci2fzWCZFsXOMbkazVOiLZZk3NkezVjtgxt6IjW+D8DsNeiWT3r10/07WL28ZGsjSH2/CG2DG3pi1+js/d98kur1e1x2N+2ZI5K3TQU8HOEc8Buyfpt6sIjqzyWqOzCGm5swktwP6S154LfItsSvdA3kqayExNlW7y8J2aoTk5qlah3XGcgcVMFdvpmtndkq4F7gXaCHUPZko6B5hjZrMIrb5+LWkeYYd7RB/u85WezkvamRAVsQvwGvDtPPPX3Kab1KO8hxCK8dPkB3kS8E1JZwG3AWeY2ZqUa2cAMwCmTp1aw1U7jtMXgiOtcmnAZnY2cHaZ+KyS86uBj1fshiVI2hz4ZjL/WuB/gW+b2dI889Q8OcLM2s1sR0Kr9V0lbQucCbyd8M4xntCoMu3amWbWYmYtEydOrNmaHcfpK6pocsRAIGmCpB8DjxAU7qXAVmb2pbwKFwYwIy1Z7B3AQWa2yAJrgIvwJpWOMygIjjRlOuoNSaMkfRV4CjgZuAXY3syON7OFfZ231tELE4FWM1sqaSTBG3mepElmtkiSgA8DD9VyXfXEcf86LpJdtOtFkWxp66hINropssiwoiPOxhrR0BrJhtEeyfYcEZdibFT8x9Ga4qd4nZR7pFy7VfOwSDZKcUhlvDpYbbED68eLDohk+49/JGW+bPuNpa0jI1m5w9KdZj1TtNKOiQl0BsFssSHwT+B0M/tbJeavtU13EnBJ8lANwNVmdqOkvyQKWYRwjRNrvC7HcapAFTLSasEjBJ/T48BJZnZ9JSevdfTCA8BOKfJ9a7kOx3FqRwEbU25JsIxsDFwo6cJexpuZZYtfxDPSHMepImbQ2lE4pfsPqhjK6krXcZyqEcwLxVK6Zvbuas7vSrfOWLxmnUh28F2npIyMx201+qVItvXwRZFsRUfsHNq4aVUkm9QU3+OBtXEm1+Ypv0X3ro0deGMVO/poiucbkWICTHPgbdm8OJJ9a8qNkezJ1rjk5WNr4r5u24yKa1FPaI6z/B58fXK8QKdbKpiRNigo1luQ4ziFooghY5J2lTQ+49jN8rZgd6XrOE4VCeaFLEcdMRs4qPOFpPUkvSpp75SxexJyCzLj5gXHcapKnfU/y0L5ghuAcYRaMf3Gla7jOFUjRC8UrgV7VXGlW2esbY9/QbcbF2eGpbFobewweseI2Dm0fmPsHFqd0aa2/bDYQdZqcb7YPvEwlnXEJSBbUwJzXkkJMRrXEF/7RkpHgvaO+Dk2bozLW7YOezmSzV8bl7Jc1h47Hcs55l/HR7JLdu0ttHNoUNDkiKriStdxnKpSQPNCVXGl6zhO1eiMXhgkVCRhwpWu4zhVpc4iE7JymqTDku+HERTuOZLKg8On5J3Yle4Ac+Adp/U6Jq2iWBq7jH0m07j3joxtuitikylPtMateUYofrNvT3n/H90Q7242aBwdyR5eGydlbNwYT9ia8hF1WlPc/mfdhrhq2bNtsc15hOIqaOMaV0ayRa3jItnUkV1bbt2/NNvf3SF3nRzJZu39k0zXFhUz0VZMpduSHKXs1c3YXDtgV7qO41SVApoXKhIa1h2udB3HqRpFtOmapYTjVBBXuo7jVJWiKd00JA0HxgLL0/o35qGQxhbHcYpBZ5xukWovdCJpU0k/k/QssBJ4EVgpab6kn0jqU3fcWrfrGQHcCQxP7n2tmZ0taTPgKkJTynuBo8ws9pIMQp5YuGGX11tNiSuFjU2pAJZG2i9uWiLEs23xj3bDxvj9d6vmODFgSXvsXGtN8SOkOc0WtcVrGZ/ytp/mNEubr91i71+j4gmbiTcmExvj53hszaSUe/S+L3nmtbg2ynazzkoZmamGyqCjiHG6kt4D3EDY3bYCjwHLk9dbAP8BfErSh8zs73nmrvVOdw2wr5ntAOwIHCRpd0IP+R+Y2ZaEPvIn1HhdjuNUATNo62jIdNQLksYCVxJqMJwErGtm25jZHma2DbAuQek2AldJGpNn/po+adLxt3O705wcBuwLXJvILyE0p3QcZxBQQPPCUcBGwKFm9gsz61L02cxWm9n5BD01Gfh0nslr/vYiqVHSXOBl4M+E9sZLzawtGbKQ8CBp186QNEfSnMWL4wLWjuPUFwW16b4fuMPM7uhpkJndDtwBfDDP5DVXumbWbmY7EjI5dgXekTasm2tnmlmLmbVMnDixmst0HKdCmCnTUUdsB9yeceztyfjMDFjImJktlXQHsDswTlJTstudAmQrq1XHTLvs2326rtyxBrDzegsyXftqW9xe5/m29SLZximZV/Na4+pm4xpix9fUpti51qxspfvS2v/0hzSnWRobNw2PZIvbY+fajiOei2RpWWrl/J5tM63jwUPOyTRusFFAR9p4suug54HMnYChxjtdSRMljUu+HwnsDzxKeLfozHM+huA1dByn4JhVzqYraWtJc0uO5ZJOKxszXdKykjFpYSS9MRrIFjIUggN6r/9ZQq13upOASyQ1EhT+1WZ2o6RHCF7AbwD3ARfUeF2O41QF0V6hyAQze5wQ9USiQ54HrksZepeZHdzf2/Xz+m6pqdI1sweAnVLkTxPsu47jDDKqZK/dD3jKzJ6txuTAtyR9McO42H7XC54G7DhO1chZe2GCpDklr2ea2cxuxh5BiKVNYw9J9xPssl8ws4ezLiDhBYJuzOqtz+WDcqVbJeZ/+sxIluZc23Jy3DamnNfbY0dQGhsMWx7JGlM+JS3riIsoNSvO7npbc+z4ejklIy0tWywraa1+sjrmsjJcaUWjYkfaCLVFsgWtvWeRrTMiWyr+Hrd8KZLNPvC8TNcWFgt23YwsMbPycooRkoYBhwDxH1nIaN3UzF6X9AHgemDLzCsAzCx3jdw81E8aiOM4g5IOlOnIwfuBe80sypk3s+WdCVhmdhPQLCluflclJI2S9OWe6jK40nUcp2pY4kjLcuTgk3RjWpC0kSQl3+9K0HGv9PtBsjMG+DqhPkMqbl5wHKeq5DAv9IqkUcABwP8tkZ0Y7mPnE0JPT5LURgj7OsKskivItsyeTrrSdRynqlQyesHMVlKWjJAo287vfwLUdQ8kV7pVYrMffy+Sibh/V3qVia5sPjJbnYk0p1laacf2lDfiCQ1xz7A0+uM0S6PSTrOsTEnJjnuhPf4ZNCt29L3W1vVnsO7w1dGYNB59vKr+mbrErGohY4XFla7jOFWlzorZDDiudB3HqSo1t6jWOa50HcepGoboqKMC5fWAK13HcaqKb3S74kq3AqQ5zbKyfG3v2WZPr8qWjbjt6Ocj2aiGOFtqfEPcI21g3Fn1xcqO+P/i9fYRfZorq9Ns05n/26f5AZ6dkaU0wAAz9Bxp7YRCPN16V13pOo5TXQbZVlfSBoS4o3lmtqL0nJktATbp6Xo3tjiOU1WK1jlC0naSPleePixpvKTfA4uAOcBiSV/OO3+PO11Jf8kxl5nZfnkX4DjO4MWAjo76UagZ+Q/gY8BPy+S/IvRDmw88AOwJfF3Sg2b2+6yT92ZeaCD7h4PC/WQdx6kyBtTRLjYjewCzzN4qgZcUsDkUeBDYzcxWS9oQuIeQklwZpWtm0/uy4u6QtAlwKaG9cQehXuaPJH0N+AzQmXr15aRC0IDSHwdZJWlQtve9Ze1x15AX29aNZFuPjOt/NKe40p5qjbPZ0so9DhbeOWxFJHvD5keyPy3t2odwizFxxuCjVDb7rBBOs24oYJzuZODXZbL9CBvLn3a2ZDezlyT9mtBiLDO1dqS1Af9lZvdKGgPcI+nPybkfmNl3a7wex3GqTfGU7jrA0jLZLoQnKe8SPI9qNqaUNFnS9yXNkfSMpG0T+WmSduvtejNbZGb3Jt+vIDSlzFB9wHGcYpLNiVZPjjRCyNfbymR7AkvN7MkyeRMQV/bvgcxKV9I2BHvGUYT2FFPhzQoumwKn5rmxpGmEfml3J6KTJT0g6UJJqX2HJM1IFP6cxYuzFYFxHGeAsYxH/TAHODqx2ZJsKLcH0gIL3kHOdj15drrfI+xMNwM+SlfH2T+A3bNOJGkd4LfAaWa2HPg54Z1lR0I4Rqox1cxmmlmLmbVMnJi1fZHjOAOGgXUo01FHnEfwOz0m6R/AbQQf1I9Sxh7MWxvHTOSx6b4b+GTSe6jc6/JSsshekdRMULiXm9nvIBikS87/Ergxx7qqxjOf+69I1h/n2ohNYkfN8pW9Zzw1EvcvS2PdxlWRLK20431r4vKM26Y4kQaz0yyNNzri7dbqjrgc5/CGuJdaxLBs/2dDg7pSqL1iZvdJOgz4NtACPA2cbWZ/Kx0n6SCCefTmPPPnUbo9/RZNIFRp75GkjcYFwKNm9v0S+SQzW5S8/AjwUI51OY5Tz9SX6SATZnY9oallT2NuBuKQoV7Io3T/BRxHejza4cDfM8yxF8Em/KCkuYnsy8AnJe1I+O+ZT0krDsdxCk4BlW41yaN0vw7cKukW4ArCj3J/SacSdqf79DZBsj1P+6wx4DG5juNUgQImR0hqAkYBr5tZR9m5jwEnEMwKDwPfMLNH8syf2ZFmZn8FPkxwpF1IUJ7nAnsDHzazXMZkx3GGBqFlT+9HHfEVYAnBbPomyQbzauAgYDvgCOCuJBIrM7mSI8zsD8AfJG0BbAC8YmaP55mj6DRNWlnV+TdaN3ZovbRmbKZrd1nnmUj2Rkq5wo2a4nusrLPf+krSanGfs7TebMs6miNZmiNynca4XGY5Y8bHoZsrXu17f7lNLz4vkj177Jf6PF9Nqa/IhCzsDdxmZi93CiSNBP6HEJN7OPBPwif884EvAp/NOnlmpSupoXOrbWbzCJkYpefXM7PXss7nOM7QIGMWez2xFcHhX8q+wFjg3MSBBnCRpP2B/fNMnidO966k6EOEpOnA/Xlu7DjOECBrYkR9KeYJwHNlst0Iq/xDmfyfkK/QRh6lOwG4X9InOgWSGiV9E7iVkDjhOI5TgoIjLctRPywDxpXJdiN0hbinTB7bn3ohj9LdCbgBuDJJ1d2ekIn2BeAMM3tf3ps7jjMEKN5Odx5wSOcLSesSai/c21lhrIQpwMvkILNN18xWAsdKuhm4iFDO7Clgj84iNoONLa/5ekXna2uLnTdT1i8vZhSzaGXsSNt5/IJIluY0e6U9ziob1xg7A5vJkGVVUNZYayZZe8Y9SLNix1wWKu1cKwzFS877JcFeexWh3sLhhBCyy1LG7kPOT/m5ohckjQLeBwwHXiGUNJsGDEql6zhOPylgnC5BuR5CULaHJ7I/EiIV3kTSZsB7gdPzTJ6nyti7gPsIsbpHEOJ1/wBcI+kXSUiF4zhOF2TZjl7nkbaWNLfkWC7ptLIxkvRjSfOSqoU7512vmbWb2ccIRbyOAvY0sw+aWfnHwWbgSNJ3wN2SZ6f7D0LJswPN7NlEdnRibvgZIbbtnXlu7jjOEKBC9tokJ2BHCE58Qt3b68qGvR/YMjl2I1Qw7LXWdzf3+xeh/EF3558AnpCUqdhXJ3kcad8B9ilRuJ03voLgZOvdOOk4jlMZ9gOeKtdHhD5ml1rgn8A4SZMqffMkcutQSbOA8jX0SB5H2ld7OPeMpL3z3LgIPPnx+JEr7VzLwm7rz880bkVHXCZyTEO5sxVGa3A4zVZ2rI1koxriUoyrUzLSWlMy8NrJ5oh87I0Nu7xu64gdpNtOfDGSzX68vBnB0CBHcsQESXNKXs80s5ndjD0CuDJFPhko9TIvTGSLUsbmRtLbgeMJZocNCG7Cv+aZo2I90ko7ZzqO4wBJD/bMjrQlZtbS2yBJwwiOrjPTTnezij4jaTRByR/PW80aZgP/DdxgZnGn1x7IG73wPuBEYGsg2laZ2eZ55nMcZwhQ+Rjc9xNiZl9KObcQ2KTk9RRyttPpRNJehIpihxGaVT4IfB/4T+D7nU0Y8pIneuEDhBKMo4C3A48RUuU2oQ9bbMdxhgaVil4o4ZOkmxYAZhEc/JK0O7CspEFCtvVKp0t6DLgL+BAhL6HFzHYghI31KwYujyPtq8BPgQ8kr79iZtOBbYBGQhxbj0jaRNLtkh6V9HBSKg1J4yX9WdKTydfUxpSO4xSQCmakJbkCBwC/K5GdKOnE5OVNhPY68whJDv/RhxWfS9CNHwU2NrNTK5kAlse88HbgLMKu1jqvNbMnJH2NoJSv7mWONuC/zOxeSWOAeyT9GTiWUErtXElnAGcABalbl87Y0bEDK43lq7taaTZcJy67+MCyuEv9LuvFDtMnVsaRKweOizsfzW+L39Maml6NZP1x+S5pj7OvJjRWNvsqzWmWxoiUMo4rOmJn4oiUTLP2lE3NlBFdA3WeXjkhGpPGtE2ydbCev2CQNV2toHkhyYxdv0x2fsn3Ro4yi93wMrAFoUfa1pIuM7Pn+znnm+TZ6XYAbclDLSa0YO/kBeI+8RFmtqjzHcPMVhDS5yYTwjwuSYZdQkjAcByn4GQ1LdRZ+cfJhF3uPOAbwHxJt0j6FMG82i/y7HQfJ6T8QkiSOE3S30l2r4TeZplJqq3vRGhfvGGn3cXMFknaoJtrZgAzAKZOTa0y6ThOvVGwIuZJJNb1wPVJ4sNxhE/jlwFrCHv3zSQp2YTmIs9O93LgHcn3ZxNsuQuBFwkFfs/KOpGkdQht2E8zs+VZrzOzmWbWYmYtEycOso9gjjNIKeBO903M7EUz+7aZbQ1MJ5hQVxGSxV6UdL6kA/PMmSc54qcl398jaTtC6MZI4NaszdkkNRMU7uUlIRcvdbZhT7JHcpVKcxynjqlThZoXM7sTuFPSyYQIihMIn7w/QwgmyESfkyPMbCHBO5gZSSK0wXjUzL5fcmoWoVTkucnXG/q6rmqTlqX2rpv+u2Lzb7NutuiWF9asG8k2H7kkkr3SFmdUNadkpE1rijstPbo2LgG5RXOctZXWb6w/TrNnWuO60Js1x8+RxusdsQNznYY4U2+dlM94aeUe57fGslGNXTPhxjT13jMN4CXGZBqX5nD7637fzXRt3VHHu9i+kvijZgIzJW1LSJrITG6lK+m9wB4EY/PzwD/M7I6Ml+9FSJ97UNLcRPZlgrK9WtIJhNjfj+ddl+M4dcogU7qlmNlDhGSJzORpTDkeuIZg1zDgNWC9cEp3AB83szjuqOsC/0b3gcX7ZV2L4zjFQQUrYp5EKeTBzKy7ZI2IPDvdHwO7EHaq15hZa2KfPZxQ2vFHyTnHcZwicxlv7c97Cr2w5LzRfYZcRB6l+yHgzKSUY7ijWStwebIL/kaOuRzHGSoU07ywmpD19kdCQ8qKkUfptgNPdnPucSq8sMHI6tb4xz2iuatT69Hl2eohb75O7DRrSPkcN7ohdvKMaIidQwvaypufwgjF4zZvjudrT2nH0qg80Yhdyeo0SyPNaZaV1pRCeY0pP9NlbV2bpLRZtmd9fU3shBz0FNOR9iVCXO6RhHDYS4GLkqLl/SbPX8YNwCe6OXcEIZjYcRynKwXrBmxm/2tm2xAc/38ATgIelfQ3SccneQZ9Jo/S/T2wv6Q/SDpW0vuTrzcRnGCzJO3befRnUY7jDCIKpnQ7MbPZZjaDUIbkOKCVECa7SNKFkvrUniyPeeHa5OsmhKSIcn6bfO00LGcOFnYcZ3Aiihe9UI6ZrSKYGC5NOgCfT8gnmA+ck3e+PEr3vXkndxxniFNMm25EUoPhKMKO9+2E9j9ze7yoG/Io3buAhtI2xEkniW2Bv5jZfX1ZwGDgng98M5KlZamNHBY7pspZf3hcEjGNCc1x1lZHikPnxbY4c60x5bPclsPjnl6bNsVlMZZlLF6yNqUOyJSmfpnCqk5DirXt7c2xw3LBiK5Z6svbsjnvNhkbZ/0NCQqqdCU1EaK2jgfeR3iSG4EvAn80sz7t4fMo3SsJFXaOThZ0IiE+F6BV0gfN7Na+LMJxnEFMwZSupO0JO9ojgQmENj2nA5eZWfwunJM8Snd3uhYW/yLwK0JZx5mEJm2udB3H6UIBzQtzCZXEriO06vl35wlJY9MuyFMtMY/S3YBQawFJWwCbAT8xsxWSLgKu6Olix3GGKMVTuhCqJ34yOXrjzU46WcijdJfzVpuM6YR2yQ8kr9tJ6Q7sOM4QxwoZvXB5NSfPo3T/AZwhqQ04jdAArpMtCAXNnYQ059oBt3++YvP/67VpkWzSyGWRbKd1notkmwx7JZKllXt8tSPuQbYyJWtrSspv0SgVq1sAQGPKmt+w+OE2au76c36jLc40W7I6W2nLW6b/MOPqCkyGYMxhAAAW8ElEQVTBdrpmVtUaMnmSI04HxhNq344AvlZy7hPA7Moty3GcwUKRO0fkRdK4pJ/aDt2NydM54klgK0nrm1n5VulUQtsex3GcrgwShZqR4cD+lHUsLiV3EfMUhYuZPZh3HsdxhgB1muI7kPS9FFQfSPKVX5b0UInsa5KelzQ3OT5QyzU5jlM9xNAyL2Shzz3S+sjFwE8Iecyl/MDMCtoEKp09b/lSijR2TJWzqr050/xbj3kpkq3XFPc0G9cYy0ZrbSQblVICsj2lfvPExtgV3ZryB7MyJSOtOaVUZFNKiY7+lIXsDys74vWl/W+Mb+yaDbjh8DhEsyOl3OVQZSgp1CzU9Lc76abZY0sfx3EGGQWtMlYtBmZLEXOypAcS88N63Q2SNEPSHElzFi+OO6Y6jlOHuNLtQj0o3Z8DbwN2JFTu+V53A81sppm1mFnLxIkTa7U+x3H6SkZ7blYTRBKSda2kxyQ9KmmPsvPTJS0r8RGdVY3H6g+1tulGmNmbxklJvyRU8XEcZ7BQ2V3sj4CbzewwScOAUSlj7jKzgyt61woy4EpX0iQzW5S8/AjwUE/ji8I/Djwv07hP331Cr2PGD4udYWmsScmeeqMjzpZ6w2KHXkNKlbr1G1ZFsvlt8bVbNsUOqNUpf2nLOmIH3gaNcebWso74vus2jIxklWZ1ys9gcUdcjnJBa9cQzHeOeiEas7wtXu8Fu1zc98UVmEqlASfFZvYh9C/DzNYC8S/VwLKKkEa8qLsBtQ4Zu5KQuba1pIWSTgC+I+lBSQ8QCqVXLlfWcZwBJ4d5YUKnzyY5ZpRNtTmwGLhI0n2SfiUpLd96D0n3S/qjpG1yr1faNelwnmXsZpKO7nxtZsvN7Cgze7S7a2q60zWztIo9F9RyDY7j1JB8TrIlZtbSw/kmYGfgFDO7W9KPgDOAr5aMuRfY1MxeT2L+rwe2zLnq2YQuEVcAJM79p4BDzeyusrF7Eso/lofBdks9ONIcxxnMVC56YSGw0MzuTl5fS1DCb90q7DRfT76/CWiWNCHnisuDrBuAcaSHbefGla7jOFWjkhlpZvYisEDS1oloP+CRLveTNpJCuThJuxJ0XFxWbwAZcEfaUGfqyN77ZqX14BreGJdinDwsnmtycywbk+IgW23xm/hqi7PFWlNki1McZC+2xw6ozZvivm6vtcdOwuGq7K9lmmOuNcVp1pxS2vGF1jhsvNw5uaR1TKZ1nPDvYyPZUHCuqaOi4QunAJcnkQtPA8clrcMws/OBw4CTkhK0q4AjzFLSIwcQV7qO41SPCic+mNlcoNzue37J+Z8QSg3ULa50HcepKoOo9kJFnsSVruM41aWYSvc0SYcl3w8jPMU5ksrrD0zJO7ErXcdxqkpBd7otxGaMvboZm+sJXekOMN/a/nddXn/1wQ9HY1a1x1lgGw2Lywkua4+zoNqbY+dQmtNsXIpzLa2047iUEpCtFgfBTGyIHWRpJSDbM/5Frkxx1o1q6L1UJsCaFKdZGv9cHdfzWJviOHxuTbdNAd5kKDjIMlM8pVuR0LDucKXrOE71KGA3YLOU7qsVxJWu4zhVozNOt2hI2h9Yk5KBVjpmb2CYmd2WZ25PjnAcp7qYZTvqBEnTgT8B03oZuilwi6TubL2puNJ1HKeqFLBH2vHAw2b2654GmdllwIPAZ/JM7uaFASbNcVbOFqNejmT3Ld8kku07/rFI1pjixWhPcXytTikL2Z7ynpw6X4rDbePG2OGWZtpb3B5fu35jSqnIlJ3QG21xVt7ohvguC9rikqt7jci232hM0QblGYI/2fmKTHMNSYrZFeLdwGUZx84CjswzuStdx3GqStEcacDGwDMZxz6TjM+MK13HcapKAZXuWiBr1fyRQPzRrAfcpus4TvUwCudIA+YDu2Ycu2syPjO17hxxoaSXJT1UIhsv6c+Snky+dtsN2HGc4lFAR9otwOGStuhpUHL+cODmPJPX2rxwMaECUGmV9TOA28zsXElnJK+/VON1DRhf3+76Xsec89CHItkWo+MW9IvWjotkwxQ7m8Y0rI5kHSnvv80p177SEXdHGdcYZ58tac+W1LM25b4r4tumlpRsTXH+NXfEF6/oiD8pnr80rmu9rD12uN392rRIdv27fxYv0Ome+lKoWfgBMAP4s6Tjzez28gFJWNmFQBuhWWZmat2u505J08rEhwLTk+8vAe5gCCldxxnMFDE5wsyel/Rp4CrgVknPAHOB5cAYYCdgM2ANcLiZPZ9n/npwpG3Y2Q3YzBZJ2qC7gUmjuhkAU6dOrdHyHMfpM2aVLmJeE8xslqTdgG8ABxGaYnbSCtwIfNXMHsg7dz0o3cyY2UxgJkBLS0vx/icdZyhS0L9UM3sQOFTSSGBrYCxht/u4mcUVojJSD0r3JUmTkl3uJCDOBBjiNKT81r7ePjySTR3+aiR7qXXdSDZ71dsi2fbrLIxkGzUti2RpiRBzVm4eybYduSCSNabEDqXZZV9uGxtfm5JakVYtrVlxrZK0a1d0xC2Qnl8T28R3Hhc/R7mN/axtfx+Ncd6iaOaFchIFOzftnKRmQpfga7POVw8hY7OAY5LvjwFuGMC1OI5TSQzosGxHgZC0naQfAC8Av8lzbU13upKuJDjNJkhaCJwNnAtcLekE4Dng47Vck+M4VaZY+rRbJI0FPkWozfAugp/wfuo8euGT3Zzar5brcByndhTdvJCEh50AfJSQgWbAr4DzzOzpvPPVg03XcZxBTBGjFyRNBo4FjiOEh70KXADcRTAn/KkvChdc6RaCr2x7YyQ7+d5PRbJX2+LEhTTGNsXJEa9lvDaNMY3xfAta45Y2aQkOmzTHzr8Oi511S9pih+ArrfGa13RkS8pIazE0vCFOrHAnWT8pYJUxSTcBBxBW/kfgi8CNZtYqKfZC58SVruM4VSMkRxRM64a43HnAJ8zsvkpPXg/RC47jDGY6Mh71w/WErhD/kvRHSUdIimMM+4jvdB3HqSpF2+ma2UclTSCEsB4HXAEsl3Q1cGd/5/edruM41cNyHHWEmS0xs++Z2bbAXsA1wCcIxboMOFjSNn2Z23e6BaXSLWJOvOeoSNYyKi6e/4HNH4pkn7/vE5HsBzvF8eJfvD8OwU6rUPbcmtgJl0ZaVl5WftlyaST7zJyj+zyf0x3FrL1QipnNBmZLOpWgeI8nRDYcI2ke8Fsz+3LW+Xyn6zhOdalgEXNJ4yRdK+kxSY9K2qPsvCT9WNI8SQ9I2rlyj2ErzewiM9sbeDvwXUI9hlxVEV3pOo5TPSy068lyZORHwM1m9nZgB+DRsvPvB7ZMjhnAzyv0JF0wsyfM7EvAFOAjea51pes4TnWp0E43ScPdh5CkgJmtNbOlZcMOBS61wD+BcUkhrcxIOlxS3G479dGs3cxm5Znfla7jONUluyNtgqQ5JceMspk2BxYDF0m6T9KvJJVnyEwGSkvDLUxkebgS2LvzhaSxkv4h6V0550nFHWkOABOHrYhkaU6zNNKcZmmklV08fqu/Z7q2P2R1kKU515z+o47MtoMlZtbSw/kmYGfgFDO7W9KPCO29vlp6u5Tr8nryyudoBnYH4rTIPuA7XcdxqodRyeSIhcBCM7s7eX0tQQmXjyk1DUwhlF+sG1zpOo5TNYQhy3b0hpm9CCyQtHUi2g94pGzYLODoJIphd2BZZzuwesHNC47jVJfKZqSdAlwuaRjwNHCcpBPDbex84CbgA4TaCSsJGWV1hStdx3GqSwWVrpnNBcrtvueXnDfgsxW41dHJThlgBMFQcrKkD6cvy07NOrGsTvKiJc0HVgDtQFsvBnVaWlpszpw5tVia4wxJJN3T299hb6w7amPbfevPZBp7y9xz+n2/SiDliBoOmFlK3dJuqLed7nvNbMlAL8JxnMqRI3qhXtismpPXm9J1HGdQkT3Ft14ws2erOX89RS8YcIuke1KCogGQNKMzcHrx4sU1Xp7jOLkxKlp7od5IoiTialE9UE9Kdy8z25mQO/1ZSfuUDzCzmWbWYmYtEydOrP0KHcfJT/GKmPdKomw/Raj9cHGea+tG6ZrZC8nXl4HrgF0HdkWO41SCSsXp1hJJe0u6QdIjkv4m6f+WnHsf8BDwa2AScF6euevCppvkTzeY2Yrk+wOBcwZ4WY7jVII6U6i9IWkv4FZC+m8neyS6aQTwDWAp8HXghylFd3qkLpQusCFwnSQIa7rCzG4e2CU5jtNvzKC9YLaDUB93DXAYcBuwBaFjxFeAMcAvgDPzKttO6kLpJv3jdxjodTiOUwUKttMFdgN+YWa/T14/IOkLwC3AJWZ2Un8mrwul6zjOIKZ4Snd94OEyWefrG/o7uStdx3GqhwHF65HWAKwtk3W+Xt7fyV3pOo5TRQyscDZdgNGSxpe87vx+TJkcADN7NevErnQdx6keRhEdaRCK6JyfIv9diszIoUtd6TqOU12KZ9O9pJqTu9J1HKe6FEzpmllVa/C60nUcp4oUt65CtXCl6zhO9TCgeKUdq4orXcdxqkvBdrpJEfOeFr2K0CroeuA7ZvZ6nvld6TqOU0UKmQZ8Jz0r3dHAloS04I9I2tPMVmSd3JWu4zjVw8AKFqdrZtN7GyOpEfg/wM8JtRq+knX+uint6DjOIKXDsh0FwszazewXwNXAR/Jc60rXcZzqMog7RwD/AKblucDNC47jVA8zj14ow5Wu4zjVpbi72CzsAczPc0HdmBckHSTpcUnzJJ0x0OtxHKcSGNbenukoEkmPtM8AhxNCxzJTFzvdxBP4U+AAYCHwb0mzzOyRgV2Z4zj9ooClHSX9pZchowghY+OAx4Bz88xfF0qX0IRyXtJBAklXAYcCrnQdp+gULGQMmN7L+dXAM4QqZOflidGF+lG6k4EFJa8XElpmdEHSDGAGwNSpU2uzMsdx+owBVsGdrqT5wAqgHWgzs5ay89MJ3R2eSUS/M7NcTW7NrKpm13pRukqRRf9TZjYTmAnQ0tJSrM8sjjMUsaoUMX+vmS3p4fxdZnZwpW8KIKmJ8Ml8MvCImZW39emVenGkLQQ2KXk9BXhhgNbiOE4FKZojTdJ0ST+WtFGZfDPgHuAu4CpCw8oLc89vdRDOkbx7PAHsBzwP/Bv4VE/vIpIWA88CE4Ce3vWKwmB4jsHwDDA4nqMSz7CpmU3szwSSbk7WkoURBHtpJzOTT7el8z0DvEb4JPyLlPPTgd8SNnIvAF/IuxuVdDFhN71pmfxO4N3A34G7gfcB7wSON7PMhc/rQukCSPoA8EOgEbjQzL6Z8bo55XadIjIYnmMwPAMMjucYDM+QhqSNzewFSRsAfwZOMbM7S86PBTrM7PVEp/zIzLbMeY+HgdvN7OQS2dsJjv07O2szSBoJ3Ac8b2b7ZZ2/XswLmNlNZraVmb0tq8J1HGdoYWYvJF9fBq4j2FdLzy/vLLVoZjcBzZKy7rQ7mUT45F3KdMLu+lcl91oFXAFsn2fyulG6juM4PSFptKQxnd8DBwIPlY3ZSJKS73cl6LhXct5qOKFmbim7JF//WiZfAKybZ/J6iV7oDzN7H1IIBsNzDIZngMHxHIPhGcrZELgu0alNwBVmdrOkEwHM7HzgMOAkSW0ExXmE5behPgdsUyZ7N/CymS0ok48CluaZvG5suo7jOPWApJkE5f0eM3tQ0kcIzrmLzez4srHnA7uZ2U6Z53el6ziO8xZJaNhDhGiKV4D1gVbgXaWREEn5ggXAb83slKzzu03XcRynBDN7BngPcBNB6f4RmJ4Sevbe5PwNeeYvtNItamUySRdKelnSQyWy8ZL+LOnJ5Ot6A7nG3pC0iaTbJT0q6WFJpybywjyHpBGS/iXp/uQZ/ieRbybp7uQZfiNp2ECvtTckNUq6T9KNyevCPUM9YWZzzOxDZvYOMzvYzP6ZMuZWM9vOzG7NM3dhlW5JZbL3EwKUPynpnQO7qsxcDBxUJjsDuC2JKbwteV3PtAH/ZWbvAHYHPpv8/Iv0HGuAfc1sB2BH4CBJuwPnAT9InuE14IQBXGNWTgUeLXldxGcYEhRW6VJSmczM1hLS8g4d4DVlIgnmfrVMfCjQmdVyCfDhmi4qJ2a2yMzuTb5fQfiDn0yBnsMCne2zm5PDgH2BaxN5XT8DgKQpwAdJYkiTkKlCPcNQoshKN60y2eQBWksl2NDMFkFQaMAGA7yezEiaBuxESI0s1HMkH8vnAi8TMpyeApaaWVsypAi/Vz8ETgc6K8usT/GeYchQZKWbqTKZU10krUMIpznNzJYP9HryknR13ZFQZGlX4B1pw2q7quxIOpgQP3pPqThlaN0+w1CjyMkRg60y2UuSJpnZIkmTCDuvukZSM0HhXm5mv0vEhXsOADNbKukOgn16nKSmZKdY779XewGHJHUGRgBjCTvfIj3DkKLIO91/A1smXtphwBHArAFeU3+YBRyTfH8MOcNQak1iN7wAeNTMvl9yqjDPIWmipHHJ9yOB/Qm26dsJwfFQ589gZmea2RQzm0b4G/iLmR1JgZ5hqFHo5Ii+ViYbaCRdSSigMQF4CTib0NzuamAqIQ3x42ZW7myrGyS9m1BX9EHesiV+mWDXLcRzSNqe4GRqJGxArjazcyRtTnDMjidUkfq0ma0ZuJVmIylr+AUzO7iozzAUKLTSdRzHKRpFNi84juMUDle6juM4NcSVruM4Tg1xpes4jlNDXOk6juPUEFe6TuGRdEeS2OA4dY8rXcdxnBriStdxHKeGuNJ1AJC0laTrkuLqqyU9J+kaSU2SpksySR+TdLGk1yQtl3S5pPXL5jlZ0mxJr0paKumfkj5YNuYbktZK2qVENjopSD9bUlOJ/D2SbpO0QtIbkv4kadvq/0Qcpzq40nU6uZFQ/u8k4H2E4uNr6Po78kNCtapPAv8NHMJbNVs7mUao6/px4BPAHOBGSe8vGfO1RH5FUqUMQkH6jYBPdZYkTJT1bcDrwKeBTwFjgLsklRY7cpzC4GnADpImAIuBQ80sKhqU5PTfDvzJzA4qkR8JXAbsb2a3pVzXQFDaNwGrzOzQknPTgLmEQiw3A1cAR5rZFSVj5gHPmtl+JbKxwNPAZWZ2WiK7A8DMpvfl+R2nlvhO14HQXO9p4FxJn5G0ZTfjri57fQ2h2M0enQJJ75J0o6SXCC19WoEDgK1LLzSz+cCJwNHARcClZQp3S+BtwOWJiaMpMTusBGYD+/T1YR1nIHGl62Dh484BhI/83waekPS0pJPKhr5Udt1aQv+tyRCaVRLMAeOBU4A9gV0IO9kRKbf+A0HhDwd+UHaus+PEBQTFXXocTOiO4DiFo8hFzJ0KYmZPA0cndXJ3AE4GfiZpPrAqGbZh6TVJHeP1gOcT0UHAusDhZrawZNyobm77U0JZxaeAmZL2MrPW5NwrydczgbRuq2uzP53j1A++03W6kDRrnAv8ZyIqjRQ4vGz4xwm/Q7OT153KtVNxImkrQneDLkj6FHAUMIPgcNsROKdkyOPAfGCbpB12+fFAX57PcQYa3+k6ncW8fwT8BphH2H0eS7DJ/oUQMQCwjaSLCMWxtwK+Cfy1xIl2a3LNpZK+B0wC/odQzPzNN3hJmwE/By4ws2sS2X8TbMq3mNntZmaSPgvckOyorwaWEHbbewLPlXWscJxC4DtdB+BFgmL8T0K7nSuBjYGDyxoenkpoevgb4FuEMLPOljCY2cPAkcCmyTynE0LP7uwckzjDrkjueWrJ3N8l2IN/3Rn7a2Y3ERxmowlhaH8CvkMILZuN4xQQDxlzeqUkZOwAM0uzrzqOkxHf6TqO49QQV7qO4zg1xM0LjuM4NcR3uo7jODXEla7jOE4NcaXrOI5TQ1zpOo7j1BBXuo7jODXk/wMm/Hy0Dw+KQQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 不加mask,只去除无效数据\n",
    "masked_array = np.ma.array(maps,mask=(maps<-666))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "p = ax.imshow(masked_array , origin='lower',vmin=5.5,vmax=8.5)\n",
    "ax.set_xlabel('spaxel',fontsize=16)\n",
    "ax.set_ylabel('spaxel',fontsize=16)\n",
    "cb = fig.colorbar(p)\n",
    "cb.set_label(prop[0:len(prop)-8],fontsize=18,labelpad=10)\n",
    "#plt.xlim((0,53))    # x轴坐标范围\n",
    "#plt.ylim((0,53))    # y轴坐标范围"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}