DIID ProblemSet 1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem Set 1, due March 1 by midnight\n",
    "\n",
    "The goal of this problem set is to replicate and extend the results of  Jean et al.'s 2016 paper, \"Combining satellite imagery and machine learning to predict poverty.\" This problem set will be challenging and time-consuming, so I suggest you start immediately. Your first step should be to carefully read <a href=\"https://pdfs.semanticscholar.org/1b3a/c4b4187a3dbc9373869e7774b1dc63f748d2.pdf\">the original paper</a>  as well as the <a href=\"http://science.sciencemag.org/content/sci/suppl/2016/08/19/353.6301.790.DC1/Jean.SM.pdf\">supplementary materials</a>.\n",
    "\n",
    "For this assignment, we will focus on the country of Rwanda. You will need to download three distinct datasets, including DHS data, satellite data from the Google Maps API, as well as nighttime luminosity data. The DHS data requires registration (which can take several days to be approved), and the Google Maps API is rate-limited, so it will necessarily take you several days to download the requisite data, so make sure to **get started on those steps asap**. The deep learning section may also take several hours to compute (or days, if you have a slow computer), so don't save it until the last minute.\n",
    "\n",
    "## Overview of the problem set\n",
    "\n",
    "These are the key steps in the problem set:\n",
    "\n",
    "1. Download satellite night lights images from NOAA\n",
    "2. Download DHS data for Rwanda\n",
    "3. Test whether night lights data can predict wealth, as observed in DHS\n",
    "4. Download daytime satellite imagery from Google Maps\n",
    "5. Test whether basic features of daytime imagery can predict wealth\n",
    "6. Extract features from daytime imagery using deep learning libraries\n",
    "7. Replicate final model and results of Jean et al (2016)\n",
    "8. Construct maps showing the predicted distribution of wealth in Rwanda\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Download nightlights for Rwanda\n",
    "\n",
    "- **INPUT**:\n",
    " - None\n",
    "- **OUTPUT**: \n",
    " - `F182010.v4d_web.stable_lights.avg_vis.tif`: Single image file giving nightlights intensity around the world\n",
    "\n",
    "Go to the [DMSP-OLS website](https://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html) and download the satellite nighttime luminosity data (roughly 400MB). We will use the one from 2010. The archive they provide constains several files. Feel free to explore these files. We will only be using the file F182010.v4d_web.stable_lights.avg_vis.tif.\n",
    "\n",
    "A code snippet to get you started is below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'F182010.v4.tar'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import wget\n",
    "\n",
    "night_image_url = 'https://ngdc.noaa.gov/eog/data/web_data/v4composites/F182010.v4.tar'\n",
    "wget.download(night_image_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Download Rwandan DHS and construct cluster-level aggregates\n",
    "\n",
    "- **INPUT**: \n",
    "  - `rwanda_clusters_location.csv`: Coordinates of the centroid of each cluster\n",
    "- **OUTPUT**: \n",
    "  - `rwanda_cluster_avg_asset_2010.csv`: Comma-delimited file indicated average wealth of each cluster \n",
    "\n",
    "[Demographic and Health Surveys (DHS)](http://dhsprogram.com/What-We-Do/Survey-Types/DHS.cfm) are nationally-representative household surveys that provide data for a wide range of monitoring and impact evaluation indicators in the areas of population, health, and nutrition. For this assignment, you will need to download the [2010 Rwandan DHS data](http://dhsprogram.com/what-we-do/survey/survey-display-364.cfm). **This requires registration, so start early!** Do not forget to request for the GPS dataset. Make sure you understand the structure of the data before starting.\n",
    "\n",
    "Your immediate goal is to take the raw survey data, covering 12,540 households, and compute the average household wealth for each survey cluster (think of a cluster as a village). Refer to the file `Recode6_DHS_22March2013_DHSG4.pdf` for information on these data.\n",
    "\n",
    "Save your output as `rwanda_cluster_avg_asset_2010.csv` and check that it matches the file that we have provided. You will use this file as input to the next step in the assignment.\n",
    "\n",
    "Hints:\n",
    "- `Household Recode` contains all the attributes of each household. It provides datasets with different formats. Feel free to explore the data. You can use `RWHR61FL.DAT` file in Flat ASCII data (.dat) format.\n",
    "- `RWHR61FL.DCF` describes the attributes and the location of each attribute.\n",
    "- Geographic Datasets: `rwge61fl.zip` contains the location of each cluster in Rwanda. It is in the format of shapefile, which needs QGIS or other GIS softwares to open. For those who are not familiar with GIS tools or who want a shortcut, you can also sue the file `rwanda_clusters_location.csv` provided with the problem set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For reference, the cluster locations, overlaid on the nightlights data, are shown in the figure below.\n",
    "<img src=\"figure/map1.png\" alt=\"Map\" style=\"width: 600px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From RWHR61FL.DAT file,\n",
    "* Cluster number: startPos=16, len=8\n",
    "* Wealth Index: startPos=230, len=1\n",
    "\n",
    "\n",
    "\n",
    "Label=Wealth index\n",
    "Name=HV270_VS1\n",
    "Value=1;Poorest\n",
    "Value=2;Poorer\n",
    "Value=3;Middle\n",
    "Value=4;Richer\n",
    "Value=5;Richest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cluster_wealth_idx = {}\n",
    "# Read .DAT file\n",
    "surveyfile = open('RW_2010_DHS_02132017_1736_101884/rwhr61fl/RWHR61FL.DAT','r')\n",
    "for line in surveyfile:\n",
    "    line = line\n",
    "    # Just store the required fields here: cluster number & wealth index\n",
    "    clusterNum = line[16:23].strip()\n",
    "    wealthIdx = line[231:239].strip()\n",
    "    \n",
    "    \n",
    "    try:\n",
    "        cluster_wealth_idx[clusterNum].append(float(wealthIdx.strip()))\n",
    "    except:\n",
    "        cluster_wealth_idx[clusterNum] = [float(wealthIdx)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from statistics import median\n",
    "from math import pow\n",
    "\n",
    "def calc_avg_wealth_idx(cluster_dict):\n",
    "    avg_wealthidx_dict = {}\n",
    "    for cluster in cluster_dict:\n",
    "#         mean = sum(cluster_dict[cluster]) / float(len(cluster_dict[cluster]))\n",
    "        m = median(cluster_dict[cluster])\n",
    "        avg_wealthidx_dict[cluster] = m * pow(10,-5)\n",
    "    \n",
    "    return avg_wealthidx_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "median_idx = calc_avg_wealth_idx(cluster_wealth_idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Join median wealth and geo coords\n",
    "import pandas as pd\n",
    "\n",
    "medianWealth_df = pd.DataFrame()\n",
    "medianWealth_df['DHSCLUST'] = median_idx.keys()\n",
    "medianWealth_df['DHSCLUST'] = medianWealth_df['DHSCLUST'].astype('float')\n",
    "medianWealth_df['median_wealth'] = median_idx.values()\n",
    "geocoords_data = pd.read_csv('provided/rwanda_clusters_location.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Merge datasets\n",
    "rwanda_clusters_df = medianWealth_df.merge(geocoords_data[['DHSCLUST','X','Y']],on='DHSCLUST')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rwanda_clusters_df.columns = ['cluster','wlthindf','longitude','latitude']\n",
    "rwanda_clusters_df['cluster'] = rwanda_clusters_df['cluster'].astype('int') \n",
    "rwanda_clusters_df.sort_values(by='cluster',inplace=True)\n",
    "rwanda_clusters_df.to_csv('rwanda_cluster_avg_asset_2010.csv',index=False,columns=['cluster','wlthindf','latitude','longitude'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Test whether night lights data can predict wealth, as observed in DHS\n",
    "\n",
    "Now that you have \"ground truth\" measures of average cluster wealth, your goal is to understand whether the nightlights data can be used to predict wealth. First, merge the DHS and nightlights data, and then fit a model of wealth on nightlights.\n",
    "\n",
    "## 3.1 Merge nightlights and DHS data at cluster level\n",
    "- **INPUT**: \n",
    " - `F182010.v4d_web.stable_lights.avg_vis.tif`: Nightlights data, from Step 1\n",
    " - `rwanda_cluster_avg_asset_2010.csv`: DHS cluster averages, from Step 2\n",
    "- **OUTPUT**: Merged datasets\n",
    " - `DHS_nightlights.csv`: Merged dataset with 492 rows, and 6 columns (one indicates average cluster wealth, 5 nightlights features)\n",
    " - Scatterplot of nightlights vs. DHS wealth\n",
    "\n",
    "Perform a \"spatial join\" to compute the average nighttime luminosity for each of the DHS clusters. To do this, you should take the average of the luminosity values for the nightlights locations surrounding the cluster centroid.\n",
    "\n",
    "Save your output as `DHS_nightlights.csv` and check that it is the same as the file we have provided.\n",
    "\n",
    "Create a scatterplot showing the relationship between average cluster wealth (y-axis) and average nighttime luminosity (x-axis). Your scatterplot should have one dot for each of the 492 DHS clusters. Report the R^2 of the regression line.\n",
    "\n",
    "Hints:\n",
    " - The resolution of each pixel in the nightlight image is about 1km. Use 10 pixels X 10 pixels to average the luminosity of each cluster.\n",
    " - Start by just taking the **Mean** of the luminosity in the 100 pixels and comparing this to cluster average wealth. If you like, you could also compute other luminosity characteristics of each cluster, such as the **Max**, **Min**, **Standard Deviation** of the 100 pixel values, but this step is not required. Note that the file we provide (`DHS_nightlights.csv`) has these added features.\n",
    " - To read the raw raster (nightlights) files, we recommend using the GDAL library. Use `conda install gdal` to install the GDAL library. We have provided some helper code for this below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import os\n",
    "import os.path\n",
    "from osgeo import gdal, ogr, osr\n",
    "from scipy import ndimage\n",
    "from scipy import misc\n",
    "from io import BytesIO\n",
    "gdal.UseExceptions()\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "%matplotlib inline\n",
    "import urllib\n",
    "import numpy as np\n",
    "\n",
    "# Helper function to read a raster file\n",
    "def read_raster(raster_file):\n",
    "    \"\"\"\n",
    "    Function\n",
    "    --------\n",
    "    read_raster\n",
    "\n",
    "    Given a raster file, get the pixel size, pixel location, and pixel value\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    raster_file : string\n",
    "        Path to the raster file\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    x_size : float\n",
    "        Pixel size\n",
    "    top_left_x_coords : numpy.ndarray  shape: (number of columns,)\n",
    "        Longitude of the top-left point in each pixel\n",
    "    top_left_y_coords : numpy.ndarray  shape: (number of rows,)\n",
    "        Latitude of the top-left point in each pixel\n",
    "    centroid_x_coords : numpy.ndarray  shape: (number of columns,)\n",
    "        Longitude of the centroid in each pixel\n",
    "    centroid_y_coords : numpy.ndarray  shape: (number of rows,)\n",
    "        Latitude of the centroid in each pixel\n",
    "    bands_data : numpy.ndarray  shape: (number of rows, number of columns, 1)\n",
    "        Pixel value\n",
    "    \"\"\"\n",
    "    raster_dataset = gdal.Open(raster_file, gdal.GA_ReadOnly)\n",
    "    # get project coordination\n",
    "    proj = raster_dataset.GetProjectionRef()\n",
    "    bands_data = []\n",
    "    # Loop through all raster bands\n",
    "    for b in range(1, raster_dataset.RasterCount + 1):\n",
    "        band = raster_dataset.GetRasterBand(b)\n",
    "        bands_data.append(band.ReadAsArray())\n",
    "        no_data_value = band.GetNoDataValue()\n",
    "    bands_data = np.dstack(bands_data)\n",
    "    rows, cols, n_bands = bands_data.shape\n",
    "\n",
    "    # Get the metadata of the raster\n",
    "    geo_transform = raster_dataset.GetGeoTransform()\n",
    "    (upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = geo_transform\n",
    "    \n",
    "    # Get location of each pixel\n",
    "    x_size = 1.0 / int(round(1 / float(x_size)))\n",
    "    y_size = - x_size\n",
    "    y_index = np.arange(bands_data.shape[0])\n",
    "    x_index = np.arange(bands_data.shape[1])\n",
    "    top_left_x_coords = upper_left_x + x_index * x_size\n",
    "    top_left_y_coords = upper_left_y + y_index * y_size\n",
    "    # Add half of the cell size to get the centroid of the cell\n",
    "    centroid_x_coords = top_left_x_coords + (x_size / 2)\n",
    "    centroid_y_coords = top_left_y_coords + (y_size / 2)\n",
    "\n",
    "    return (x_size, top_left_x_coords, top_left_y_coords, centroid_x_coords, centroid_y_coords, bands_data)\n",
    "\n",
    "\n",
    "# Helper function to get the pixel index of the point\n",
    "def get_cell_idx(lon, lat, top_left_x_coords, top_left_y_coords):\n",
    "    \"\"\"\n",
    "    Function\n",
    "    --------\n",
    "    get_cell_idx\n",
    "\n",
    "    Given a point location and all the pixel locations of the raster file,\n",
    "    get the column and row index of the point in the raster\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    lon : float\n",
    "        Longitude of the point\n",
    "    lat : float\n",
    "        Latitude of the point\n",
    "    top_left_x_coords : numpy.ndarray  shape: (number of columns,)\n",
    "        Longitude of the top-left point in each pixel\n",
    "    top_left_y_coords : numpy.ndarray  shape: (number of rows,)\n",
    "        Latitude of the top-left point in each pixel\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    lon_idx : int\n",
    "        Column index\n",
    "    lat_idx : int\n",
    "        Row index\n",
    "    \"\"\"\n",
    "    lon_idx = np.where(top_left_x_coords < lon)[0][-1]\n",
    "    lat_idx = np.where(top_left_y_coords > lat)[0][-1]\n",
    "    return lon_idx,lat_idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# this illustrates how you can read the nightlight image\n",
    "raster_file = 'data/nighttime_image/F182010.v4d_web.stable_lights.avg_vis.tif'\n",
    "x_size, top_left_x_coords, top_left_y_coords, centroid_x_coords, centroid_y_coords, bands_data = read_raster(raster_file)\n",
    "\n",
    "# save the result in compressed format - see https://docs.scipy.org/doc/numpy/reference/generated/numpy.savez.html\n",
    "np.savez('nightlight.npz', top_left_x_coords=top_left_x_coords, top_left_y_coords=top_left_y_coords, bands_data=bands_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of each pixel:  0.008333333333333333\n"
     ]
    }
   ],
   "source": [
    "print('Size of each pixel: ',x_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['top_left_x_coords', 'top_left_y_coords', 'bands_data']\n"
     ]
    }
   ],
   "source": [
    "npzfile = np.load('nightlight.npz')\n",
    "print(npzfile.files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Look through 10 pixels x 10 pixels around the centroid cluster\n",
    "rwanda_avgwealth = pd.read_csv('rwanda_cluster_avg_asset_2010.csv')\n",
    "\n",
    "# Loop through each cluster and find the location of cluster centroid\n",
    "lonidx_arr = []\n",
    "latidx_arr = []\n",
    "for index,row in rwanda_avgwealth.iterrows():\n",
    "    lon_ix,lat_ix = get_cell_idx(row['longitude'],row['latitude'],\n",
    "                                 npzfile['top_left_x_coords'],npzfile['top_left_y_coords'])\n",
    "    lonidx_arr.append(lon_ix)\n",
    "    latidx_arr.append(lat_ix)\n",
    "\n",
    "rwanda_avgwealth['lat_idx'] = latidx_arr\n",
    "rwanda_avgwealth['lon_idx'] = lonidx_arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def area_luminosity(X,Y,aggfunc=np.mean):\n",
    "    return aggfunc(bands_data[X-5:X+5,Y-5:Y+5,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Proxima/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:2: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "rwanda_avgwealth['lum_mean'] = rwanda_avgwealth.apply(lambda item: area_luminosity(item.lat_idx,item.lon_idx),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rwanda_avgwealth[['cluster','lum_mean','wlthindf']].to_csv('DHS_nightlights.csv',index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import linear_model\n",
    "from sklearn.cross_validation import train_test_split\n",
    "\n",
    "reg = linear_model.LinearRegression()\n",
    "# Split dataset into test and train\n",
    "X_train, X_test, y_train, y_test = train_test_split(rwanda_avgwealth['lum_mean'], rwanda_avgwealth['wlthindf']\n",
    "                                                    , test_size=0.2, random_state=0)\n",
    "reg.fit(X_train.to_frame(), y_train.to_frame())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11f77f2e8>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X903HWd7/HnOzRpU9oUKqFoCwnIYiuKpIdUuXU1qS0g\nHmU97gGLPxbIRaGyq6CrCGelwOFecK/2qCsuQrmru9uWdXddfqwK7ZKg4MVk07Qg/UFRUwShGRUj\nhdCkzef+8Z1JJpOfM9/vzHe+n3k9zslJ8s1k5vPNJK/55P358TXnHCIi4o+quBsgIiLRUrCLiHhG\nwS4i4hkFu4iIZxTsIiKeUbCLiHgmsmA3syoz225m90V1nyIikr8oe+yfBnZFeH8iIlKASILdzJYA\n5wN3RXF/IiJSuKh67BuAvwa0jFVEJGahg93M3gcccM7tACz9JiIiMbGwe8WY2f8CPgocBmqB+cC/\nO+c+nnM79eZFRArgnMurwxy6x+6cu845d5Jz7hTgw8DDuaGedVtv32644YbY26Dz07np/Px7K4Tm\nsYuIeGZWlHfmnHsEeCTK+xQRkfyoxx6RlpaWuJtQVD6fn8/nBjq/ShR68HTGD2TmSvVYIiK+MDNc\nqQdPRUSkvCjYRUQ8o2AXEfGMgl1ExDMKdhERzyjYRUQ8o2AXEfGMgl1ExDMKdhERzyjYRUQ8o2AX\nEfGMgl1ExDMKdhERzyjYRUQ8o2AXEfGMgl2kCFKpFF1dXaRSqbibIhVIwS4Ssc2b76GhYSlr1lxB\nQ8NSNm++J+4mSYXRFZREIpRKpWhoWMrAQDtwBvAEtbWt7N+/h/r6+ribJwmkKyiJxKy3t5eamkaC\nUAc4g+rqBnp7e+NrlFQcBbtIhBobGxkc7AWeSB95gqGh/TQ2NsbXKKk4CnaRCNXX17Nx4+3U1rZS\nV7ec2tpWNm68XWUYKSnV2EWKIJVK0dvbS2Njo0JdQimkxq5gFxEpYxo8FfGc5sfLTCjYRRJC8+Nl\nplSKEUkAzY+vXLGUYsxstpn9zMx6zOxJM7sh7H2KyFiaHy/5CB3szrlDQKtzrgk4E3ivma0I3TIR\nGaH58ZKPSGrszrlX0x/OBmYBqrmIREjz4yUfkdTYzawK6AbeCHzTOffFCW6jGrtISJofPz3ffkax\nTXd0zg2nSzFLgLeb2ZujuF8RGau+vp7m5mYvAqsYNHMoEPmsGDP7G+AV59xXc467G24YHVdtaWmh\npaUl0scWkcpV1JlD//3f0NwMt90Gn/98FM2dVEdHBx0dHSOf33jjjaVfeWpmxwFDzrl+M6sFHgRu\ndc79IOd2KsWISNF0dXWxZs0V9Pd3jxyrq1vOtm130NzcXNidPv44nH32yKdP3nMPJ7S2lvQ/prhK\nMa8H2s1sB/Az4MHcUBcRKbZIZw795CdgNhLqP7z1NubWLuRPP3FbIko8WqAkIt7YvPke2trWUV3d\nwNDQfjZuvJ21ay+a+R08/DC85z2jnz/1FKn6+lgXhxXSY59VrMaIiOQq9oyVtWsvYvXqVfk/xoMP\nwnnnjX6+dy+cdhoAvV1d1NQ0MjAwfnFYuQ5ia68YkRKq5E28SjVjJa+ZQw88EJRcMqH+i1+AcyOh\nDgldHOacK8lb8FAilWvTpi2utnahW7BguautXeg2bdoSd5Py0tfX5zo7O11fX19B31tbu9DBThck\n505XW7uwoPuKxEc/6tINca662rne3ilvnnnu6uqaSv7cpbMzv7zN9xsKfVOwSyUru2DLU9gXpc7O\nTjd/ftNIloJzdXVNrrOzs0gtnsQHP+jGNOJnP5vxt4Z5YQujkGBXKUYkpJmUV5K8iVcqlaKtbR0D\nA+3093czMNBOW9u6vMpJ27fv4OWX9xBbOeOcc4KSy/e/n2lQEO0rZr6tVZIWhynYJVZJrznPtG6c\nyDptWtgXpVQqxdVXXwusB1qBtwHvYMOGW4sfkitXBoG+dWvw+VNPBYHe1FTcx41bvl38Qt9QKUZy\n+FBzzqe8EledNmwJIWwZqbOz0y1YsDz9vX0OOt28eW8pbhkmu9wCzj39dPEeq8hQjV2SIuk1Z+dy\nA2tmdeNS12mjevEM86JU0uc6N9B/9avoH6PEFOySGIWEYrkp9xenqNsX5kWp6P+t5AZ6T0+09x+j\nQoJdC5QkFmNrzsFqvqTUnDMye6S3tbWOWelYLoNrmdp4VAtr6uvrCz63iRYORbJYyXIWZO7ZA296\nU2H35ZN8XwkKfUM9dskR59zgKMU1DW465fwfRagS0fDw+B76L39ZvMbGjAJ67NorxlNJudhAUtqZ\nVKH3TimCgrfXdQ6qcibyPfccLF5czObGTnvFCDD6x1xTE5Q7yuGPeTJh/r2X6RW8d0oR5V0iOnIE\nZuVEVV8fqfR9NdbUlMV5lRP12D1T1IsNiERgxr+jQ0NQUzP2m3//ezj22ER1XsKK7dJ4Uj6SvMKx\n3CV9MdVMFfs8p70w92uvBYOi2aH+8stBKebYYyNZCes7BbtnkrzCsZxVyrU0S3GeqVSKU089he7u\nR9m27Q7279/D6tWr6P7xj4NAr60dvfGrrwaBPm/eyCF1XmYg39HWQt/QrJiS8WW2Sbko59klUSrF\neU40G+Z7d909fpbL4GCs7SwnaIGSZJTrFLwkinIxVTk/L8VeNJYbyK+jY1ygV7F9RiE9WeelnH++\nhVKwixRBVD3Ect8bp9g94cwLxxKeHRfoxpG8X0xyQ7zcf76FUrCLFEnY8lZSygfFLOP97rHHxgX6\nnNkLIvm5JOXnW4hCgl3z2EVmIOx88KiX9xdLUea9d3bC29/OwqxDC+qaGBraz90bbwcIvS1DUn6+\npaJ57CJTiGplbEWuL2hvh1WrxhxK9fXR09MDQFNT08i5h/05+/zzLWQeu0oxIpOIumZbMbOV7r9/\nXMnFueLXwKf7+SZ1YBXV2JMhqb9glaRYNVuvn/tNmyYMdOdKVwOf7Oeb5IHVQoJdC5RKrFIWuiRd\nsRbBlPK6mSVbKfvtbwcLiy6+ePRYJtrTSrWoqL6+nsbGRnp7e0fOuyJXqub7SlDoG+qxez1y75sk\nP1d9fX3u5ptvcXPmHFPcHuqXvzxpD32iNpXi5zlRzzzpF3VBpZjylvRfsEqTxJp4ps1wqoNjHWyJ\nPkSvu27GgT5R24r185zsxWPXrl2JfZF2TsFe9pLcC6xUSaqJT/T7BQsd9EXTgbjyyoICPbeNxfp5\nTtVxSuKLdEYhwa557CVU7pdSk/GStF/8RHO5oQHYGm4juIsvhs2bxx5zhU1dLubPc6rLLTY3N5fd\nvvTFFHoeu5ktAb4LLAKGgTudc1+f4HYu7GP5QlcNkmKYaC43nM2cOTXcffff579f+bnnwkMPjT1W\n5n/D5XjFqLAKmcceRbCfAJzgnNthZvOAbuAC59yenNsp2EWKLDfYrrvus3zyk5fn14E46yzo7h57\nLEF/u7t376azs5MVK1awbNmyuJsTWizBPkEj/gP4hnPuv3KOK9hFSqDg/whPPhlypx8m7G/Wxysr\nxR7sZtYIdABvcc4dzPmagl2kHB19dHBBi2wJ/Fv1dVuBWC9mnS7D/Cvw6dxQz1i/fv3Ixy0tLbS0\ntET18CKSL5sgKxIY6Bm9vb3MmtXARIugkhTsHR0ddHR0hLqPSHrsZjYLeAD4oXPua5PcRj12kXKQ\nG+ivex389rfxtCVCd9xxJ1dc8WngcSq9xx7VlgJ3A7smC3URKQNmY0N92bKgh16iUC/mFgepVIqr\nr74WWA+0Am8D3sGGDbcmOtQLFTrYzWwl8BFglZn1mNl2MzsvfNNEJBK5gb5yZRDou3aVrAnF3iNp\ndC+azwN7gLuYN++NLF9+ZqSPkxTaj13EV7kll/e/H+67r+TNKMWgpq8DpxBvKUZEykVuD/2SS4Ie\neoShnk9ZpRQ7O2ZWddfWtlJXt5za2taKXtWtHruIL3J76J/5DGzYEPnD5DtXvJS9aR9Xdcc+j33K\nB1KwixRHbqDfeCN86UtFeahCQ9rHpf6lEus8dhEpsdxA37Ah6KUXUaEXjS7KRbJlUgp2kaTJDfS7\n74ZLLy3JQ0+1g+J0krRTZtJp8FQkKXIHRf/lX4JB0RKFOmiQMilUYxcpd7k99C1b4KJ469M+DlKW\nKw2eivcqKlByA/0//xPOPz+etpRQRT3HM6B57OK1Yq9eLAvOjS+5PPJIcLwCQr0inuMSUI9dEsHn\nlYVAENxVOf2s7m5Yvjye9sTA++e4QOqxi7dKsXoxFkeOBL3z7FDfvTsI+goKdfD4OY6Bgl0SYew0\nO8hnml1ZGhwMAn1W1ozj3t4g0Jcuja1ZcfLuOY6Rgl0SwZtpdq+8EgT67Nmjx154IQj0hob42lUG\nvHmOy4Bq7JIoiZ0x8Yc/wLHHjj320ktwzDHxtKeMJfY5LhJNdxQpNwcOwAknjD128GBwnVGRGdDg\nqUi5ePbZoOSSHeqvvRaUXBTqUmQK9jJUzEuISZHt3RsEena9fGgoCPTsurpIESnYy8zmzfdw0kmn\n0dr6MU466TQt0EiKnp4g0LNntAwPB4E+S3vtSWmpxl5GUqkUixe/kaGhWcDJwC846qhBnnyym2XL\nlsXdPJnIT34C73rX2GPDw+O3A5CS82UQVjX2hOvp6WFo6AjQQXBR3iqOHHkDTU3/Qz33cvOjHwXh\nnR3qzo1uCSCxqvStCdRjLyMPPfQQ5577KeCnwFJAS6vLzve+BxdeOPaYfq/Lyky3JkhKj1499oRr\namqipiYFbAUa0dLqMnLXXUFPPDvUMz10KSsz2ZrA9x69gr2M1NfX8w//cAdz5nwK2IOWVpeBr3wl\nCPTLLx89pkAva9NtTZBKpWhrW8fAQDv9/d0MDLTT1rbOq1loCvYIRDk9ce3ai3j22ae5+ebrtbQ6\nTtdfHwT65z43ekyBngjTbU1QCZuNqcYeUubq6zU1QS8hyquvJ6UG6JV16+Bb3xp7zMPf20ow2d9P\n0rYH1pYCJZa0XxCZwoUXBgOj2Tz7fZVRmQ5ZdXUDQ0P7I+2QRa2QYNfKiRAy/9INDIz/l07BnhCr\nVkF7+9hjCnTvrV17EatXr/L2P2IFewhjB2mCHrsvg5zel4He8hZ46qmxxxToFaW+vt7P320iGjw1\ns41mdsDMnpj+1v7wdf9or6eCHX98MCiaHeoaFBXPRFJjN7N3AgeB7zrnzpjkNt7V2DN86t16O24w\n0WpQT38fxS+x1didc4+aWcVe/qVY/9LF8YLh3biBAl0qkOaxl6mw5ZBC59Z7c91Js/GhrpKLVIiS\nDp6uX79+5OOWlhZaWlpK+fCJkb0yLug5P0FbWyurV6+aUa85zNz6zLhBW1vrmKlgiemtq4cuCdfR\n0UFHR0eo+4hsHnu6FHN/JdbYp5NvSaWrq4s1a66gv7975Fhd3XK2bbuD5ubmaR8rihp54sYNFOji\nqbg3AbP0m2QppKQSphwS1XLp+vp6mpubyz/Uc0suc+ao5CIVL6rpjpsI9po9zcyeNbNLo7jfpCt0\ns6Ew0yi9qZFPJzfQGxuDMB8YiK1JIuUiqlkxF0dxP74JM8Ok0JVxia+RTye35NLcDJ2d8bRFpExp\nr5gIlONmQ4mrkU8nN9Df9z544IF42iJSQnHX2CvSVDX0OFemJqZGPp3cksvHPx6UXBTqIpNSjz0E\n3y7BVVZye+jXXBNc9CJGeh4lDuqxl9joDJTXA13A6yecgeJN77kUcnvoN90U9NBjDnWv988R76jH\nHkIqlWLx4jcyNDQLOBn4FdXVQzz//C8V4vnK7aF/4xtw1VXxtCWHt/vnSCKoxx4DsyqgA+gGOjA7\nKt4GJU1uD/273w166GUS6lAZl1ITvyjYQ+jt7aW29o1k/8HPmXOK/uBnIjfQ7703CPSPfSy+Nk2i\nYtYGiDcU7CHoD74AuYF+331BoH/gA/G1aRq+7rsv/lKNPaTMhltVVUsYHn5uwg23Kn42hXNQldOH\nePhhaG2Npz0FqvjnUWKhGntMnBsGDqXfj1XRsymGh4PeeXaod3UFQZ+wUAfNbpLkUI89hOlmS1Ts\nbIqhIaipGXts1y5Ytiye9ogkmHrsJTbdbImKm03x6qtBDz071Ht7gx66Ql2kZBTsIUw3eFoxg6v9\n/UGgH3306LEXXwwCvaFir5goEhsFewjTzZbwfjZFX18Q6MccM3rspZeCQF+0KL52yRiFXiZRkks1\n9ghMN1vCu9kU+/cH+59ne+UVmDs3lubI5MJcJlHKQyE1dgW7zNyePeNr5YcOjR8olbJQsYP3ntHg\nqRRHd3dQcskO9cOHg5KLQr1sVdzgvYxQsMcgMTXPjo4g0M86a/TY8HAQ6EdpT5xyVzGD9zKOgr3E\nErFg6YEHgkDPXkSUuUB07i6MUra8H7yXSanGXkJlX/P8p38avwlXhT9nPvBu8L7CFFJjj+Ri1jIz\nYS5uXVTf/Ob4bXIV6N6or69XoFcYlWJKqOxqnjfdFJRWskM9U3IRkcRSsEdgpoOhZVPz/MxngkC/\n4YbRYwp0EW+oxh7S6La9JzI8/OsZLQCJreb5kY/Apk1jj3n4nIj4RAuUSiyVSrFkyZ8wOPhjMoOh\nNTXv4rnn9pVXTfOcc2Dr1rHHPHsuRHylBUol1tPTw+BgPdkLQAYHj6OnpyfOZo0644yg5JId6iq5\niHhPs2JC+w3BxayPBl4BXoi1NUCwAVdf39hjCnORiqFSTAipVIoTTmhgeNgBxwN9VFXBiy8+G08p\nZqLFQ579zEUqTWylGDM7z8z2mNnTZvaFKO4zKYLcrAWOA2rJ8+cfjdwLRGcaplAXqUihg93MqoC/\nA84FTgfWmtnSsPebBO3t7ekg7wC6gQ6cM9rb24v2mGOmVk4S6Km+vmTsRSMiRRFFj30FsM85t985\nNwRsAS6I4H7L3oEDB4DXkz14Cq9n586dRQnWzD4zzStWUH/88WO/mO6hJ2IvGhEpqtA1djP7EHCu\nc+4T6c8/Cqxwzv1Vzu28q7E/9thjvPOda4DHyUx3hHcAw8yZ8yeY/SayCxukUqnxYQ6k+vpG6vll\nvxeNaN8WyVvZ7xWzfv36kY9bWlpoaWkp5cNH7pVXXgEWAK3AEuA5oA74Kq+9djHwBG1traxevSrc\nH7EZud9tOOrqlrMta5+Zst2LRgBdzUhmpqOjg46OjnB34pwL9UbQRf1R1ufXAl+Y4HbONw8++KCD\nuQ7aHXSm38910OcytZH58890nZ2dhT3A6BDoyBvsTH+409XWLnR9fX0jN+/r63O1tQunvI3EQ8+N\nFCqdnXnlchQ19i7gVDNrMLMa4MPAfRHcb9k78cQTgUPA+cBa4IMEwxaZuexPMDjYy7x58/KruU8y\nKLp505Zp95m57rrPxr8XjYyjqxlJSeX7SjDRG3AesBfYB1w7yW2K/LpWenfeeaeD2Q4WOHiDgzkO\nFjmodfBmB7Xusss+4WprF7oFC5a72tqFbtOmLZPfYW4Pfe7ccTfp6+tznZ2d43p6mzZtGXmcOXOO\ncTfffMu43vxE3yeloR67FIoCeuyRBPuMHsjDYL/kkkvSpZevO2hzcIuD2a6q6hg3e3ad+9u//crM\n/phzA/2UU/Jqx3ShkR360764SNFknoe6uiY9DzJjhQS7thQIYXh4GHDAFwgGT4OdE9/61ga2bt06\n5WAmMH6Wy8qV8Oijebdjusdpa1vHwEB7+usRDehK3tauvYjVq1dpVowUnYI9hJ///Ofpj8ZOd9y5\ncyeQe2GN4OtDQ/tpXrFizP386m1n8ts7vx38sRfQjskep7GxUTNlyoyuZiSloN0dQ9i3bx+wmLEL\nlBYDjARn9oU1HG/j1YHfj3z/HXwC4zZO2bmX97znkwUvKJrqAh5ld9UmESk6bQIWQlVVFc7NYfwC\npQF27drFsmXLghvmzHD56uxFfPbQi0AKWApEs6BossUvmfnT1dUNDA3t1/xpkQTRhTZKzMyARcAQ\nowuUqoGX+NKXruXGm24ac/uD69fz07PP5oILLuK11x4hmCp5ObBj5DZ1dcvZtu2OkTJKVLVYrXgU\nSaayX3nqpz8ANen3g8DLOIaCC0Vn3HUXm+fOS686vI/Dh48AzQRlmxfJrY1v376Dd7/7vBmtUJxp\nYKu2K1I5VGOPRBVQh+NgEOppH+IrGHO4yzkuu+wKBgba6e/v5vDhR4G5wN8DnyMo35wKnM3hw4Nc\nddXVDAx8g/7+HzEw0E5b27oJFzdpwy8RmYiCPbQqHP04fj5y5Bz+EmMB/858oIpPfeqrvPbaILA7\nfYszgEbgWOAv08duAnoZGvobDh8eBv4PQf1994QrFFOp1Mg0xv7+7ilfAESksqgUE9obgF8AsIr/\nop1PEIT2ccA1wP9jcDAzsNoCrAKeIgj5nxAs1l0MXJw+dgvwrfR9LwAuY3DQjZvFommMIjIZBXto\nL2DsZHRWzAvAI+n3xzN2KuTC9Ps/EIT59Rx1lHHUUTUMDn4Z+N/p77kSOIZM/f766z8/Lqynmrsu\nIpVNpZjQDhP0xJen3x8G7ieYKXOA7PnjQdj/keBSesek31dxyy1fAtYTvCDsI5g+eQj4ATBIa+u7\nxz3qVHPXRaSyqcceiWHglfR7gKMBI5gt86cEpZl9HHWU48iRKoJL6QW97CNHzmb27Grmz1/Kyy9n\n9+4b0vezhGeeeYaVK1eOe1QtUReRiSjYQzsJ+CnQSxDgZwP9BAOjzwKHqan5JV//+tfo7e3l1lvv\nIfdSeosWLeLw4f1kl1VgP8GLxXOsyNmCIJumMYpILpViQnueoMTSnH7/G4IyywGCy7/W8M1vbuCT\nn7yca675DNXVKbLLMzU1v6W1tXWkrFJT82aC6Y+zgfO56qrLR1ewiojMgHrsoc0mCOLMytMa4NfA\naUAbsICBgQEg6F1/5zvfpq2tlaqqJQwPP8fGjXdQX18/pqwyODjIM888w4oVKxTqIpI3bSkQQrCl\nQC3BIOfRBKWT84HbCMI++PzRR7eOqZFreb+IzJS2FIjFYoLZMBlvAK4j6LH3UlW1kJqamjHfobq4\niBSTauwhzJo1i6DGnj2l8TcE0x27gXaGh3/PvHnzYmqhiFQi9dhDqKoK9oiBVoLpib8iGDhtSd/i\nDObMOYWDBw/G00ARqUgK9hBmzZrF4OAfGa2xP0+wNcDotEWzF7QaVERKSsEewllnncWPf7wP+BBB\nj30/mVkyc+eeinPPazWoiJScZsWE8Nhjj/HOd64hd1bMrFlV/PM/301ra6tCXURCKWRWjAZPQ1i5\nciXnnNNCMMXx4vT7Iaqrj+eSS65k27aHC7rfVCpFV1eXtuAVkYIo2EN68MEf8OijW7nyyjVUV1cB\nWxkY+GXB+6Pr4hkiEpZKMRHp6upizZor6O/vHjmWuX5pc3PzjO4jlUrR0LCUgYFoLm4tIsmnUkyM\nxu6PDoXsj565eEb2JmETXT1JRGQqCvaIRLE/ehQvDiIiKsVEKJVK0dPTA0BTU1NB5ZPNm++hrW0d\n1dUNDA3tZ+PG21m79qKomyoiCVFIKSZUsJvZnxNc+mcZ0Oyc2z7Fbb0O9kwg19QEve4NG25l+fIz\nC9roS5uEiUhGHMH+JoLLBt0BfK5Sg338oOeXgfXMn7+Uw4cn73UrwEVkOiUfPHXO7XXO7SO4DlzF\nGjvomSLYtvdxXn55+6TTHjWtUUSKRYOnERg76NkLnMhUM1tSqRRtbesYGGinv7+bgYF2Lr30Cnbv\n3l3ilouIj6YNdjPbamZPZL09mX7//lI0MAmyZ8TMm3cZ8DRTzWyZaFrjoUP1NDW9Qz13EQlt2k3A\nnHNronqw9evXj3zc0tJCS0tLVHcdu+xL223fvoOrr24dM7Mlu4Y+toefuXj17zh06F7a2j7E6tWr\nVHMXqVAdHR10dHSEuo9IpjuaWTvB4Gn3FLfxdvB0ItMNjG7efA+XXnoFhw7VA78Dbgcuynu1qoj4\nLY5ZMX8GfAM4DvgDsMM5995JbltRwT4Tu3fvpqnpHRw6dC/BxTm0hYCIjFXyYM/rgRTsE9KCJBGZ\nioI9oTSfXUQmo2AXEfGMdncUEREFu4iIbxTsIiKeUbCLiHhGwS4i4hkFu4iIZxTsIiKeUbCLiHhG\nwS4i4hkFu4iIZxTsIiKeUbCLiHhGwS4i4hkFu4iIZxTsIiKeUbCLiHhGwS4i4hkFu4iIZxTsIiKe\nUbCLiHhGwS4i4hkFu4iIZxTsIiKeUbCLiHhGwS4i4hkFu4iIZ0IFu5l92cx2m9kOM/s3M6uLqmEi\nIlKYsD32h4DTnXNnAvuAL4ZvUjJ1dHTE3YSi8vn8fD430PlVolDB7pzb5pwbTn/6OLAkfJOSyfdf\nLp/Pz+dzA51fJYqyxn4Z8MMI709ERAowa7obmNlWYFH2IcAB1zvn7k/f5npgyDm3qSitFBGRGTPn\nXLg7MLsEuBxY5Zw7NMXtwj2QiEiFcs5ZPreftsc+FTM7D/hr4F1ThXohDRMRkcKE6rGb2T6gBvhd\n+tDjzrl1UTRMREQKE7oUIyIi5aWoK0/N7M/N7OdmdsTMlud87Ytmti+9wOmcYrajmMzsPDPbY2ZP\nm9kX4m5PWGa20cwOmNkTWceONbOHzGyvmT1oZgvibGMYZrbEzB42s6fM7Ekz+6v0cS/O0cxmm9nP\nzKwnfX43pI97cX4AZlZlZtvN7L705z6dW6+Z7Uw/f53pY3mfX7G3FHgS+CDwSPZBM1sGXAgsA94L\n3G5miavBm1kV8HfAucDpwFozWxpvq0L7vwTnk+1aYJtz7k3AwyR7Idph4Brn3OnA2cCn0s+ZF+eY\nHutqdc41AWcC7zWzFXhyfmmfBnZlfe7TuQ0DLc65JufcivSxvM+vqMHunNvrnNtHMEUy2wXAFufc\nYedcL8Gq1RW5358AK4B9zrn9zrkhYAvBuSWWc+5R4KWcwxcA30l//B3gz0raqAg55150zu1If3wQ\n2E2wsM6nc3w1/eFsggkSDk/Oz8yWAOcDd2Ud9uLc0ozxuZz3+cW1Cdhi4NdZnz+fPpY0uefxHMk8\nj+kc75w7AEEwAsfH3J5ImFkjQa/2cWCRL+eYLlX0AC8CW51zXfhzfhsIZuJlDw76cm4QnNdWM+sy\ns/+ZPpbidaWPAAAB10lEQVT3+YWa7ggzW8Ak3kn8iLuZzQP+Ffi0c+7gBOssEnuO6W0+mtKb8n3f\nzE5n/Pkk7vzM7H3AAefcDjNrmeKmiTu3LCudcy+YWT3wkJntpYDnLnSwO+fWFPBtzwMnZn2+JH0s\naZ4HTsr6PKnnMZ0DZrbIOXfAzE4A+uJuUBhmNosg1P/ROXdv+rBX5wjgnPujmXUA5+HH+a0EPmBm\n5wO1wHwz+0fgRQ/ODQDn3Avp9ykz+w+Ccm/ez10pSzHZdfb7gA+bWY2ZnQycCnSWsC1R6QJONbMG\nM6sBPkxwbklnjH++Lkl//BfAvbnfkDB3A7ucc1/LOubFOZrZcZlZE2ZWC6whGEdI/Pk5565zzp3k\nnDuF4G/tYefcx4D7Sfi5AZjZ3PR/kpjZ0cA5BBNQ8n/unHNFeyMo8v8aGABeAH6Y9bUvAs8Q/NKd\nU8x2FPkczwP2EgwAXxt3eyI4n03Ab4BDwLPApcCxwLb0eT4EHBN3O0Oc30rgCLAD6AG2p5/DhT6c\nI/DW9DntAJ4gKIniy/llnee7gft8Ojfg5KzfyyczeVLI+WmBkoiIZ3RpPBERzyjYRUQ8o2AXEfGM\ngl1ExDMKdhERzyjYRUQ8o2AXEfGMgl1ExDP/HwWZlI587m2JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b843ac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "figure,axes = plt.subplots()\n",
    "plt.scatter(X_train,y_train)\n",
    "plt.plot(X_train,reg.predict(X_train.to_frame()),color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe calculate score for train data first?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R squared for the regression is:  0.727340621199\n"
     ]
    }
   ],
   "source": [
    "print('R squared for the regression is: ',reg.score(X_test.to_frame(), y_test.to_frame()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2. Fit a model of wealth as a function of nightlights\n",
    "- **INPUT**: \n",
    " - `DHS_nightlights.csv`, from Step 3.1\n",
    "- **OUTPUT**: \n",
    " - R^2 of model\n",
    " \n",
    "Above, you fit a regression line to illustrate the relationship between cluster average wealth and corresponding cluster nightlights. Now, use [cross-validation](https://en.wikipedia.org/wiki/Cross-validation_%28statistics%29) to get a better sense of out of sample accuracy.\n",
    "\n",
    "There are two options for this. The basic way, for those new to machine learning, is to randomly divide your dataset into a training and a test dataset. Randomly select 80% of your clusters and fit a model of cluster-average DHS wealth (your response/dependent variable) on nightlights (your predictor/independent variables). You can use a regression or any other model you prefer. Then, use that model to predict the wealth of the remaining 20% of your data, and compare the predicted values to the actual values, and report the R^2 on these 20%.\n",
    "\n",
    "The preferred way is to use 10-fold cross-validation, where you repeat the above procedure 10 times, so that you have 10 different and non-overlapping test sets. Then, you report the cross-validated R^2 of your model (i.e., the average R^2 of your 10 test folds).\n",
    "\n",
    "Hints:\n",
    " - The scikit learn library has built-in functions for [cross-validation](http://scikit-learn.org/stable/modules/cross_validation.html) that make this quite easy.\n",
    " \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nightlights_df = pd.read_csv('DHS_nightlights.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy:  0.751521364962\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cross_validation import KFold\n",
    "from sklearn.cross_validation import train_test_split\n",
    "# from sklearn.linear_model import Ridge\n",
    "\n",
    "# Split dataset into test and train\n",
    "X_train, X_test, y_train, y_test = train_test_split(nightlights_df['lum_mean'], nightlights_df['wlthindf']\n",
    "                                                    , test_size=0.2, random_state=9)\n",
    "\n",
    "# Use cross validation on training data to ensure we aren't overfitting\n",
    "cv = KFold(n=X_train.shape[0],n_folds=10)\n",
    "reg = linear_model.LinearRegression()\n",
    "val_scores = 0\n",
    "\n",
    "for train,test in cv:\n",
    "    reg.fit(X_train.iloc[train].to_frame(),y_train.iloc[train].to_frame())\n",
    "    val_scores += reg.score(X_train.iloc[test].to_frame(),y_train.iloc[test].to_frame())\n",
    "    \n",
    "val_scores /= 10\n",
    "\n",
    "print('Accuracy: ',val_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Score:  0.647575378611\n"
     ]
    }
   ],
   "source": [
    "# Predict on test data\n",
    "test_score = reg.score(X_test.to_frame(),y_test.to_frame())\n",
    "print('Score: ',test_score)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Download daytime satellite imagery \n",
    "- **INPUT**: \n",
    " - Google Maps API key\n",
    " - `Sector_Boundary_2012.shp`: Rwandan shapefile\n",
    "- **OUTPUT**: \n",
    " - Thousands of satellite images (store in directory `google_image/`)\n",
    "\n",
    "We will use the Google Static Maps API to download satellite images. Refer [Google Static Maps introduction](https://developers.google.com/maps/documentation/static-maps/intro) and [Google Static Maps API Usage Limits](https://developers.google.com/maps/documentation/static-maps/usage-limits). You must apply for an API key before downloading. ** Note that it may take you several days to download the required images, so start early!**\n",
    "\n",
    "Download the images from Google at zoom level 16 (pixel resolution is about 2.5m). Set the image size to be 400 pixels X 400 pixels, so that each image you download will cover 1 square kilometer. In this way, each daytime image you download will correspond to a single pixel from the nighttime imagery from Step 1 above.\n",
    "\n",
    "Hints:\n",
    " - You will need to tell Google the locations for which you wish to download images. One way to do this is to use a [shapefiles](https://en.wikipedia.org/wiki/Shapefile) that specifies the borders of Rwanda. We have provided this shapefile (`Sector_Boundary_2012.shp`) as well as a helper function to read in the shapefile.\n",
    " - The function we provide below does not limit the maximum number of images downloaded per day. Note that if you attempt to  download more than the daily limit, Google will return blank images instead of an error.\n",
    " - You can organize the files however you like. However, for later analysis (Steps 6 and beyond), it may help if you organize these daytime images into 64 folders, with one folder indicating the nightlight intensity of the pixel corresponding to the daytime image. In other words, if you download a daytime image for which the corresponding nighttime pixel has value 32, store that daytime image in a folder labeled '32'. This way, all the satellite images within each folder will have the same nightlight intensity. The file name is columnIndex_rowIndex.jpg, in which row index and column index are the index in the nightlight image (See the diagram below)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](figure/data_description.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Helper function to read a shapefile\n",
    "def get_shp_extent(shp_file):\n",
    "    \"\"\"\n",
    "    Function\n",
    "    --------\n",
    "    get_shp_extent\n",
    "\n",
    "    Given a shapefile, get the extent (boundaries)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    shp_file : string\n",
    "        Path to the shapefile\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    extent : tuple\n",
    "        Boundary location of the shapefile (x_min, x_max, y_min, y_max)\n",
    "    \"\"\"\n",
    "    inDriver = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
    "    inDataSource = inDriver.Open(inShapefile, 0)\n",
    "    inLayer = inDataSource.GetLayer()\n",
    "    extent = inLayer.GetExtent()\n",
    "    # x_min_shp, x_max_shp, y_min_shp, y_max_shp = extent\n",
    "    return extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "200\n",
      "300\n",
      "400\n",
      "500\n",
      "600\n",
      "700\n",
      "800\n",
      "900\n",
      "1000\n",
      "1100\n",
      "1200\n",
      "1300\n",
      "1400\n",
      "1500\n",
      "1600\n",
      "1700\n",
      "1800\n",
      "1900\n",
      "2000\n",
      "2100\n",
      "2200\n",
      "2300\n",
      "2400\n",
      "2500\n",
      "2600\n",
      "2700\n",
      "2800\n",
      "2900\n",
      "3000\n",
      "3100\n",
      "3200\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Proxima/anaconda/lib/python3.5/site-packages/PIL/Image.py:870: UserWarning: Palette images with Transparency   expressed in bytes should be converted to RGBA images\n",
      "  'to RGBA images')\n"
     ]
    }
   ],
   "source": [
    "# Helper functions to download images from Google Maps API\n",
    "import urllib.request\n",
    "\n",
    "from retrying import retry\n",
    "\n",
    "# @retry(wait_exponential_multiplier=1000, wait_exponential_max=3600000)\n",
    "def save_img(url, file_path, file_name):\n",
    "    \"\"\"\n",
    "    Function\n",
    "    --------\n",
    "    save_img\n",
    "\n",
    "    Given a url of the map, save the image\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    url : string\n",
    "        URL of the map from Google Map Static API\n",
    "    file_path : string\n",
    "        Folder name of the map\n",
    "    file_name : string\n",
    "        File name\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    None\n",
    "    \"\"\"\n",
    "    with urllib.request.urlopen(url) as url_handle:\n",
    "        a = url_handle.read()\n",
    "#     a = urllib.urlopen(url).read()\n",
    "    b = BytesIO(a)\n",
    "    image = ndimage.imread(b, mode='RGB')\n",
    "    # when no image exists, api will return an image with the same color. \n",
    "    # and in the center of the image, it said'Sorry. We have no imagery here'.\n",
    "    # we should drop these images if large area of the image has the same color.\n",
    "    if np.array_equal(image[:,:10,:],image[:,10:20,:]):\n",
    "        pass\n",
    "    else:\n",
    "        misc.imsave(file_path + file_name, image[50:450, :, :])\n",
    "\n",
    "# Now read in the shapefile for Rwanda and extract the edges of the country\n",
    "inShapefile = \"provided/Sector_Boundary_2012/Sector_Boundary_2012.shp\"\n",
    "x_min_shp, x_max_shp, y_min_shp, y_max_shp = get_shp_extent(inShapefile)\n",
    "\n",
    "left_idx, top_idx = get_cell_idx(x_min_shp, y_max_shp, top_left_x_coords, top_left_y_coords)\n",
    "right_idx, bottom_idx = get_cell_idx(x_max_shp, y_min_shp, top_left_x_coords, top_left_y_coords)\n",
    "\n",
    "\n",
    "key = 'AIzaSyAuMuxiWaz7QJ_2PeZlyVY3iUU3OWGQ9CE'\n",
    "m = 1\n",
    "for i in range(25752, 26831):\n",
    "    for j in range(11884, 11887):\n",
    "        lon = centroid_x_coords[i]\n",
    "        lat = centroid_y_coords[j]\n",
    "        \n",
    "        url = 'https://maps.googleapis.com/maps/api/staticmap?center=' + str(lat) + ',' + \\\n",
    "               str(lon) + '&zoom=16&size=400x500&maptype=satellite&key=' + key\n",
    "#         print('url: ',url)    \n",
    "        lightness = bands_data[j, i, 0]\n",
    "        file_path = 'google_image/' + str(lightness) + '/'\n",
    "        if not os.path.isdir(file_path):\n",
    "            os.makedirs(file_path)\n",
    "        file_name = str(i) + '_' + str(j) +'.jpg'\n",
    "        save_img(url, file_path, file_name)\n",
    "        if m % 100 == 0:\n",
    "            print (m)\n",
    "        m += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 5. Test whether basic features of daytime imagery can predict wealth\n",
    "In step 3, you tested whether nightlight imagery could predict the wealth of Rwandan villages. You will now test whether daytime imagery can predict village wealth. Start by extracting simple metrics from the daytime imagery; in step 6 you will use more sophsticated methods to engineer these features from the images. **You don't need to do this step if you are able to do step 6.**\n",
    "\n",
    "## 5.1. Extract \"basic\" features from daytime imagery\n",
    "- **INPUT**: \n",
    " - `google_image/...`: Raw images, from Step 4\n",
    "- **OUTPUT**: \n",
    " - `google_image_features_basic.csv`: Image features \n",
    "\n",
    "Convert the raw data from the satellite imagery into a set of features that can be used in a machine learning algorithm. A simple way to do this is to take the raw R/G/B values for each pixel and average them for the image. Thus, if an image has 100 pixels, you will have an average R value, an average G value, and an average B value. Create more features by also computing the min, max, median, and standard deviation of R, G, and B for each image. This process will convert each image into a vector of 15 features.\n",
    "\n",
    "Feel free to be creative if you wish to generate additional features from the imagery -- this is similar to the process described in section 2.3 of the paper's supplementary materials. But don't waste too much time, and don't expect these features to be terribly useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Loop through all image files in google_image\n",
    "import os\n",
    "from PIL import Image\n",
    "\n",
    "# index = range(0,50001)\n",
    "color_df = pd.DataFrame(columns=['filename','r_mean','g_mean','b_mean'\n",
    "                                              ,'r_max','g_max','b_max'\n",
    "                                             ,'r_min','g_min','b_min'\n",
    "                                             ,'r_median','g_median','b_median'\n",
    "                                             ,'r_std','g_std','b_std'])\n",
    "for subdir,dirs,files in os.walk('google_image/'):\n",
    "    for file in files:\n",
    "        filepath = subdir+os.sep+file\n",
    "        \n",
    "        if filepath.endswith(\".jpg\"):\n",
    "            # Start processing image here\n",
    "            img = Image.open(filepath)\n",
    "            rgb_arr = np.array(img)\n",
    "            \n",
    "            # Access RGB values here\n",
    "            r_mean = np.mean(rgb_arr[:,:,0])\n",
    "            g_mean = np.mean(rgb_arr[:,:,1])\n",
    "            b_mean = np.mean(rgb_arr[:,:,2])\n",
    "#             a_mean = np.mean(rgb_arr[:,:,3])\n",
    "            \n",
    "            r_max = np.amax(rgb_arr[:,:,0])\n",
    "            g_max = np.amax(rgb_arr[:,:,1])\n",
    "            b_max = np.amax(rgb_arr[:,:,2])\n",
    "#             a_max = max(rgb_arr[:,:,3])\n",
    "            \n",
    "            r_min = np.amin(rgb_arr[:,:,0])\n",
    "            g_min = np.amin(rgb_arr[:,:,1])\n",
    "            b_min = np.amin(rgb_arr[:,:,2])\n",
    "#             a_min = min(rgb_arr[:,:,3])\n",
    "            \n",
    "            r_median = np.median(rgb_arr[:,:,0])\n",
    "            g_median = np.median(rgb_arr[:,:,1])\n",
    "            b_median = np.median(rgb_arr[:,:,2])\n",
    "#             a_median = np.median(rgb_arr[:,:,3])\n",
    "            \n",
    "            r_std = np.std(rgb_arr[:,:,0])\n",
    "            g_std = np.std(rgb_arr[:,:,1])\n",
    "            b_std = np.std(rgb_arr[:,:,2])\n",
    "#             a_std = np.std(rgb_arr[:,:,3])\n",
    "            \n",
    "            color_df.loc[len(color_df)] = [filepath.split('/')[-1],r_mean,g_mean,b_mean\n",
    "                                           ,r_max,g_max,b_max\n",
    "                                          ,r_min,g_min,b_min\n",
    "                                          ,r_median,g_median,b_median\n",
    "                                          ,r_std,g_std,b_std]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Export dataframe to csv so that we can merge later\n",
    "color_df.to_csv('daylights1.csv',index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2. Merge daytime images with DHS data\n",
    "\n",
    "- **INPUT**: \n",
    " - `google_image_features_basic.csv`: Satellite imagery features, from Step 5.1\n",
    " - `rwanda_cluster_avg_asset_2010.csv`: DHS cluster averages, from Step 2\n",
    "- **OUTPUT**: Merged datasets\n",
    " - `data/model/DHS_daytime.csv`: Merged dataset with 492 rows, and 16 columns (one indicates average cluster wealth, 15 daytime image features)\n",
    "\n",
    "Now that you have feature vectors for each image, you should merge these with the DHS data indicated average cluster wealth. Follow a similar procedure as you did with 3.1, i.e., determine which image feature vectors are associated with each cluster, and then calculate, for each cluster, the average value of each feature.\n",
    "\n",
    "Save your output as `DHS_daytime.csv` and check that it is roughly the same as the file we have provided. There may be slight differences if you chose to calculate a different set of features than those described in 5.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dhs_df = pd.read_csv('rwanda_cluster_avg_asset_2010.csv')\n",
    "daytime_df = pd.read_csv('daylights1.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def extract_row_cols(name,dim):\n",
    "    name = name.split('.')[0]\n",
    "    if dim == 'col':\n",
    "        return int(name.split('_')[0])\n",
    "    if dim == 'row':\n",
    "        return int(name.split('_')[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "daytime_df['row'] = daytime_df.apply(lambda item: extract_row_cols(item.filename,'row'),axis=1)\n",
    "daytime_df['col'] = daytime_df.apply(lambda item: extract_row_cols(item.filename,'col'),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def avg_feature_cluster(lat,lon,day_df,feature,top_left_x_coords=npzfile['top_left_x_coords']\n",
    "                        ,top_left_y_coords=npzfile['top_left_y_coords'],bands_data=npzfile['bands_data']):\n",
    "    \n",
    "    lon_idx,lat_idx = get_cell_idx(lon,lat,top_left_x_coords,top_left_y_coords)\n",
    "    \n",
    "    return np.mean(day_df[day_df['row'].between(lon_idx-5,lon_idx+5)\n",
    "                          & day_df['col'].between(lat_idx-5,lat_idx+5)][feature])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def merge_daytime_dhs(day_df,night_df):\n",
    "    \n",
    "    night_df['mean_r_mean'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                              ,item.longitude,day_df,'r_mean'),axis=1)\n",
    "    night_df['mean_g_mean'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                              ,item.longitude,day_df,'g_mean'),axis=1)\n",
    "    night_df['mean_b_mean'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                              ,item.longitude,day_df,'b_mean'),axis=1)\n",
    "    night_df['mean_r_max'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'r_max'),axis=1)\n",
    "    night_df['mean_g_max'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'g_max'),axis=1)\n",
    "    night_df['mean_b_max'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'b_max'),axis=1)\n",
    "    night_df['mean_r_min'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'r_min'),axis=1)\n",
    "    night_df['mean_g_min'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'g_min'),axis=1)    \n",
    "    night_df['mean_b_min'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'b_min'),axis=1)\n",
    "    night_df['mean_r_median'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                                ,item.longitude,day_df,'r_median'),axis=1)\n",
    "    night_df['mean_g_median'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                                ,item.longitude,day_df,'g_median'),axis=1)    \n",
    "    night_df['mean_b_median'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                                ,item.longitude,day_df,'b_median'),axis=1)\n",
    "    night_df['mean_r_std'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'r_std'),axis=1)\n",
    "    night_df['mean_g_std'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'g_std'),axis=1)    \n",
    "    night_df['mean_b_std'] = night_df.apply(lambda item: avg_feature_cluster(item.latitude\n",
    "                                                                             ,item.longitude,day_df,'b_std'),axis=1)    \n",
    "    # So on for other features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "merge_daytime_dhs(daytime_df,dhs_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3. Fit a model of wealth as a function of basic daytime features\n",
    "- **INPUT**: \n",
    " - `data/model/DHS_daytime.csv`, from Step 5.2\n",
    "- **OUTPUT**: \n",
    " - R^2 of model\n",
    " \n",
    "As in 3.2, use 10-fold cross-validation to fit a model of cluster-level DHS wealth (your response/dependent variable) as a function of the nightlights data (your predictor/independent variables). Since you have a reasonably large number of predictor variables, you should use a model that incorporates some form of regularization (e.g., ridge regression, lasso regression, or a tree-based method).  Report the cross-validated R^2 of your model (i.e., the average R^2 of your 10 test folds).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Extract features from daytime imagery using deep learning libraries\n",
    "\n",
    "This is where things get interesting. You will use existing libraries to extract more meaningful features from the daytime imagery, similar to what is shown in Fig. 2 of the paper.\n",
    "\n",
    "## 6.1. Use the keras library to use a basic CNN to extract features of the daytime images \n",
    " \n",
    "- **INPUT**: \n",
    " - `google_image/...`: Raw images, from Step 4\n",
    "- **OUTPUT**: \n",
    " - `google_image_features_cnn.csv`: Image features \n",
    "\n",
    "Begin by using a Convolutional Neural Network that has been pre-trained on ImageNet to extract features from the images. We recommend using the [`Keras` library](https://keras.io/), which provides a very straightforward interface to [TensorFlow](https://www.tensorflow.org/).\n",
    "\n",
    "Hints:\n",
    " - This [short intro](https://github.com/fchollet/deep-learning-models/blob/master/README.md) will help you get started with extracting features from the CNN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K.image_dim_ordering: tf\n"
     ]
    }
   ],
   "source": [
    "# import the necessary packages\n",
    "from keras.preprocessing import image as image_utils\n",
    "from imagenet_utils import decode_predictions\n",
    "from imagenet_utils import preprocess_input\n",
    "from vgg16 import VGG16\n",
    "import numpy as np\n",
    "\n",
    "model = VGG16(weights='imagenet', include_top=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# model.layers.pop()\n",
    "# model.summary(line_length=150)\n",
    "from keras.models import Model\n",
    "# We need to access the output of layer 'fc2'\n",
    "layer_name = 'fc2'\n",
    "intermediate_layer_model = Model(input=model.input,\n",
    "                                 output=model.get_layer(layer_name).output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 4096)\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "from keras.preprocessing import image\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "columns = [i for i in range(0,4096)]\n",
    "columns.append('filename')\n",
    "\n",
    "cnn_features = pd.DataFrame(columns=columns)\n",
    "\n",
    "for root, dirs, files in os.walk('google_image/4'):\n",
    "    for name in files:\n",
    "        if name.endswith('.jpg'):\n",
    "            file_path = os.path.join(root, name)\n",
    "            img = image.load_img(file_path, target_size=(224, 224))\n",
    "            x = image.img_to_array(img)\n",
    "            x = np.expand_dims(x, axis=0)\n",
    "            x = preprocess_input(x)\n",
    "\n",
    "            features = intermediate_layer_model.predict(x)\n",
    "            print(features.shape)\n",
    "            # Convert numpy array into dataframe\n",
    "            arr_to_df = pd.DataFrame(features)\n",
    "            arr_to_df['filename'] = name\n",
    "            \n",
    "            cnn_features = cnn_features.append(arr_to_df,ignore_index=True)\n",
    "            \n",
    "            break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cnn_features.to_csv('foo.csv',index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6.2. Test whether these new features of satellite imagery can predict wealth\n",
    "- **INPUT**: \n",
    " - `google_image_features_cnn.csv`: Satellite imagery features, from Step 6.1\n",
    " - `rwanda_cluster_avg_asset_2010.csv`: DHS cluster averages, from Step 2\n",
    "- **OUTPUT**: \n",
    " - `data/model/DHS_daytime.csv`: Merged dataset with 492 rows, and 4097 columns (one indicates average cluster wealth, 4096 CNN-based features)\n",
    " - R^2 of model\n",
    " \n",
    "Calculate the average value of each feature for each of the DHS clusters. As in Step 3.1 and 5.2, you will want to aggregate over images near the cluster centroid by taking the average value for each feature. Create a scatterplot showing the relationship between average cluster wealth (y-axis) and the first principal component of all of your image features (x-axis) - in other words, run PCA on your 4096 image features and plot the first PC on the x-axis. Your scatterplot should have one dot for each of the 492 DHS clusters.\n",
    "\n",
    "Use 10-fold cross-validation to fit a model of cluster-level DHS wealth (your response/dependent variable) as a function of the \"deep\" features (your predictor/independent variables). Use a model that incorporates some form of regularization (e.g., ridge regression, lasso regression, or a tree-based method).  Report the cross-validated R^2 of your model (i.e., the average R^2 of your 10 test folds)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "cnn_features_df = pd.read_csv('foo.csv')\n",
    "dhs_avgs_df = pd.read_csv('rwanda_cluster_avg_asset_2010.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cnn_features_df['row'] = cnn_features_df.apply(lambda item: extract_row_cols(item.filename,'row'),axis=1)\n",
    "cnn_features_df['col'] = cnn_features_df.apply(lambda item: extract_row_cols(item.filename,'col'),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def merge_cnn_dhs(cnn,dhs):\n",
    "    dhs['mean_0'] = dhs.apply(lambda item: avg_feature_cluster(item.latitude,item.longitude,cnn_features_df,'0'),axis=1)\n",
    "    \n",
    "    return dhs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "(0, 'occurred at index 0')",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/indexes/base.py\u001b[0m in \u001b[0;36mget_loc\u001b[0;34m(self, key, method, tolerance)\u001b[0m\n\u001b[1;32m   2133\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2134\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2135\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/index.pyx\u001b[0m in \u001b[0;36mpandas.index.IndexEngine.get_loc (pandas/index.c:4443)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/index.pyx\u001b[0m in \u001b[0;36mpandas.index.IndexEngine.get_loc (pandas/index.c:4289)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/src/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas.hashtable.PyObjectHashTable.get_item (pandas/hashtable.c:13733)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/src/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas.hashtable.PyObjectHashTable.get_item (pandas/hashtable.c:13687)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 0",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-60-ead679996c9a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmerge_cnn_dhs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcnn_features_df\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdhs_avgs_df\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-43-3eda50d523bb>\u001b[0m in \u001b[0;36mmerge_cnn_dhs\u001b[0;34m(cnn, dhs)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mmerge_cnn_dhs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcnn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdhs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0mdhs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mean_0'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdhs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mavg_feature_cluster\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlatitude\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlongitude\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcnn_features_df\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mdhs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36mapply\u001b[0;34m(self, func, axis, broadcast, raw, reduce, args, **kwds)\u001b[0m\n\u001b[1;32m   4161\u001b[0m                     \u001b[0;32mif\u001b[0m \u001b[0mreduce\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4162\u001b[0m                         \u001b[0mreduce\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4163\u001b[0;31m                     \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_apply_standard\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreduce\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mreduce\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4164\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4165\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_apply_broadcast\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_apply_standard\u001b[0;34m(self, func, axis, ignore_failures, reduce)\u001b[0m\n\u001b[1;32m   4257\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4258\u001b[0m                 \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseries_gen\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4259\u001b[0;31m                     \u001b[0mresults\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4260\u001b[0m                     \u001b[0mkeys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4261\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-43-3eda50d523bb>\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(item)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mmerge_cnn_dhs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcnn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdhs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0mdhs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'mean_0'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdhs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mlambda\u001b[0m \u001b[0mitem\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mavg_feature_cluster\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlatitude\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlongitude\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcnn_features_df\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mdhs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-59-9aef3db86299>\u001b[0m in \u001b[0;36mavg_feature_cluster\u001b[0;34m(lat, lon, day_df, feature, top_left_x_coords, top_left_y_coords, bands_data)\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m     return np.mean(day_df[day_df['row'].between(lon_idx-5,lon_idx+5)\n\u001b[0;32m---> 10\u001b[0;31m                           & day_df['col'].between(lat_idx-5,lat_idx+5)][feature])\n\u001b[0m",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   2057\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_multilevel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2058\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2059\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2060\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2061\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_getitem_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_getitem_column\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   2064\u001b[0m         \u001b[0;31m# get column\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2065\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_unique\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2066\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_item_cache\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2067\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2068\u001b[0m         \u001b[0;31m# duplicate columns & possible reduce dimensionality\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36m_get_item_cache\u001b[0;34m(self, item)\u001b[0m\n\u001b[1;32m   1384\u001b[0m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1385\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mres\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1386\u001b[0;31m             \u001b[0mvalues\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_data\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1387\u001b[0m             \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_box_item_values\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalues\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1388\u001b[0m             \u001b[0mcache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/core/internals.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self, item, fastpath)\u001b[0m\n\u001b[1;32m   3539\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3540\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misnull\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3541\u001b[0;31m                 \u001b[0mloc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mitem\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3542\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3543\u001b[0m                 \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0misnull\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/Proxima/anaconda/lib/python3.5/site-packages/pandas/indexes/base.py\u001b[0m in \u001b[0;36mget_loc\u001b[0;34m(self, key, method, tolerance)\u001b[0m\n\u001b[1;32m   2134\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2135\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2136\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_cast_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2137\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2138\u001b[0m         \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtolerance\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtolerance\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/index.pyx\u001b[0m in \u001b[0;36mpandas.index.IndexEngine.get_loc (pandas/index.c:4443)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/index.pyx\u001b[0m in \u001b[0;36mpandas.index.IndexEngine.get_loc (pandas/index.c:4289)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/src/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas.hashtable.PyObjectHashTable.get_item (pandas/hashtable.c:13733)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/src/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas.hashtable.PyObjectHashTable.get_item (pandas/hashtable.c:13687)\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: (0, 'occurred at index 0')"
     ]
    }
   ],
   "source": [
    "result = merge_cnn_dhs(cnn_features_df,dhs_avgs_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.5328179999999998"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lon_idx,lat_idx = get_cell_idx(dhs_avgs_df.iloc[1]['longitude'],dhs_avgs_df.iloc[0]['latitude']\n",
    "                               ,npzfile['top_left_x_coords'],npzfile['top_left_y_coords'])\n",
    "    \n",
    "#     row_cond = day_df['row'] in range(lon_idx-5,lon_idx+5)\n",
    "#     col_cond = day_df['col'] in range(lat_idx-5,lat_idx+5)\n",
    "    \n",
    "#     return np.mean(day_df[(day_df['row'].between(lon_idx-5,lon_idx+5)) \n",
    "#                           & (day_df['col'].between(lat_idx-5,lat_idx+5)),0][[feature]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "nan"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(cnn_features_df[cnn_features_df['row'].between(lon_idx-5,lon_idx+5) & cnn_features_df['col'].between(lat_idx-5,lat_idx+5)]['1'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}