GOES_investigate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # Import the Matplotlib package\n",
    "from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit\n",
    " \n",
    "import numpy as np # Import the Numpy package\n",
    " \n",
    "from remap import remap2 # Import the Remap function\n",
    " \n",
    "# from cpt_convert import loadCPT # Import the CPT convert function\n",
    "from matplotlib.colors import LinearSegmentedColormap # Linear interpolation for color maps\n",
    " \n",
    "import datetime # Library to convert julian day to dd-mm-yyyy\n",
    " \n",
    "from matplotlib.patches import Rectangle # Library to draw rectangles on the plot\n",
    " \n",
    "from netCDF4 import Dataset # Import the NetCDF Python interface\n",
    "\n",
    "import glob\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['C:\\\\Projects\\\\RD\\\\OPIR_GOES\\\\GOES\\\\from_jeff\\\\OR_ABI-L1b-RadC-M3C06_G16_s20172920002122_e20172920004501_c20172920004549.nc', 'C:\\\\Projects\\\\RD\\\\OPIR_GOES\\\\GOES\\\\from_jeff\\\\OR_ABI-L1b-RadC-M3C06_G16_s20172920007122_e20172920009501_c20172920009546.nc', 'C:\\\\Projects\\\\RD\\\\OPIR_GOES\\\\GOES\\\\from_jeff\\\\OR_ABI-L1b-RadC-M3C06_G16_s20172920012122_e20172920014501_c20172920014543.nc', 'C:\\\\Projects\\\\RD\\\\OPIR_GOES\\\\GOES\\\\from_jeff\\\\OR_ABI-L1b-RadC-M3C06_G16_s20172920017122_e20172920019501_c20172920019548.nc']\n"
     ]
    }
   ],
   "source": [
    "root_dir = r'C:\\Projects\\RD\\OPIR_GOES\\GOES\\ABI-L1b-RadC\\2017\\282\\05'\n",
    "root_dir = r'C:\\Projects\\RD\\OPIR_GOES\\GOES\\from_jeff'\n",
    "b6_paths = glob.glob(root_dir + r'\\OR_ABI-L1b-RadC-M3C06_G16_*')\n",
    "\n",
    "print(b6_paths)\n",
    "\n",
    "# store some variables\n",
    "time_ls = []\n",
    "data_ls = []\n",
    "grid_ls = []\n",
    "date_ls = []\n",
    "\n",
    "for path in b6_paths:\n",
    "     \n",
    "    # Getting information from the file name ==============================================================\n",
    "    # Search for the Scan start in the file name\n",
    "    Start = (path[path.find(\"_s\")+2:path.find(\"_e\")])\n",
    "    # Search for the GOES-16 channel in the file name\n",
    "    Band = int((path[path.find(\"M3C\" or \"M4C\")+3:path.find(\"_G16\")]))\n",
    "    # Create a GOES-16 Bands string array\n",
    "    Wavelenghts = ['[]','[0.47 μm]','[0.64 μm]','[0.865 μm]','[1.378 μm]','[1.61 μm]','[2.25 μm]','[3.90 μm]','[6.19 μm]','[6.95 μm]','[7.34 μm]','[8.50 μm]','[9.61 μm]','[10.35 μm]','[11.20 μm]','[12.30 μm]','[13.30 μm]']\n",
    "     \n",
    "    # Converting from julian day to dd-mm-yyyy\n",
    "    year = int(Start[0:4])\n",
    "    dayjulian = int(Start[4:7]) - 1 # Subtract 1 because the year starts at \"0\"\n",
    "    dayconventional = datetime.datetime(year,1,1) + datetime.timedelta(dayjulian) # Convert from julian to conventional\n",
    "    date = dayconventional.strftime('%d-%b-%Y') # Format the date according to the strftime directives\n",
    "    date_ls.append(date)\n",
    "    \n",
    "    time = Start [7:9] + \":\" + Start [9:11] + \":\" + Start [11:13] + \" UTC\" # Time of the Start of the Scan\n",
    "    time_ls.append(time)\n",
    "    \n",
    "    # Get the unit based on the channel. If channels 1 trough 6 is Albedo. If channels 7 to 16 is BT.\n",
    "    if Band <= 6:\n",
    "        Unit = \"Reflectance\"\n",
    "    else:\n",
    "        Unit = \"Brightness Temperature [°C]\"\n",
    "\n",
    "    # Open the file using the NetCDF4 library\n",
    "    nc = Dataset(path)\n",
    "     \n",
    "    # Get the latitude and longitude image bounds\n",
    "    geo_extent = nc.variables['geospatial_lat_lon_extent']\n",
    "    min_lon = float(geo_extent.geospatial_westbound_longitude)\n",
    "    max_lon = float(geo_extent.geospatial_eastbound_longitude)\n",
    "    min_lat = float(geo_extent.geospatial_southbound_latitude)\n",
    "    max_lat = float(geo_extent.geospatial_northbound_latitude)\n",
    "     \n",
    "    # Choose the visualization extent (min lon, min lat, max lon, max lat)\n",
    "    #extent = [-85.0, -5.0, -60.0, 12.0]\n",
    "    extent = [min_lon, min_lat, max_lon, max_lat]\n",
    "     \n",
    "    # Choose the image resolution (the higher the number the faster the processing is)\n",
    "    resolution = 2.0\n",
    "     \n",
    "    # Calculate the image extent required for the reprojection\n",
    "    H = nc.variables['goes_imager_projection'].perspective_point_height\n",
    "    x1 = nc.variables['x_image_bounds'][0] * H\n",
    "    x2 = nc.variables['x_image_bounds'][1] * H\n",
    "    y1 = nc.variables['y_image_bounds'][1] * H\n",
    "    y2 = nc.variables['y_image_bounds'][0] * H\n",
    "     \n",
    "    # Call the reprojection function\n",
    "    raw,grid = remap2(path, extent, resolution, x1, y1, x2, y2)\n",
    "    grid_ls.append(grid)\n",
    "     \n",
    "    # Read the data returned by the function\n",
    "    if Band <= 6:\n",
    "        data = grid.ReadAsArray()\n",
    "    else:\n",
    "        # If it is an IR channel subtract 273.15 to convert to ° Celsius\n",
    "        data = grid.ReadAsArray() - 273.15\n",
    "        \n",
    "    data_ls.append(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAACxCAYAAADEfcp5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADV5JREFUeJzt3X+s3XV9x/Hna6WU+YPRihDSNgOW/iFLNmQ3pQnLwmQrpVtWlmiCWUbDSLpsmmi2ZKkzGZv+45bMLSQO001iSVRkqKExbHipLv4FtChCEWuvyOSuDY0rIosJir73x3lfPZbbH/dHzz3tfT6Sk+/3+z6f7zmf7yc9fd3v5/s9OakqJEn6haXugCRpPBgIkiTAQJAkNQNBkgQYCJKkZiBIkoAlCIQkW5IcTDKVZOeo31+SNLuM8nsISVYA3wR+F5gG9gHvrKqvj6wTkqRZjfoMYSMwVVXPVtUPgXuBbSPugyRpFqMOhLXA80Pb012TJC2x80b8fpml9nNzVkl2ADsAVrDiN17HhaPolySdM17mxe9W1Zvnut+oA2EaWD+0vQ44PNygqnYBuwAuzJq6NjeMrneSdA54uO7/7/nsN+opo33AhiRXJDkfuAXYM+I+SJJmMdIzhKp6Ncm7gYeAFcDdVfX0KPsgSZrdqKeMqKoHgQdH/b6SpJPzm8qSJMBAkCQ1A0GSBBgIkqRmIEiSAANBktQMBEkSYCBIkpqBIEkCDARJUjMQJEmAgSBJagaCJAkwECRJzUCQJAEGgiSpGQiSJMBAkCQ1A0GSBBgIkqRmIEiSAANBktQMBEkSYCBIktqCAiHJc0meSvJEkv1dW5NkMsmhXq7uepLcmWQqyZNJrlmMA5AkLY7FOEP47aq6uqomensnsLeqNgB7exvgJmBDP3YAdy3Ce0uSFsmZmDLaBuzu9d3AzUP1e2rgEeCiJJedgfeXJM3DQgOhgC8keTzJjq5dWlVHAHp5SdfXAs8P7TvdNUnSGDhvgftfV1WHk1wCTCb5xknaZpZavabRIFh2AFzA6xbYPUnS6VrQGUJVHe7lUeBzwEbghZmpoF4e7ebTwPqh3dcBh2d5zV1VNVFVEytZtZDuSZLmYN6BkOT1Sd44sw5sBg4Ae4Dt3Ww78ECv7wFu7buNNgEvzUwtSZKW3kKmjC4FPpdk5nU+WVX/mWQfcF+S24HvAO/o9g8CW4Ep4AfAbQt4b0nSIpt3IFTVs8Cvz1L/X+CGWeoFvGu+7ydJOrP8prIkCTAQJEnNQJAkAQaCJKkZCJIkwECQJDUDQZIEGAiSpGYgSJIAA0GS1AwESRJgIEiSmoEgSQIMBElSMxAkSYCBIElqBoIkCTAQJEnNQJAkAQaCJKkZCJIkwECQJDUDQZIEGAiSpHbKQEhyd5KjSQ4M1dYkmUxyqJeru54kdyaZSvJkkmuG9tne7Q8l2X5mDkeSNF+nc4bwcWDLcbWdwN6q2gDs7W2Am4AN/dgB3AWDAAHuAK4FNgJ3zISIJGk8nDIQqurLwLHjytuA3b2+G7h5qH5PDTwCXJTkMuBGYLKqjlXVi8Akrw0ZSdISmu81hEur6ghALy/p+lrg+aF20107UV2SNCbOW+TXyyy1Okn9tS+Q7GAw3cQFvG7xeiZJOqn5niG80FNB9PJo16eB9UPt1gGHT1J/jaraVVUTVTWxklXz7J4kaa7mGwh7gJk7hbYDDwzVb+27jTYBL/WU0kPA5iSr+2Ly5q5JksbEKaeMknwKuB64OMk0g7uFPgTcl+R24DvAO7r5g8BWYAr4AXAbQFUdS/JBYF+3+0BVHX+hWpK0hFI161T+WLgwa+ra3LDU3ZCks8rDdf/jVTUx1/38prIkCTAQJEnNQJAkAQaCJKkZCJIkwECQJDUDQZIEGAiSpGYgSJIAA0GS1AwESRJgIEiSmoEgSQIMBElSMxAkSYCBIElqBoIkCTAQJEnNQJAkAQaCJKkZCJIkwECQJDUDQZIEGAiSpHbKQEhyd5KjSQ4M1f42yf8keaIfW4eee1+SqSQHk9w4VN/StakkOxf/UCRJC3E6ZwgfB7bMUv+nqrq6Hw8CJLkKuAX41d7nX5KsSLIC+AhwE3AV8M5uK0kaE+edqkFVfTnJ5af5etuAe6vqFeDbSaaAjf3cVFU9C5Dk3m779Tn3WJJ0RizkGsK7kzzZU0qru7YWeH6ozXTXTlSXJI2J+QbCXcCvAFcDR4B/7HpmaVsnqb9Gkh1J9ifZ/yNemWf3JElzNa9AqKoXqurHVfUT4F/52bTQNLB+qOk64PBJ6rO99q6qmqiqiZWsmk/3JEnzMK9ASHLZ0OYfAjN3IO0BbkmyKskVwAbgMWAfsCHJFUnOZ3Dhec/8uy1JWmynvKic5FPA9cDFSaaBO4Drk1zNYNrnOeBPAarq6ST3MbhY/Crwrqr6cb/Ou4GHgBXA3VX19KIfjSRp3lI161T+WLgwa+ra3LDU3ZCks8rDdf/jVTUx1/38prIkCTAQJEnNQJAkAQaCJKkZCJIkwECQJDUDQZIEGAiSpGYgSJIAA0GS1AwESRJgIEiSmoEgSQIMBElSMxAkSYCBIElqBoIkCTAQJEnNQJAkAQaCJKkZCJIkwECQJDUDQZIEGAiSpHbKQEiyPsmXkjyT5Okk7+n6miSTSQ71cnXXk+TOJFNJnkxyzdBrbe/2h5JsP3OHJUmaq9M5Q3gV+MuqeguwCXhXkquAncDeqtoA7O1tgJuADf3YAdwFgwAB7gCuBTYCd8yEiCRp6Z0yEKrqSFV9pddfBp4B1gLbgN3dbDdwc69vA+6pgUeAi5JcBtwITFbVsap6EZgEtizq0UiS5m1O1xCSXA68FXgUuLSqjsAgNIBLutla4Pmh3aa7dqK6JGkMnHYgJHkD8BngvVX1/ZM1naVWJ6kf/z47kuxPsv9HvHK63ZMkLdBpBUKSlQzC4BNV9dkuv9BTQfTyaNengfVDu68DDp+k/nOqaldVTVTVxEpWzeVYJEkLcDp3GQX4GPBMVX146Kk9wMydQtuBB4bqt/bdRpuAl3pK6SFgc5LVfTF5c9ckSWPgvNNocx3wx8BTSZ7o2l8DHwLuS3I78B3gHf3cg8BWYAr4AXAbQFUdS/JBYF+3+0BVHVuUo5AkLViqXjONPzYuzJq6NjcsdTck6azycN3/eFVNzHU/v6ksSQIMBElSMxAkSYCBIElqBoIkCRjzu4ySvAwcXOp+jIGLge8udSfGgOMw4DgMOA4nHoNfrqo3z/XFTud7CEvp4HxunTrXJNnvODgOMxyHAcdh8cfAKSNJEmAgSJLauAfCrqXuwJhwHAYchwHHYcBxWOQxGOuLypKk0Rn3MwRJ0oiMbSAk2ZLkYJKpJDtPvcfZJcndSY4mOTBUW5NkMsmhXq7uepLc2WPxZJJrhvbZ3u0PJdk+23uNqyTrk3wpyTNJnk7ynq4vt3G4IMljSb7W4/B3Xb8iyaN9TJ9Ocn7XV/X2VD9/+dBrva/rB5PcuDRHtDBJViT5apLP9/ayG4ckzyV5KskTSfZ37cx/Lqpq7B7ACuBbwJXA+cDXgKuWul+LfIy/BVwDHBiq/QOws9d3An/f61uB/2Dwq3ObgEe7vgZ4tpere331Uh/bHMbgMuCaXn8j8E3gqmU4DgHe0OsrGfxE7SbgPuCWrn8U+LNe/3Pgo71+C/DpXr+qPyurgCv6M7RiqY9vHuPxF8Angc/39rIbB+A54OLjamf8czGuZwgbgamqeraqfgjcC2xb4j4tqqr6MnD870FsA3b3+m7g5qH6PTXwCHBR/0rdjcBkVR2rqheBSWDLme/94qiqI1X1lV5/GXiGwe9sL7dxqKr6v95c2Y8C3gbc3/Xjx2FmfO4HbugfstoG3FtVr1TVtxn8JsnGERzCokmyDvg94N96OyzDcTiBM/65GNdAWAs8P7Q93bVz3aU1+HU5enlJ1080HufMOPXp/lsZ/HW87Mahp0meYPBTtJMM/qr9XlW92k2Gj+mnx9vPvwS8iXNgHIB/Bv4K+Elvv4nlOQ4FfCHJ40l2dO2Mfy7G9ZvKmaW2nG+HOtF4nBPjlOQNDH6z+71V9f3BH3mzN52ldk6MQ1X9GLg6yUXA54C3zNasl+fkOCT5feBoVT2e5PqZ8ixNz+lxaNdV1eEklwCTSb5xkraLNg7jeoYwDawf2l4HHF6ivozSC32qRy+Pdv1E43HWj1OSlQzC4BNV9dkuL7txmFFV3wP+i8Fc8EVJZv5oGz6mnx5vP/9LDKYfz/ZxuA74gyTPMZgmfhuDM4blNg5U1eFeHmXwB8JGRvC5GNdA2Ads6LsLzmdwwWjPEvdpFPYAM3cCbAceGKrf2ncTbAJe6lPGh4DNSVb3HQebu3ZW6PnejwHPVNWHh55abuPw5j4zIMkvAr/D4HrKl4C3d7Pjx2FmfN4OfLEGVxH3ALf03TdXABuAx0ZzFAtXVe+rqnVVdTmDz/wXq+qPWGbjkOT1Sd44s87g3/MBRvG5WOqr6Se5yr6VwV0n3wLev9T9OQPH9yngCPAjBkl+O4P5z73AoV6u6bYBPtJj8RQwMfQ6f8LgotkUcNtSH9ccx+A3GZzCPgk80Y+ty3Acfg34ao/DAeBvun4lg//IpoB/B1Z1/YLenurnrxx6rff3+BwEblrqY1vAmFzPz+4yWlbj0Mf7tX48PfP/3yg+F35TWZIEjO+UkSRpxAwESRJgIEiSmoEgSQIMBElSMxAkSYCBIElqBoIkCYD/B3BOJvnIzBIZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make the list into an array (time will be dim=0)\n",
    "data_arr = np.array(data_ls)    \n",
    "\n",
    "# make an array of zeros to match\n",
    "shp = data_ls[0].shape\n",
    "targets = np.zeros((shp[0]*shp[1],1))\n",
    "\n",
    "# reshape\n",
    "data_arr = np.rollaxis(data_arr, 0, 3)\n",
    "data_vec = np.reshape(data_arr, (shp[0]*shp[1], -1))\n",
    "\n",
    "# make into a pandas data frame\n",
    "import pandas as pd\n",
    "        \n",
    "dts = []\n",
    "for i,t in enumerate(time_ls):\n",
    "    s = 'UTC_' + t.replace(':','').replace(' ','').replace('UTC', '')\n",
    "    dt = (s, 'float32')    \n",
    "    dts.append(dt)\n",
    "    \n",
    "# restructure the numpy array for use in a pandas dataframe\n",
    "myRecord = np.core.records.array(list(tuple(data_vec.transpose())), dtype=dts)\n",
    "df = pd.DataFrame(myRecord)\n",
    "\n",
    "\n",
    "absmax = data_vec.max(axis=1).max()\n",
    "loc = np.where(data_vec.max(axis=1) == absmax)[0]\n",
    "\n",
    "targets[loc]=1\n",
    "plt.imshow(np.reshape(targets, (shp[0], shp[1])))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADtZJREFUeJzt3H+M5Hddx/Hni54tQbC/7gq112NrekQPTAQnBeKvamm5ktAj2pirIRymeglaE0GNJcQUCn8ASmqIVTxp49lEWmyibERyKS0NhtDaOYrIobXL8aNrG3p4tUnTQD14+8d8S/azzt3O7Xxvp3v3fCSXne93Prvz/tzu9bnfmd2mqpAk6VnPm/UAkqTnFsMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEmNDbMeYDU2btxYc3Nzsx5DktaV/fv3f7uqNq20bl2GYW5ujuFwOOsxJGldSfKNSdb5VJIkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1eglDku1JHkqykOT6MfefkeSO7v77k8wtu39LkqeS/H4f80iSVm/qMCQ5DbgZuBLYBlyTZNuyZdcCT1TVxcBNwAeW3X8T8KlpZ5EkTa+PK4ZLgIWqOlhVzwC3AzuWrdkB7O1u3wlcliQASd4EHAQO9DCLJGlKfYThAuCRJceL3bmxa6rqCPAkcG6SHwb+EHhPD3NIknrQRxgy5lxNuOY9wE1V9dSKD5LsTjJMMjx06NAqxpQkTWJDDx9jEbhwyfFm4NGjrFlMsgE4EzgMvBq4OskHgbOA7yf5TlX92fIHqao9wB6AwWCwPDySpJ70EYYHgK1JLgL+C9gJ/NqyNfPALuDzwNXAPVVVwM89uyDJu4GnxkVBkrR2pg5DVR1Jch2wDzgNuLWqDiS5ERhW1TxwC3BbkgVGVwo7p31cSdKJkdE37uvLYDCo4XA46zEkaV1Jsr+qBiut8zefJUkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDV6CUOS7UkeSrKQ5Pox95+R5I7u/vuTzHXnL0+yP8m/dW9/qY95JEmrN3UYkpwG3AxcCWwDrkmybdmya4Enqupi4CbgA935bwNvrKqfBHYBt007jyRpOn1cMVwCLFTVwap6Brgd2LFszQ5gb3f7TuCyJKmqB6vq0e78AeD5Sc7oYSZJ0ir1EYYLgEeWHC9258auqaojwJPAucvW/ArwYFV9t4eZJEmrtKGHj5Ex5+p41iR5OaOnl6446oMku4HdAFu2bDn+KSVJE+njimERuHDJ8Wbg0aOtSbIBOBM43B1vBv4eeEtVffVoD1JVe6pqUFWDTZs29TC2JGmcPsLwALA1yUVJTgd2AvPL1swzenEZ4GrgnqqqJGcBnwTeWVWf62EWSdKUpg5D95rBdcA+4N+Bj1fVgSQ3JrmqW3YLcG6SBeAdwLM/0nodcDHwR0m+2P05b9qZJEmrl6rlLwc89w0GgxoOh7MeQ5LWlST7q2qw0jp/81mS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWr0EoYk25M8lGQhyfVj7j8jyR3d/fcnmVty3zu78w8leX0f80iSVm/qMCQ5DbgZuBLYBlyTZNuyZdcCT1TVxcBNwAe6990G7AReDmwH/rz7eJKkGenjiuESYKGqDlbVM8DtwI5la3YAe7vbdwKXJUl3/vaq+m5VfQ1Y6D6eJGlG+gjDBcAjS44Xu3Nj11TVEeBJ4NwJ31eStIb6CEPGnKsJ10zyvqMPkOxOMkwyPHTo0HGOKEmaVB9hWAQuXHK8GXj0aGuSbADOBA5P+L4AVNWeqhpU1WDTpk09jC1JGqePMDwAbE1yUZLTGb2YPL9szTywq7t9NXBPVVV3fmf3U0sXAVuBf+lhJknSKm2Y9gNU1ZEk1wH7gNOAW6vqQJIbgWFVzQO3ALclWWB0pbCze98DST4OfAU4Avx2VX1v2pkkSauX0Tfu68tgMKjhcDjrMSRpXUmyv6oGK63zN58lSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNaYKQ5JzktyV5OHu7dlHWberW/Nwkl3duRck+WSS/0hyIMn7p5lFktSPaa8YrgfurqqtwN3dcSPJOcANwKuBS4AblgTkT6rqx4FXAj+T5Mop55EkTWnaMOwA9na39wJvGrPm9cBdVXW4qp4A7gK2V9XTVfUZgKp6BvgCsHnKeSRJU5o2DC+uqscAurfnjVlzAfDIkuPF7twPJDkLeCOjqw5J0gxtWGlBkk8DLxlz17smfIyMOVdLPv4G4GPAh6vq4DHm2A3sBtiyZcuEDy1JOl4rhqGqXne0+5J8K8n5VfVYkvOBx8csWwQuXXK8Gbh3yfEe4OGq+tMV5tjTrWUwGNSx1kqSVm/ap5LmgV3d7V3AJ8as2QdckeTs7kXnK7pzJHkfcCbwu1POIUnqybRheD9weZKHgcu7Y5IMknwUoKoOA+8FHuj+3FhVh5NsZvR01DbgC0m+mOQ3ppxHkjSlVK2/Z2UGg0ENh8NZjyFJ60qS/VU1WGmdv/ksSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1pgpDknOS3JXk4e7t2UdZt6tb83CSXWPun0/y5WlmkST1Y9orhuuBu6tqK3B3d9xIcg5wA/Bq4BLghqUBSfLLwFNTziFJ6sm0YdgB7O1u7wXeNGbN64G7qupwVT0B3AVsB0jyQuAdwPumnEOS1JNpw/DiqnoMoHt73pg1FwCPLDle7M4BvBf4EPD0lHNIknqyYaUFST4NvGTMXe+a8DEy5lwl+Sng4qp6e5K5CebYDewG2LJly4QPLUk6XiuGoaped7T7knwryflV9ViS84HHxyxbBC5dcrwZuBd4LfDTSb7ezXFeknur6lLGqKo9wB6AwWBQK80tSVqdaZ9Kmgee/SmjXcAnxqzZB1yR5OzuRecrgH1V9RdV9aNVNQf8LPCfR4uCJGntTBuG9wOXJ3kYuLw7JskgyUcBquowo9cSHuj+3NidkyQ9B6Vq/T0rMxgMajgcznoMSVpXkuyvqsFK6/zNZ0lSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSI1U16xmOW5JDwDdW+e4bgW/3OM564J5PDafank+1/cL0e35pVW1aadG6DMM0kgyrajDrOdaSez41nGp7PtX2C2u3Z59KkiQ1DIMkqXEqhmHPrAeYAfd8ajjV9nyq7RfWaM+n3GsMkqRjOxWvGCRJx3DShiHJ9iQPJVlIcv2Y+89Ickd3//1J5tZ+yv5MsN93JPlKki8luTvJS2cxZ59W2vOSdVcnqSTr/idYJtlzkl/tPtcHkvztWs/Ytwm+trck+UySB7uv7zfMYs6+JLk1yeNJvnyU+5Pkw93fx5eSvKr3IarqpPsDnAZ8Ffgx4HTgX4Fty9b8FvCR7vZO4I5Zz32C9/uLwAu6229bz/uddM/duhcBnwXuAwaznnsNPs9bgQeBs7vj82Y99xrseQ/wtu72NuDrs557yj3/PPAq4MtHuf8NwKeAAK8B7u97hpP1iuESYKGqDlbVM8DtwI5la3YAe7vbdwKXJckaztinFfdbVZ+pqqe7w/uAzWs8Y98m+RwDvBf4IPCdtRzuBJlkz78J3FxVTwBU1eNrPGPfJtlzAT/S3T4TeHQN5+tdVX0WOHyMJTuAv6mR+4Czkpzf5wwnaxguAB5ZcrzYnRu7pqqOAE8C567JdP2bZL9LXcvoO471bMU9J3klcGFV/eNaDnYCTfJ5fhnwsiSfS3Jfku1rNt2JMcme3w28Ocki8E/A76zNaDNzvP/ej9uGPj/Yc8i47/yX//jVJGvWi4n3kuTNwAD4hRM60Yl3zD0neR5wE/DWtRpoDUzyed7A6OmkSxldFf5zkldU1f+c4NlOlEn2fA3w11X1oSSvBW7r9vz9Ez/eTJzw/3adrFcMi8CFS4438/8vL3+wJskGRpegx7p8ey6bZL8keR3wLuCqqvruGs12oqy05xcBrwDuTfJ1Rs/Fzq/zF6An/br+RFX9b1V9DXiIUSjWq0n2fC3wcYCq+jzwfEb/T6GT1UT/3qdxsobhAWBrkouSnM7oxeX5ZWvmgV3d7auBe6p7ZWcdWnG/3dMqf8koCuv9eWdYYc9V9WRVbayquaqaY/S6ylVVNZzNuL2Y5Ov6Hxj9oAFJNjJ6aungmk7Zr0n2/E3gMoAkP8EoDIfWdMq1NQ+8pfvppNcAT1bVY30+wEn5VFJVHUlyHbCP0U813FpVB5LcCAyrah64hdEl5wKjK4Wds5t4OhPu94+BFwJ/173G/s2qumpmQ09pwj2fVCbc8z7giiRfAb4H/EFV/ffspp7OhHv+PeCvkryd0VMqb13H3+SR5GOMngrc2L1ucgPwQwBV9RFGr6O8AVgAngZ+vfcZ1vHfnyTpBDhZn0qSJK2SYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLU+D+elOgAosSSTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "arr_vars = np.var(data_vec, axis=1)\n",
    "hist, bins, _ = plt.hist(arr_vars[arr_vars>0], bins=255)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bins[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABI4AAAH8CAYAAAC3ueT0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHqVJREFUeJzt3W3M5fVd5/HPV2i7rroLtTdBwAXdqZFudrGdUJJGU+8KZY3Ujd2l2VjSJRk1bdImJit1H+DWfeDuWrtpVjGjktJNBVm1lphqHbFrn0gL0yI3pcjQYpkygSi1ranBgN99cP5Tj/Cdmatz3cwFvF7JyTnnd37nnP958IOLN/+b6u4AAAAAwFN93aneAAAAAAB2J+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABjteDiqqkur6r6qOlRVV+/09wMAAACwMdXdO/dlVacl+fMkP5jkcJLbkryxuz+1YxsBAAAAwIbs9B5HFyU51N2f6e6/S3Jjkst3eBsAAAAA2ICdDkdnJ3lo7fnhZQwAAACAXeb0Hf6+GsaedqxcVe1Lsm95+spt3SIAAACA556/7O4Xn2jSToejw0nOXXt+TpKHnzqpu/cn2Z8kVbVzJ2ECAAAAeG74i41M2ulD1W5Lsqeqzq+q5ye5IsnNO7wNAAAAAGzAju5x1N1PVNVbk3w4yWlJruvue3ZyGwAAAADYmOre3UeCOVQNAAAAYMsd7O69J5q004eqAQAAAPAMIRwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAo5MOR1V1blV9pKrurap7qupty/jPVtXnq+qO5XbZ2nveUVWHquq+qrpkK34AAAAAANvj9E2894kkP9Xdn6iqb0pysKoOLK+9u7t/YX1yVV2Q5IokL0/yLUn+qKpe1t1PbmIbAAAAANgmJ73HUXcf6e5PLI+/nOTeJGcf5y2XJ7mxux/v7s8mOZTkopP9fgAAAAC215ac46iqzkvyXUk+tgy9tarurKrrqurMZezsJA+tve1wjh+aAAAAADiFNh2Oquobk/x2krd395eSXJvk25NcmORIkncdnTq8vY/xmfuq6vaqun2z2wcAAADAydlUOKqq52UVjd7f3b+TJN39SHc/2d1/n+RX8w+Hox1Ocu7a289J8vD0ud29v7v3dvfezWwfAAAAACdvM1dVqyS/nuTe7v7FtfGz1qb9SJK7l8c3J7miql5QVecn2ZPk4yf7/QAAAABsr81cVe3VSX4syV1Vdccy9jNJ3lhVF2Z1GNqDSX48Sbr7nqq6Kcmnsroi21tcUQ0AAABg96ru8TRDu0ZV7e4NBAAAAHjmObiRUwRtyVXVAAAAAHj2EY4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACA0abDUVU9WFV3VdUdVXX7MvbCqjpQVfcv92cu41VV76mqQ1V1Z1W9YrPfDwAAAMD22Ko9jr63uy/s7r3L86uT3NLde5LcsjxPktcl2bPc9iW5dou+HwAAAIAttl2Hql2e5Prl8fVJXr82/r5euTXJGVV11jZtAwAAAACbsBXhqJP8YVUdrKp9y9hLu/tIkiz3L1nGz07y0Np7Dy9jAAAAAOwyp2/BZ7y6ux+uqpckOVBVnz7O3BrG+mmTVgFq3zAXAAAAgB2y6T2Ouvvh5f7RJB9IclGSR44egrbcP7pMP5zk3LW3n5Pk4eEz93f33rVzJgEAAACwwzYVjqrqG6rqm44+TvLaJHcnuTnJlcu0K5N8cHl8c5I3LVdXuzjJF48e0gYAAADA7rLZQ9VemuQDVXX0s36ju/+gqm5LclNVXZXkc0nesMz/UJLLkhxK8pUkb97k9wMAAACwTar7aacY2lWqandvIAAAAMAzz8GNnCJoK66qBgAAAMCzkHAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjE46HFXVd1TVHWu3L1XV26vqZ6vq82vjl6295x1Vdaiq7quqS7bmJwAAAACwHaq7N/8hVacl+XySVyV5c5K/6e5feMqcC5LckOSiJN+S5I+SvKy7nzzBZ29+AwEAAABYd7C7955o0lYdqvb9SR7o7r84zpzLk9zY3Y9392eTHMoqIgEAAACwC21VOLoiq72JjnprVd1ZVddV1ZnL2NlJHlqbc3gZe5qq2ldVt1fV7Vu0fQAAAAB8jTYdjqrq+Ul+OMn/XYauTfLtSS5MciTJu45OHd4+HobW3fu7e+9GdpkCAAAAYHtsxR5Hr0vyie5+JEm6+5HufrK7/z7Jr+YfDkc7nOTctfedk+ThLfh+AAAAALbBVoSjN2btMLWqOmvttR9Jcvfy+OYkV1TVC6rq/CR7knx8C74fAAAAgG1w+mbeXFX/NMkPJvnxteH/UVUXZnUY2oNHX+vue6rqpiSfSvJEkrec6IpqAAAAAJw61b27r3ZfVbt7AwEAAACeeQ5u5NzSW3VVNQAAAACeZYQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADAaEPhqKquq6pHq+rutbEXVtWBqrp/uT9zGa+qek9VHaqqO6vqFWvvuXKZf39VXbn1PwcAAACArbLRPY7em+TSp4xdneSW7t6T5JbleZK8Lsme5bYvybXJKjQluSbJq5JclOSao7EJAAAAgN1nQ+Gouz+a5LGnDF+e5Prl8fVJXr82/r5euTXJGVV1VpJLkhzo7se6+wtJDuTpMQoAAACAXWIz5zh6aXcfSZLl/iXL+NlJHlqbd3gZO9Y4AAAAALvQ6dvwmTWM9XHGn/4BVfuyOswNAAAAgFNkM3scPbIcgpbl/tFl/HCSc9fmnZPk4eOMP0137+/uvd29dxPbBwAAAMAmbCYc3Zzk6JXRrkzywbXxNy1XV7s4yReXQ9k+nOS1VXXmclLs1y5jAAAAAOxCGzpUrapuSPKaJC+qqsNZXR3t55PcVFVXJflckjcs0z+U5LIkh5J8Jcmbk6S7H6uqn0ty2zLvnd391BNuAwAAALBLVPd4mqFdo6p29wYCAAAAPPMc3MgpgjZzqBoAAAAAz2LCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADA6YTiqquuq6tGquntt7H9W1aer6s6q+kBVnbGMn1dVf1tVdyy3X1l7zyur6q6qOlRV76mq2p6fBAAAAMBW2MgeR+9NculTxg4k+Vfd/a+T/HmSd6y99kB3X7jcfmJt/Nok+5LsWW5P/UwAAAAAdpEThqPu/miSx54y9ofd/cTy9NYk5xzvM6rqrCT/rLv/tLs7yfuSvP7kNhkAAACAnbAV5zj6T0l+f+35+VX1yar6k6r67mXs7CSH1+YcXsYAAAAA2KVO38ybq+q/JHkiyfuXoSNJvrW7/6qqXpnkd6vq5Umm8xn1cT53X1aHtQEAAABwipx0OKqqK5P8UJLvXw4/S3c/nuTx5fHBqnogycuy2sNo/XC2c5I8fKzP7u79SfYv33PMwAQAAADA9jmpQ9Wq6tIkP53kh7v7K2vjL66q05bH35bVSbA/091Hkny5qi5erqb2piQf3PTWAwAAALBtTrjHUVXdkOQ1SV5UVYeTXJPVVdRekOTAqgPl1uUKat+T5J1V9USSJ5P8RHcfPbH2T2Z1hbavz+qcSOvnRQIAAABgl6nlKLNdy6FqAAAAAFvuYHfvPdGkrbiqGgAAAADPQsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMDphOKqq66rq0aq6e23sZ6vq81V1x3K7bO21d1TVoaq6r6ouWRu/dBk7VFVXb/1PAQAAAGArbWSPo/cmuXQYf3d3X7jcPpQkVXVBkiuSvHx5zy9X1WlVdVqSX0ryuiQXJHnjMhcAAACAXer0E03o7o9W1Xkb/LzLk9zY3Y8n+WxVHUpy0fLaoe7+TJJU1Y3L3E99zVsMAAAAwI7YzDmO3lpVdy6Hsp25jJ2d5KG1OYeXsWONAwAAALBLnWw4ujbJtye5MMmRJO9axmuY28cZH1XVvqq6vapuP8ntAwAAAGCTTnio2qS7Hzn6uKp+NcnvLU8PJzl3beo5SR5eHh9rfPr8/Un2L59/zMAEAAAAwPY5qT2Oquqstac/kuToFdduTnJFVb2gqs5PsifJx5PclmRPVZ1fVc/P6gTaN5/8ZgMAAACw3U64x1FV3ZDkNUleVFWHk1yT5DVVdWFWh5s9mOTHk6S776mqm7I66fUTSd7S3U8un/PWJB9OclqS67r7ni3/NQAAAABsmere3UeCOVQNAAAAYMsd7O69J5q0mauqAQAAAPAsJhwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAoxOGo6q6rqoeraq718Z+s6ruWG4PVtUdy/h5VfW3a6/9ytp7XllVd1XVoap6T1XV9vwkAAAAALbC6RuY894k/zvJ+44OdPd/OPq4qt6V5Itr8x/o7guHz7k2yb4ktyb5UJJLk/z+177JAAAAAOyEE+5x1N0fTfLY9Nqy19C/T3LD8T6jqs5K8s+6+0+7u7OKUK//2jcXAAAAgJ2y2XMcfXeSR7r7/rWx86vqk1X1J1X13cvY2UkOr805vIwBAAAAsEtt5FC143lj/vHeRkeSfGt3/1VVvTLJ71bVy5NM5zPqY31oVe3L6rA2AAAAAE6Rkw5HVXV6kn+X5JVHx7r78SSPL48PVtUDSV6W1R5G56y9/ZwkDx/rs7t7f5L9y/ccMzABAAAAsH02c6jaDyT5dHd/9RC0qnpxVZ22PP62JHuSfKa7jyT5clVdvJwX6U1JPriJ7wYAAABgm50wHFXVDUn+NMl3VNXhqrpqeemKPP2k2N+T5M6q+rMkv5XkJ7r76Im1fzLJryU5lOSBuKIaAAAAwK5Wq4uc7V4OVQMAAADYcge7e++JJm32qmoAAAAAPEsJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMDohOGoqs6tqo9U1b1VdU9VvW0Zf2FVHaiq+5f7M5fxqqr3VNWhqrqzql6x9llXLvPvr6ort+9nAQAAALBZG9nj6IkkP9Xd35nk4iRvqaoLklyd5Jbu3pPkluV5krwuyZ7lti/JtckqNCW5JsmrklyU5JqjsQkAAACA3eeE4ai7j3T3J5bHX05yb5Kzk1ye5Ppl2vVJXr88vjzJ+3rl1iRnVNVZSS5JcqC7H+vuLyQ5kOTSLf01AAAAAGyZ07+WyVV1XpLvSvKxJC/t7iPJKi5V1UuWaWcneWjtbYeXsWONT9+zL6u9lQAAAAA4RTYcjqrqG5P8dpK3d/eXquqYU4exPs740we79yfZv3zvOAcAAACA7bWhq6pV1fOyikbv7+7fWYYfWQ5By3L/6DJ+OMm5a28/J8nDxxkHAAAAYBfayFXVKsmvJ7m3u39x7aWbkxy9MtqVST64Nv6m5epqFyf54nJI24eTvLaqzlxOiv3aZQwAAACAXWgjh6q9OsmPJbmrqu5Yxn4myc8nuamqrkryuSRvWF77UJLLkhxK8pUkb06S7n6sqn4uyW3LvHd292Nb8isAAAAA2HLVvbtPIeQcRwAAAABb7mB37z3RpA2d4wgAAACA5x7hCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADASjgAAAAAYCUcAAAAAjIQjAAAAAEbCEQAAAAAj4QgAAACAkXAEAAAAwEg4AgAAAGAkHAEAAAAwEo4AAAAAGAlHAAAAAIyEIwAAAABGwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACMhCMAAAAARsIRAAAAACPhCAAAAICRcAQAAADASDgCAAAAYCQcAQAAADA6/VRvwAb8TZL7TvVGwHPAi5L85aneCHgOsNZg51hvsDOsNdgZW73W/sVGJj0TwtF93b33VG8EPNtV1e3WGmw/aw12jvUGO8Nag51xqtaaQ9UAAAAAGAlHAAAAAIyeCeFo/6neAHiOsNZgZ1hrsHOsN9gZ1hrsjFOy1qq7T8X3AgAAALDLPRP2OAIAAADgFNi14aiqLq2q+6rqUFVdfaq3B56Jquq6qnq0qu5eG3thVR2oqvuX+zOX8aqq9yxr7s6qesXae65c5t9fVVeeit8Cu1VVnVtVH6mqe6vqnqp62zJurcEWq6p/UlUfr6o/W9bbf13Gz6+qjy1r5zer6vnL+AuW54eW189b+6x3LOP3VdUlp+YXwe5WVadV1Ser6veW59YabLGqerCq7qqqO6rq9mVsV/0duSvDUVWdluSXkrwuyQVJ3lhVF5zarYJnpPcmufQpY1cnuaW79yS5ZXmerNbbnuW2L8m1yeofWkmuSfKqJBclueboP7iAJMkTSX6qu78zycVJ3rL8O8tag633eJLv6+5/k+TCJJdW1cVJ/nuSdy/r7QtJrlrmX5XkC939L5O8e5mXZY1ekeTlWf178peXvz+Bf+xtSe5de26twfb43u6+sLv3Ls931d+RuzIcZfVDD3X3Z7r775LcmOTyU7xN8IzT3R9N8thThi9Pcv3y+Pokr18bf1+v3JrkjKo6K8klSQ5092Pd/YUkB/L0GAXPWd19pLs/sTz+clZ/YJ8daw223LJu/mZ5+rzl1km+L8lvLeNPXW9H1+FvJfn+qqpl/Mbufry7P5vkUFZ/fwKLqjonyb9N8mvL84q1BjtlV/0duVvD0dlJHlp7fngZAzbvpd19JFn9B2+Slyzjx1p31iNs0LJr/ncl+VisNdgWy6EzdyR5NKs/jB9I8tfd/cQyZX3tfHVdLa9/Mck3x3qDjfhfSf5zkr9fnn9zrDXYDp3kD6vqYFXtW8Z21d+Rp2/VB22xGsZc/g2217HWnfUIG1BV35jkt5O8vbu/tPofrfPUYcxagw3q7ieTXFhVZyT5QJLvnKYt99YbnISq+qEkj3b3wap6zdHhYaq1Bpv36u5+uKpekuRAVX36OHNPyVrbrXscHU5y7trzc5I8fIq2BZ5tHll2Z8xy/+gyfqx1Zz3CCVTV87KKRu/v7t9Zhq012Ebd/ddJ/l9W5xY7o6qO/g/R9bXz1XW1vP7PszqE23qD43t1kh+uqgezOm3I92W1B5K1Blusux9e7h/N6n+IXJRd9nfkbg1HtyXZs5y1//lZnVDt5lO8TfBscXOSo2fZvzLJB9fG37Scqf/iJF9cdov8cJLXVtWZywnWXruMAfnqOR9+Pcm93f2Lay9Za7DFqurFy55GqaqvT/IDWZ1X7CNJfnSZ9tT1dnQd/miSP+7uXsavWK4EdX5WJxn9+M78Ctj9uvsd3X1Od5+X1X+L/XF3/8dYa7Clquobquqbjj7O6u+/u7PL/o7clYeqdfcTVfXWrH7oaUmu6+57TvFmwTNOVd2Q5DVJXlRVh7M60/7PJ7mpqq5K8rkkb1imfyjJZVmdtPArSd6cJN39WFX9XFZBN0ne2d1PPeE2PJe9OsmPJblrOe9KkvxMrDXYDmcluX65KtPXJbmpu3+vqj6V5Maq+m9JPplVzM1y/3+q6lBWez9ckSTdfU9V3ZTkU1ldGfEtyyFwwPH9dKw12EovTfKB5RQHpyf5je7+g6q6Lbvo78hahWAAAAAA+Md266FqAAAAAJxiwhEAAAAAI+EIAAAAgJFwBAAAAMBIOAIAAABgJBwBAAAAMBKOAAAAABgJRwAAAACM/j+LqmOaLWaT6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x1440 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vars_hi_im = np.reshape(arr_vars>bins[-254], (shp[0], shp[1]))\n",
    "vars_im = np.reshape(arr_vars, (shp[0], shp[1]))\n",
    "fig = plt.figure(figsize=(20,20))\n",
    "plt.imshow(vars_im, cmap='bone', vmin=0., vmax=80)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Write out one of the rasters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal, osr\n",
    "driver = gdal.GetDriverByName('GTiff')\n",
    "rows,cols = data_ls[0].shape\n",
    "\n",
    "fname = dts[0][0] + \"_\" + date_ls[0].replace('-','') + '_GOES_L1B_RAD_fromJeff.tif'\n",
    "dst_ds = driver.Create(fname, cols, rows, 1, gdal.GDT_Float32)\n",
    "dst_ds.SetGeoTransform(grid_ls[0].GetGeoTransform())\n",
    "#srs = osr.SpatialReference()\n",
    "#srs.ImportFromWkt(projInfo)\n",
    "dst_ds.SetProjection(grid_ls[0].GetProjection()) \n",
    "\n",
    "# Write out_arr to band\n",
    "dst_ds.GetRasterBand(1).WriteArray(data_ls[0])\n",
    "\n",
    "# Flush to disk and close file\n",
    "dst_ds.FlushCache()\n",
    "dst_ds = None\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1500, 2500), (2157, 5089))"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid2 = grid.ReadAsArray()\n",
    "raw.shape, grid.ReadAsArray().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "31.9973146667\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEIdJREFUeJzt3X2sZVV9xvHvU176gqSAXBGBcawlpGhkJDejhtbgCwgjEW3UMmkttphRA42m/UNqE7E2TeiL2rQYyCgTsFHQVlESRmVibZDEF+7QAQYRQTKW6xDmIgoSbc3or3/cPen1cs6dy9ln5t7D+n6Sk7P32uvstVZ25pk96+yzJlWFJKkdv7LSHZAkHVwGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4Jakxh650BwY59thja+3atSvdDUmaGNu3b3+kqqaWU3dVBv/atWuZmZlZ6W5I0sRI8r3l1nWqR5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGrMqf7nbx9pLbxpYvuvy1x7knkjS6uQdvyQ1Zr93/Em2AOcBe6rqhV3Zp4BTuipHAT+qqnUDPrsL+DHwc2BvVU2Pqd+SpBEtZ6rnGuAK4OP7CqrqD/ZtJ/kg8NgSn39FVT0yagclSeO13+CvqluSrB10LEmANwOvHG+3JEkHSt85/t8DHq6q+4YcL+DmJNuTbFrqREk2JZlJMjM3N9ezW5KkYfoG/0bguiWOn1FVpwPnAhcnefmwilW1uaqmq2p6ampZ/5eAJGkEIwd/kkOB3wc+NaxOVe3u3vcANwDrR21PkjQefe74Xw18u6pmBx1MckSSI/dtA2cDO3u0J0kag/0Gf5LrgK8BpySZTXJRd+gCFk3zJHlOkq3d7nHArUnuAL4J3FRVXxxf1yVJo1jOUz0bh5S/dUDZbmBDt/0AcFrP/kmSxsxf7kpSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTH7Df4kW5LsSbJzQdn7k3w/yY7utWHIZ89Jcm+S+5NcOs6OS5JGs5w7/muAcwaUf7iq1nWvrYsPJjkE+AhwLnAqsDHJqX06K0nqb7/BX1W3AI+OcO71wP1V9UBV/Qy4Hjh/hPNIksaozxz/JUnu7KaCjh5w/ATgwQX7s13ZQEk2JZlJMjM3N9ejW5KkpYwa/FcCzwfWAQ8BHxxQJwPKatgJq2pzVU1X1fTU1NSI3ZIk7c9IwV9VD1fVz6vqF8BHmZ/WWWwWOGnB/onA7lHakySNz0jBn+T4BbtvAHYOqHYbcHKS5yU5HLgAuHGU9iRJ43Po/iokuQ44Ezg2ySxwGXBmknXMT93sAt7e1X0O8LGq2lBVe5NcAnwJOATYUlV3H5BRSJKWbb/BX1UbBxRfPaTubmDDgv2twJMe9ZQkrRx/uStJjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1Zr/Bn2RLkj1Jdi4o+4ck305yZ5Ibkhw15LO7ktyVZEeSmXF2XJI0muXc8V8DnLOobBvwwqp6EfAd4C+X+PwrqmpdVU2P1kVJ0jjtN/ir6hbg0UVlN1fV3m7368CJB6BvkqQDYBxz/H8KfGHIsQJuTrI9yaYxtCVJ6unQPh9O8lfAXuATQ6qcUVW7kzwL2Jbk292/IAadaxOwCWDNmjV9uiVJWsLId/xJLgTOA/6wqmpQnara3b3vAW4A1g87X1VtrqrpqpqempoatVuSpP0YKfiTnAO8B3hdVf1kSJ0jkhy5bxs4G9g5qK4k6eBZzuOc1wFfA05JMpvkIuAK4Ejmp292JLmqq/ucJFu7jx4H3JrkDuCbwE1V9cUDMgpJ0rLtd46/qjYOKL56SN3dwIZu+wHgtF69kySNnb/claTGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSY5YV/Em2JNmTZOeCsmOSbEtyX/d+9JDPXtjVuS/JhePquCRpNMu9478GOGdR2aXAl6vqZODL3f4vSXIMcBnwEmA9cNmwvyAkSQfHsoK/qm4BHl1UfD5wbbd9LfD6AR99DbCtqh6tqh8C23jyXyCSpIOozxz/cVX1EED3/qwBdU4AHlywP9uVSZJWyIH+cjcDympgxWRTkpkkM3Nzcwe4W5LUrj7B/3CS4wG69z0D6swCJy3YPxHYPehkVbW5qqaranpqaqpHtyRJS+kT/DcC+57SuRD4/IA6XwLOTnJ096Xu2V2ZJGmFLPdxzuuArwGnJJlNchFwOXBWkvuAs7p9kkwn+RhAVT0K/A1wW/f6QFcmSVohhy6nUlVtHHLoVQPqzgBvW7C/BdgyUu8kSWPnL3clqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktSYkYM/ySlJdix4PZ7k3YvqnJnksQV13te/y5KkPg4d9YNVdS+wDiDJIcD3gRsGVP1qVZ03ajuSpPEa11TPq4DvVtX3xnQ+SdIBMq7gvwC4bsixlyW5I8kXkrxgTO1JkkbUO/iTHA68Dvi3AYdvB55bVacB/wJ8bonzbEoyk2Rmbm6ub7ckSUOM447/XOD2qnp48YGqeryqnui2twKHJTl20EmqanNVTVfV9NTU1Bi6JUkaZBzBv5Eh0zxJnp0k3fb6rr0fjKFNSdKIRn6qByDJbwBnAW9fUPYOgKq6Cngj8M4ke4GfAhdUVfVpU5LUT6/gr6qfAM9cVHbVgu0rgCv6tCFJGq9ewd+itZfeNLB81+WvPcg9kTSpVjpHXLJBkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjGv1jMlKr70hScvlHb8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY3pHfxJdiW5K8mOJDMDjifJPye5P8mdSU7v26YkaXTj+gHXK6rqkSHHzgVO7l4vAa7s3iVJK+BgTPWcD3y85n0dOCrJ8QehXUnSAOMI/gJuTrI9yaYBx08AHlywP9uV/ZIkm5LMJJmZm5sbQ7ckSYOMI/jPqKrTmZ/SuTjJyxcdz4DP1JMKqjZX1XRVTU9NTY2hW5KkQXoHf1Xt7t73ADcA6xdVmQVOWrB/IrC7b7uSpNH0Cv4kRyQ5ct82cDawc1G1G4E/7p7ueSnwWFU91KddSdLo+j7VcxxwQ5J95/pkVX0xyTsAquoqYCuwAbgf+AnwJz3blCT10Cv4q+oB4LQB5Vct2C7g4j7tSJLGx1/uSlJjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmP6/teLE2/tpTcNLN91+Wsnvt2VGtu4DOs/DB/DpI9ZOhi845ekxowc/ElOSvKVJPckuTvJuwbUOTPJY0l2dK/39euuJKmvPlM9e4G/qKrbkxwJbE+yraq+tajeV6vqvB7tSJLGaOQ7/qp6qKpu77Z/DNwDnDCujkmSDoyxzPEnWQu8GPjGgMMvS3JHki8kecE42pMkja73Uz1JngF8Bnh3VT2+6PDtwHOr6okkG4DPAScPOc8mYBPAmjVr+nZLkjRErzv+JIcxH/qfqKrPLj5eVY9X1RPd9lbgsCTHDjpXVW2uqumqmp6amurTLUnSEvo81RPgauCeqvrQkDrP7uqRZH3X3g9GbVOS1F+fqZ4zgLcAdyXZ0ZW9F1gDUFVXAW8E3plkL/BT4IKqqh5tSpJ6Gjn4q+pWIPupcwVwxahtSJLGz1/uSlJjmlmrZ6l1X6RJNcp6RpJ3/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5Ia08ySDRrdsGUBXBJAmkze8UtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1JhewZ/knCT3Jrk/yaUDjv9qkk91x7+RZG2f9iRJ/Y0c/EkOAT4CnAucCmxMcuqiahcBP6yq3wY+DPzdqO1Jksajzx3/euD+qnqgqn4GXA+cv6jO+cC13fa/A69Kkh5tSpJ66hP8JwAPLtif7coG1qmqvcBjwDN7tClJ6ilVNdoHkzcBr6mqt3X7bwHWV9WfLahzd1dnttv/blfnBwPOtwnY1O2eAtzbbR8LPDJSJ1cvxzQZHNNkcEzznltVU8up2GeRtlngpAX7JwK7h9SZTXIo8JvAo4NOVlWbgc2Ly5PMVNV0j36uOo5pMjimyeCYnro+Uz23AScneV6Sw4ELgBsX1bkRuLDbfiPwHzXqPzEkSWMx8h1/Ve1NcgnwJeAQYEtV3Z3kA8BMVd0IXA38a5L7mb/Tv2AcnZYkja7XevxVtRXYuqjsfQu2/wd4U582GDD98zTgmCaDY5oMjukpGvnLXUnSZHLJBklqzKoO/v0tCTGJkuxKcleSHUlmVro/o0iyJcmeJDsXlB2TZFuS+7r3o1eyj0/VkDG9P8n3u2u1I8mGlezjU5XkpCRfSXJPkruTvKsrn9hrtcSYJvZaJfm1JN9Mckc3pr/uyp/XLXVzX7f0zeFja3O1TvV0S0J8BziL+cdCbwM2VtW3VrRjPSXZBUxX1cQ+d5zk5cATwMer6oVd2d8Dj1bV5d1f0kdX1XtWsp9PxZAxvR94oqr+cSX7NqokxwPHV9XtSY4EtgOvB97KhF6rJcb0Zib0WnWrGRxRVU8kOQy4FXgX8OfAZ6vq+iRXAXdU1ZXjaHM13/EvZ0kIrYCquoUn/x5j4fIc1zL/h3FiDBnTRKuqh6rq9m77x8A9zP+afmKv1RJjmlg174lu97DuVcArmV/qBsZ8nVZz8C9nSYhJVMDNSbZ3v1Z+ujiuqh6C+T+cwLNWuD/jckmSO7upoImZElmsWxn3xcA3eJpcq0Vjggm+VkkOSbID2ANsA74L/Khb6gbGnH+rOfgHLea2Ouelnpozqup05lc1vbibYtDqdCXwfGAd8BDwwZXtzmiSPAP4DPDuqnp8pfszDgPGNNHXqqp+XlXrmF8BYT3wO4Oqjau91Rz8y1kSYuJU1e7ufQ9wA/MX+eng4W7+dd887J4V7k9vVfVw9wfyF8BHmcBr1c0Zfwb4RFV9tiue6Gs1aExPh2sFUFU/Av4TeClwVLfUDYw5/1Zz8C9nSYiJkuSI7gspkhwBnA3sXPpTE2Ph8hwXAp9fwb6Mxb5w7LyBCbtW3ZeGVwP3VNWHFhya2Gs1bEyTfK2STCU5qtv+deDVzH938RXml7qBMV+nVftUD0D3SNY/8f9LQvztCneplyS/xfxdPsz/avqTkzimJNcBZzK/guDDwGXA54BPA2uA/wbeVFUT82XpkDGdyfzUQQG7gLfvmxufBEl+F/gqcBfwi674vczPiU/ktVpiTBuZ0GuV5EXMf3l7CPM345+uqg90eXE9cAzwX8AfVdX/jqXN1Rz8kqTxW81TPZKkA8Dgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMf8Hn0v5iIIE5ksAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(raw.mean())\n",
    "plt.hist(grid2[grid2>1].ravel(), bins=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x2bc8e998a58>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAADxCAYAAAAEJzaTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFjFJREFUeJzt3X+s1fV9x/HnS/xVbR0i1TCwk66ss20qtUTZzBqrosiaYhddcD8khoxlwa1dl6x2yeLarosm3dzMnMldZeLSopTWSDqiRWpnmkwFqiJICbfW6i0MiqDV2SLc+94f38/Br5dz7v2ec8+vz+X1SL459/s53/s9n489fd8P7+/nhyICMzPLywm9roCZmTXPwdvMLEMO3mZmGXLwNjPLkIO3mVmGHLzNzDLk4G1m1gGSVkraJ2lbqWyapA2SdqXXM1O5JN0haVDSVkkXjnd/B28zs864B1g4quxmYGNEzAE2pnOAq4E56VgO3DXezR28zcw6ICIeAw6MKl4MrEo/rwKuKZXfG4XHgamSZox1/xPbWVkzs5xd9fHT4+UDw5Wu3bL10Hbgl6WigYgYGOfXzomIPQARsUfS2al8JvBS6bqhVLan0Y0cvM3MkpcPDPPkw++pdO2UGbt+GRHz2vTRqlM25tolDt5mZkkAI4x08iP2SpqRet0zgH2pfAg4t3TdLGD3WDdyztvMLAmCwzFc6WjROmBp+nkp8GCp/IY06mQ+8GotvdKIe95mZiXt6nlLWg1cCkyXNATcAtwKrJG0DHgRuC5dvh5YBAwCbwA3jnd/B28zsyQIhtu0THZEXN/grcvrXBvAimbu7+BtZlYyMvZzwr7h4G1mlgQw7OBtZpYf97zNzDITwOFMtoZ08DYzS4Jw2sTMLDsBw3nEbgdvM7OaYoZlHhy8zcyOEsN1lxnpPw7eZmZJ8cDSwdvMLCvFOG8HbzOz7Iy4521mlhf3vM3MMhSI4UxWynbwNjMrcdrEzCwzgXgzpvS6GpU4eJuZJcUkHadNzMyy4weWZmaZiRDD4Z63mVl2RtzzNjPLS/HAMo+wmEctzcy6wA8szcwyNZzJOO88/sSYmXVBbYZllaMKSZ+WtE3SdkmfSWXTJG2QtCu9ntlKXR28zcxKRuKESsd4JH0I+BPgIuAC4BOS5gA3AxsjYg6wMZ03zcHbzCwpFqZqW8/7fODxiHgjIo4A/w18ClgMrErXrAKuaaWuznmbmSWBOFx9evx0SZtL5wMRMVA63wZ8WdJZwC+ARcBm4JyI2AMQEXsknd1KXR28zcySCJqZpLM/IuY1vlfskHQbsAF4HXgGODLxWhacNjEzO0qMVDyqiIi7I+LCiPgYcADYBeyVNAMgve5rpaYO3mZmSVD0vKscVdRSIpLeA/wesBpYByxNlywFHmylrk6bmJmVtHkzhm+mnPdhYEVEHJR0K7BG0jLgReC6Vm7s4G1mlgRq62YMEfE7dcpeBi6f6L0dvM3MkgAOe20TM7PcyOt5m5nlJqDS7Ml+4OBtZlbinreZWWYi5J63mVluigeW3j3ezCwz3sPSzCw7xQNL57zNzLLT5hmWHePgbWaWtHuGZSc5eJuZlXgDYjOzzETA4ZE8gnfXaylpoaSdkgYltbR3m5lZJxRpk/bsYdlpXe15S5oC3AksAIaATZLWRcRz3ayHmVkjucyw7Pafj4uAwYh4PiLeBO6j2IzTzKznakMFqxy91u2c90zgpdL5EHBx+QJJy4HlAFOY8tHTOKN7tTOzbL3Gwf0R8e6J3cXT4xup9+cq3nZS7L48AHCGpsXFmvCa5WZ2HHgk1v6kHfepuj9lr3U7eA8B55bOZwG7u1wHM7O6itEmXtuknk3AHEmzgZ8CS4A/6HIdzMzq8iSdBiLiiKSbgIeBKcDKiNjezTqYmY3FaZMGImI9sL7bn2tmNh4vTGVmlqlcRpvkUUszsy6IEEfihEpHFZL+UtJ2SdskrZZ0qqTZkp6QtEvS/ZJObqWuDt5mZiXtmqQjaSbwF8C8iPgQxXO+JcBtwO0RMQc4CCxrpZ4O3mZmSQdmWJ4IvEPSicBpwB7gMmBten8VcE0rdXXO28yspInAPF3S5tL5QJpkCEBE/FTSV4AXgV8A3wG2AK9ExJF02RDFzPOmOXibmSVNjvPeHxHzGr0p6UyKtZtmA68A3wCurvuxLXDwNjMraeM47yuAH0fEzwAkfQv4bWCqpBNT77vlWebOeZuZJRFwZOSESkcFLwLzJZ0mScDlwHPAo8C16ZqlwIOt1NXB28yspF0PLCPiCYoHkz8AnqWItwPA54DPShoEzgLubqWeTpuYmSXtXtskIm4BbhlV/DzF3gYT4uBtZlYSnh5vZpYfL0xlZpaZCC9MZWaWITFcbSRJzzl4m5mVOOdtZpYZr+dtZpajKPLeOZhQckfSC5KelfR0bYEWSdMkbUhr1W5I8/tR4Q5Jg5K2SrqwHQ0wM2unEVTp6LV2ZOY/HhFzSwu03AxsTGvVbkznUCzIMicdy4G72vDZZmZtE+mBZZWj1zpRg8UUa9TC29eqXQzcG4XHKRZnmdGBzzcza1lEtaPXJhq8A/iOpC2SlqeycyJiD0B6PTuVzwReKv1uy+vYmpl1SoQqHb020QeWl0TEbklnAxsk/XCMa+u19pi/X+mPwHKAUzltgtUzM6uu6FX3PjBXMaGed0TsTq/7gAcoFlvZW0uHpNd96fIh4NzSr9ddxzYiBiJiXkTMO4lTJlI9M7OmtXkbtI5pOXhLOl3Su2o/A1cC24B1FGvUwtvXql0H3JBGncwHXq2lV8zM+kUuOe+JpE3OAR4o1hjnRODrEfGQpE3AGknLKBYjvy5dvx5YBAwCbwA3TuCzzczaLhAjfTCSpIqWg3dEPA9cUKf8ZYodI0aXB7Ci1c8zM+uGPuhUV+IZlmZmNRk9sHTwNjMry6Tr7eBtZlbinreZWWYCGBlx8DYzy0sA7nmbmeWnH8ZwV5HHgEYzs26Jisc4JL0/LZddO34u6TONls1uloO3mdlR1RalqvJQMyJ2puWy5wIfpZic+ACNl81uioO3mVlZm3reo1wO/CgifkLjZbOb4py3mVlNQHRmtMkSYHX6+W3LZqdVWZvmnreZ2duo4sF0SZtLx/K6d5NOBj4JfKOdtXTP28ysrHpKZH9p+8exXA38ICL2pvO9kmakXnd52eymuOdtZlbW/pz39byVMoHGy2Y3xcHbzKymNkmnylGBpNOABcC3SsW3Agsk7Urv3dpKVZ02MTMraecknYh4AzhrVFndZbOb5eBtZlbmtU3MzPKjTKbHO3ibmdW0NgGnJ8Z9YClppaR9kraVyurOzU+bC98haVDSVkkXln5nabp+l6Sl9T7LzKy3Kj6s7IOVB6uMNrkHWDiqrNHc/KuBOelYDtwFRbAHbgEuBi4Cbml1MRYzs47qzPT4ths3eEfEY8CBUcWN5uYvBu6NwuPA1DQI/SpgQ0QciIiDwAaO/YNgZtZ7IxWPHms1591obv5M4KXSdUOprFH5MdIU0+UAp3Jai9UzM2vBcbwZQ71WxxjlxxZGDAADAGdoWh/848TMjie5jDZpdYbl3pQOYdTc/CHg3NJ1s4DdY5SbmfWXyZLzbqDR3Px1wA1p1Ml84NWUXnkYuFLSmelB5ZWpzMzMWjBu2kTSauBSiuUPhyhGjdwKrJG0DHgRuC5dvh5YBAxS7BpxI0BEHJD0JWBTuu6LETH6IaiZWc/lkjYZN3hHxPUN3jpmbn5EBLCiwX1WAiubqp2ZWTcFnh5vZpalydLzNjM7nkyatImZ2XHFwdvMLEMO3mZmeVE4bWJmliePNjEzy4973mZmOXLwNjPLjHPeZmaZyiR4t7owlZnZpKSRakele0lTJa2V9ENJOyT9VqNtJJvl4G1m1jn/AjwUEb8JXADsoPE2kk1x8DYzK2vTet6SzgA+BtwNEBFvRsQrNN5GsikO3mZmNfHWRJ3xjgreC/wM+A9JT0n6qqTTGbWNJHD2WDdpxMHbzKyses97uqTNpWP5qDudCFwI3BURHwH+jxZTJPV4tImZWVn10Sb7I2LeGO8PAUMR8UQ6X0sRvPdKmpE2by9vI9kU97zNzBLRvtEmEfG/wEuS3p+KLgeeo/E2kk0ZN3hLWilpn6RtpbK/k/RTSU+nY1Hpvc9LGpS0U9JVpfKFqWxQUtv+6WBm1jbtzXkD/DnwNUlbgbnAP1BsI7lA0i5gQTpvWpW0yT3AvwL3jiq/PSK+Ui6Q9AFgCfBB4FeBRyT9Rnr7zlTRIWCTpHUR8VwrlTYz65g2TtKJiKeBeqmVY7aRbFaVPSwfk3RexfstBu6LiEPAjyUNAhel9wYj4nkASfelax28zay/HAczLG+StDWlVWozhGYCL5WuGUpljcqPIWl57entYQ5NoHpmZs1rc9qkY1oN3ncBv06Rw9kD/GMqr7cQboxRfmxhxEBEzIuIeSdxSovVMzNrUZsm6XRaS0MFI2Jv7WdJ/w58O50OAeeWLp0F7E4/Nyo3M+sPUX3dkl5rqeedxibWfAqojURZByyRdIqk2cAc4ElgEzBH0mxJJ1M81FzXerXNzDpksvS8Ja0GLqWYTTQE3AJcKmkuRRNeAP4UICK2S1pD8SDyCLAiIobTfW4CHgamACsjYnvbW2NmNkH9kM+uospok+vrFN89xvVfBr5cp3w9sL6p2pmZddtkCd5mZseNPkmJVOHgbWaWiEmUNjEzO544eJuZ5cjB28wsQw7eZmaZ6ZOp71U4eJuZlTl4m5nlJ5fp8Q7eZmYlTpuYmeXGk3TMzDLl4G1mlhfPsDQzy5RG8ojeDt5mZjXOeZuZ5clpEzOzHDl4m5nlp509b0kvAK8Bw8CRiJgnaRpwP3AexU5kvx8RB5u997h7WEo6V9KjknZI2i7p06l8mqQNknal1zNTuSTdIWlQ0lZJF5butTRdv0vS0mYra2bWce3fw/LjETE3Iual85uBjRExB9iYzptWZQPiI8BfRcT5wHxghaQPjFGBqyk2Hp4DLAfugiLYU+x/eTFwEXBLLeCbmfWFtHt8lWMCFgOr0s+rgGtaucm4wTsi9kTED9LPrwE7gJljVGAxcG8UHgempt3mrwI2RMSB9E+EDcDCViptZtYJtXHeVQ6KTdk3l47ldW4ZwHckbSm9f05E7IEivgJnt1LXpnLeks4DPgI8MboCkmoVmAm8VPq1oVTWqNzMrH9E5ZzI/lIqpJFLImJ3io8bJP1wYpV7S5W0CQCS3gl8E/hMRPx8rEvrlMUY5aM/Z3ntL9lhDlWtnplZWzTR8x5XROxOr/uAByhSxntTNoL0uq+VelYK3pJOogjcX4uIb6XiRhUYAs4t/fosYPcY5W8TEQMRMS8i5p3EKc20xcxsYqo+rKwQvCWdLuldtZ+BK4FtwDqgNmBjKfBgK1WtMtpEwN3Ajoj4p9JbjSqwDrghjTqZD7ya0isPA1dKOjM9qLwylZmZ9Y02PrA8B/i+pGeAJ4H/ioiHgFuBBZJ2AQvSedOq5LwvAf4YeFbS06nsb9IHrpG0DHgRuC69tx5YBAwCbwA3AkTEAUlfAjal674YEQdaqbSZWae0azOGiHgeuKBO+cvA5RO9/7jBOyK+T/18NfUqEBEBrGhwr5XAymYqaGbWNUEzDyx7yjMszcxKvLaJmVmOHLzNzPLizRja5DUOvv5IrN3Z63q02XRgf68r0WaTrU2TrT1wfLTp1yZ8xwhvxtAmOyvMYMqKpM1uU3+bbO0Bt6kpecTuvg/eZmZd5bSJmVluAnDapC0Gel2BDnCb+t9kaw+4TdXlEbv7O3hHxKT7wrlN/W+ytQfcpmY4bWJmliGPNjEzy03zW5z1TOX1vLtN0kJJO9NemC3t8dYtklZK2idpW6ks2z0+J+O+pZJOlfSkpGdSm76QymdLeiLV735JJ6fyU9L5YHr/vNK9Pp/Kd0q6qjctOlqXKZKekvTtdJ57e16Q9KykpyVtTmVd+94Vk3Si0tFrfRm8JU0B7qTYD/MDwPUq9s3sV/dw7JZuOe/xORn3LT0EXBYRFwBzgYVpyeLbgNtTmw4Cy9L1y4CDEfE+4PZ0Hem/wxLggxT/m/9b+r72yqcptiasyb09UH3D3s5870YqHj3Wl8Gb4j/4YEQ8HxFvAvdR7I3ZlyLiMWD08rbZ7vE5GfctTXV7PZ2elI4ALgPWpvLRbaq1dS1wuSSl8vsi4lBE/Jhi6eOLutCEY0iaBfwu8NV0LjJuzxi6+r1zz3tiJsN+l402Gc1qj0+NsW8pmbUppRieptj1aQPwI+CViDhSp35H657efxU4i/5q0z8Df81b/cCzyLs90NyGve3/3rVxJ51O69cHlpX2u8zUhPb47CaN2re06KjVv7ROWd+1KSKGgbmSplLsJ3h+vcvSa1+3SdIngH0RsUXSpbXiOpdm0Z6SZjbs7UCb8lnbpF973pX2u+xzHdnjs1vUxX1Luy0iXgG+R5HPnyqp1okp1+9o3dP7v0KRGuuXNl0CfFLSCxRpxcsoeuK5tgdoesPeznzvIqodPdavwXsTMCc9OT+Z4oHKuh7XqVnZ7vGZcqGTat9SSe9OPW4kvQO4giKX/yhwbbpsdJtqbb0W+G7aJWodsCSN3phN8bDsye604i0R8fmImBUR51H8/+O7EfGHZNoeaGnD3vZ/76Kte1h2VF+mTSLiiKSbKP6DTwFWRsT2HlerIUmrgUuB6ZKGKJ5057zH52Tct3QGsCqNpDgBWBMR35b0HHCfpL8HnqL4o0V6/U9JgxQ91CUAEbFd0hrgOYpROStSOqZffI5823MO8EBKz50IfD0iHpK0iW5+7/qgV12FIpOKmpl12hnvnBkXf/jPKl37yP/87ZZeLrPblz1vM7Ne0Ugf5EQq6Nect5lZ9wVtnaSjijNgW+HgbWaWiGoTdJqYpFN1BmzTHLzNzMraNFSwyRmwTXPO28ysrHqverrS4lnJwKg1xmszYN+VzseaAds0B28zs5pazrua/Y1Gm7QwA7ZpDt5mZiVtGm1SmwG7CDgVOIPSDNjU+57QbFbnvM3MjqqY7x4ntdLCDNimOXibmdUEnV7b5HPAZ9NM17N4awZs05w2MTMra/McnYj4HsVCaETE87RpvXQHbzOzkn7YaKEKB28zszIHbzOzzETAcB5rmzh4m5mVuedtZpYhB28zs8wEkMkelg7eZmZHBYRz3mZmeQn8wNLMLEvOeZuZZcjB28wsNxNat6SrHLzNzGoCyGQDYgdvM7My97zNzHLj6fFmZvkJCI/zNjPLkGdYmpllyDlvM7PMRHi0iZlZltzzNjPLTRDDw72uRCUO3mZmNV4S1swsU5kMFTyh1xUwM+sXAcRIVDrGI+lUSU9KekbSdklfSOWzJT0haZek+yWd3EpdHbzNzGoibcZQ5RjfIeCyiLgAmAsslDQfuA24PSLmAAeBZa1U1cHbzKwkhocrHePep/B6Oj0pHQFcBqxN5auAa1qpp3PeZmbJaxx8+JFYO73i5adK2lw6H4iIgfIFkqYAW4D3AXcCPwJeiYgj6ZIhYGYrdXXwNjNLImJhm+83DMyVNBV4ADi/3mWt3NtpEzOzDouIV4DvAfOBqZJqHedZwO5W7ungbWbWAZLenXrcSHoHcAWwA3gUuDZdthR4sKX7RyZTQc3MciLpwxQPJKdQdJTXRMQXJb0XuA+YBjwF/FFEHGr6/g7eZmb5cdrEzCxDDt5mZhly8DYzy5CDt5lZhhy8zcwy5OBtZpYhB28zswz9P8Ti+Vc+sZ5sAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(grid2, vmin=30, vmax=100)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-32.350648676314037"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid2.mean() - raw.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:netcdf]",
   "language": "python",
   "name": "conda-env-netcdf-py"
  },
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}