FPARPAR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from pyhdf.SD import SD, SDC\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import netCDF4 as nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.getcwd()\n",
    "path1 = 'E:\\\\SIF data\\\\FPAR_2016'\n",
    "path2 = 'E:\\\\SIF data\\\\GOME_2016\\\\GOME_2016\\\\'\n",
    "os.chdir('E:\\\\SIF data\\\\FPAR_2016')\n",
    "# os.chdir('E:\\\\SIF data\\\\GOME_2016\\\\GOME_2016\\\\')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import my_functions\n",
    "my_functions.coords_to_ind(-90,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ILLINOIS = (36.7281, -91.8524, 42.493703, -87.020001)\n",
    "ILLINOIS = (30, -120, 49, -80)\n",
    "arr_ = []\n",
    "\n",
    "\n",
    "def find_max_qc(ILLINOIS):\n",
    "    \n",
    "    qc_data = None\n",
    "    fpar_data = None\n",
    "    prefix ='FPAR_A2016'\n",
    "    suffix = '.hdf'\n",
    "    max_ = -np.inf\n",
    "    acc = []\n",
    "    for i in range(1,361,8):\n",
    "        a = str(int(i))\n",
    "        if i < 10:\n",
    "            a = '00'+ a\n",
    "        elif i < 100:\n",
    "            a = '0' + a\n",
    "\n",
    "        query = prefix + a + suffix\n",
    "        try:\n",
    "            data = SD(query, SDC.READ)\n",
    "            print('read')\n",
    "            cur_fpar = data.select('Fpar_500m')[:]\n",
    "            cur_qc = data.select('FparExtra_QC')[:]\n",
    "            left, down = my_functions.coords_to_ind(ILLINOIS[1], ILLINOIS[0])\n",
    "            right, up = my_functions.coords_to_ind(ILLINOIS[3], ILLINOIS[2])\n",
    "            print(right, left)\n",
    "            print(down, up)\n",
    "            qcs = cur_qc[up:down, left:right].flatten()\n",
    "            ratio = my_functions.good_qc(qcs)\n",
    "            acc.append([ratio, query])\n",
    "#             print(ratio)\n",
    "#             if ratio > max_:\n",
    "#                 fpar_data = cur_fpar\n",
    "#                 qc_data = cur_qc\n",
    "#                 max_ = ratio\n",
    "\n",
    "        except:\n",
    "            continue\n",
    "        \n",
    "    return fpar_data, qc_data, acc\n",
    "        \n",
    "        \n",
    "\n",
    "def run_fpar_workflow(ILLINOIS, fpar_data, qc_data):\n",
    "    arr_ = []\n",
    "    \n",
    "    for i in np.arange(ILLINOIS[0], ILLINOIS[2], 0.1):\n",
    "        for j in np.arange(ILLINOIS[1], ILLINOIS[3], 0.1):\n",
    "            cur_lat, cur_lon = i,j\n",
    "            \n",
    "            lon_ind, lat_ind = my_functions.coords_to_ind(cur_lon, cur_lat)\n",
    "            fp, qc = fpar_data[lat_ind, lon_ind], qc_data[lat_ind, lon_ind]\n",
    "            if my_functions.get_cloud(qc)[0] == '0' and my_functions.get_cloud(qc)[1] == '0':\n",
    "                arr_.append([cur_lat, cur_lon, fp])\n",
    "\n",
    "    return arr_\n",
    "\n",
    "def find_fpar_date(fpar):\n",
    "    prefix ='FPAR_A2016'\n",
    "    suffix = '.hdf'\n",
    "    suspects = []\n",
    "    for i in range(1,361,8):\n",
    "        a = str(int(i))\n",
    "        if i < 10:\n",
    "            a = '00'+ a\n",
    "        elif i < 100:\n",
    "            a = '0' + a\n",
    "\n",
    "        query = prefix + a + suffix\n",
    "        try:\n",
    "            data = SD(query, SDC.READ)\n",
    "            print('read')\n",
    "            cur_fpar = data.select('Fpar_500m')[:]\n",
    "            cur_qc = data.select('FparExtra_QC')[:]\n",
    "            if np.mean(cur_fpar) == np.mean(fpar):\n",
    "                print(query)\n",
    "                suspects.append(query)\n",
    "        except:\n",
    "            continue\n",
    "        \n",
    "    return suspects\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fp_dat, qc_dat,acc = find_max_qc(ILLINOIS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cleaned_acc = sorted(acc, key = lambda x: x[0], reverse = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cleaned_acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filtered_fpar = run_fpar_workflow(ILLINOIS, a,b)\n",
    "print(filtered_fpar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dat = SD('FPAR_A2016225.hdf', SDC.READ)\n",
    "\n",
    "fpp2 = dat.select('Fpar_500m')[:]\n",
    "qc2 = dat.select('FparExtra_QC')[:]\n",
    "print(fpp2.shape)\n",
    "arr2 = run_fpar_workflow(ILLINOIS, fpp2, qc2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filtered_data = np.array(arr2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "suspects = find_fpar_date(a)\n",
    "print(suspects)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filtered_data = np.array(filtered_fpar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filtered_data[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from dbfread import DBF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sifutil\n",
    "ILLINOIS = (36.7281, -91.8524, 42.493703, -87.020001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><b>Functions for processing fpar data, like get cloud, convert the qc index</b></p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy.linalg as la\n",
    "\n",
    "def convert_binary(num):\n",
    "    str_ = \"{0:b}\".format(num)\n",
    "    if len(str_) < 8:\n",
    "        str_ = '0'*(8-len(str_)) + str_\n",
    "        \n",
    "    return str_\n",
    "        \n",
    "def get_cloud(num):\n",
    "    str_ = convert_binary(num)\n",
    "    return str_[1], str_[2]\n",
    "\n",
    "\n",
    "def interpolation(x,y):\n",
    "    x = np.array(x)\n",
    "    matrix = np.array([x**i for i in range(len(x))]).transpose()\n",
    "    print(matrix)\n",
    "    coeffs = la.solve(matrix,y)\n",
    "    return coeffs\n",
    "\n",
    "def get_smooth_line(x,y):\n",
    "    coeffs = interpolation(x,y)\n",
    "    x_values = np.linspace(min(x), max(x), 100)\n",
    "    y_values = []\n",
    "    for i in x_values:\n",
    "        value = 0\n",
    "        for j in range(len(coeffs)):\n",
    "            value += coeffs[j]*i**j\n",
    "            \n",
    "        y_values.append(value)\n",
    "    return [list(x_values), y_values]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "time_series = []\n",
    "qc = []\n",
    "x_values = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>function for getting the whole fpar series for a year of a certain location (lat,lon)</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "def get_fpar_series(lat, lon):\n",
    "\n",
    "    time_series = []\n",
    "    qc = []\n",
    "    x_values = []\n",
    "    prefix ='FPAR_A2016'\n",
    "    suffix = '.hdf'\n",
    "\n",
    "    print(prefix+suffix)\n",
    "\n",
    "    for i in range(1,9,8):\n",
    "        a = str(int(i))\n",
    "        if i < 10:\n",
    "            a = '00'+ a\n",
    "        elif i < 100:\n",
    "            a = '0' + a\n",
    "\n",
    "        query = prefix + a + suffix\n",
    "        try:\n",
    "            data = SD(query, SDC.READ)\n",
    "            FPAR_data = data.select('Fpar_500m')[:]\n",
    "            QC_data = data.select('FparExtra_QC')[:]\n",
    "\n",
    "    #         print(QC_data.shape)\n",
    "            Q = QC_data[lat,lon]\n",
    "            time_series.append(FPAR_data[lat,lon])\n",
    "            x_values.append(i)\n",
    "            qc.append(Q)\n",
    "            print(time_series)\n",
    "            print(x_values)\n",
    "            print(qc)\n",
    "        except:\n",
    "            continue\n",
    "\n",
    "    return time_series, x_values, qc\n",
    "    \n",
    "    \n",
    "\n",
    "# data = SD('FPAR_A2016361.hdf', SDC.READ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>function for filtering the fpars with cloud, given qc arrays, x values and fpar series</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "def clean_data(time_series, x_values, qc):\n",
    "    good_fpars = []\n",
    "    good_dates = []\n",
    "    for i in range(len(time_series)):\n",
    "        if get_cloud(qc[i])[0] == '0' and get_cloud(qc[i])[1] == '0':\n",
    "            good_fpars.append(time_series[i])\n",
    "            good_dates.append(x_values[i])\n",
    "            \n",
    "    return good_dates, good_fpars\n",
    "\n",
    "def good_qc(qc):\n",
    "    cnt = 0\n",
    "    for i in range(len(qc)):\n",
    "        if get_cloud(qc[i])[0] == '0' and get_cloud(qc[i])[1] == '0':\n",
    "            cnt += 1\n",
    "            \n",
    "    return cnt/len(qc)\n",
    "#     return cnt/len(qc)\n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.plot(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = SD('FPAR_A2016225.hdf', SDC.READ).select('Fpar_500m')[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import cv2\n",
    "cv2.imwrite('225_august_fpar.jpg', data, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# from pyhdf.SD import SD, SDC\n",
    "# import sifutil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "coords_to_ind(135, 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import math\n",
    "lat_eq = 360 / 40075\n",
    "lon_eq = math.cos(self.lat_0) * 360 / 40057.\n",
    "lon_a = self.lon_0 - ((self.side_len/2) * lon_eq)\n",
    "lat_a = self.lat_0 - ((self.side_len/2) * lat_eq)\n",
    "unit = self.side_len / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from io import StringIO\n",
    "import re\n",
    "# import requests\n",
    "# from bs4 import BeautifulSoup\n",
    "# import pandas as pd\n",
    "# import warnings\n",
    "# BASE_CDL_URL = 'https://nassgeodata.gmu.edu/axis2/services/CDLService/GetCDLStat'\n",
    "# CHAMPAIGN = 17019\n",
    "\n",
    "\n",
    "# def from_csv(filepath):\n",
    "#     df = pd.read_csv(filepath)\n",
    "#     return df\n",
    "\n",
    "\n",
    "# def get_by_fips(year, fips):\n",
    "#     '''Grab CDL data for a county by FIPS code'''\n",
    "#     url = BASE_CDL_URL + '?year=' + str(year) + '&fips=' + str(fips) + \"&format=csv\"\n",
    "#     res = requests.get(url, verify = False)\n",
    "#     returnurl = BeautifulSoup(res.text, 'lxml').find('returnurl').text\n",
    "#     #print(returnurl)\n",
    "#     rawdata = requests.get(returnurl, verify = False).text\n",
    "#     raw_iter = StringIO(rawdata)\n",
    "#     df = pd.read_csv(raw_iter, sep=\" *, * \", engine='python')\\\n",
    "#            .apply(pd.to_numeric, errors='ignore')\\\n",
    "#            .set_index(\"Category\")\n",
    "#     return df\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# def get_by_box(year, llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat):\n",
    "#     '''Grab CDL data by bounding box'''\n",
    "#     x1, y1 = sifutil.convertProjection(llcrnrlon, llcrnrlat, sifutil.WGS84, sifutil.CONUS_ALBERS)\n",
    "#     x2, y2 = sifutil.convertProjection(urcrnrlon, urcrnrlat, sifutil.WGS84, sifutil.CONUS_ALBERS)\n",
    "#     url = BASE_CDL_URL + '?year=' + str(year) + '&bbox=' + str(min(x1,x2)) + \",\" +\\\n",
    "#           str(min(y1, y2)) + \",\" + str(max(x1, x2)) + \",\" + str(max(y1, y2)) + \"&format=csv\"\n",
    "#     print(url)\n",
    "\n",
    "#     with warnings.catch_warnings():\n",
    "#         warnings.simplefilter(\"ignore\")\n",
    "#         res = requests.get(url, verify = False)\n",
    "#         returnurl = BeautifulSoup(res.text, 'lxml').find('returnurl').text\n",
    "#     #print(returnurl)\n",
    "#     with warnings.catch_warnings():\n",
    "#         rawdata = requests.get(returnurl, verify = False).text\n",
    "#     raw_iter = StringIO(rawdata)\n",
    "#     df = pd.read_csv(raw_iter, sep=\" *, * \")\\\n",
    "#            .apply(pd.to_numeric, errors='ignore')\\\n",
    "#            .set_index(\"Category\")\n",
    "#     return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getNC_GOME(filepath):\n",
    "    data =  nc.Dataset(filepath, mode='r')\n",
    "    lons = data['longitude'][:]\n",
    "    lats = data['latitude'][:]\n",
    "    sif = data['SIF_740'][:]\n",
    "    data.close()\n",
    "    return (lons, lats, sif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ILLINOIS = (36.7281, -91.8524, 42.493703, -87.020001)\n",
    "bounding_box = (40.16, -89, 41, -87.5)\n",
    "I_STATES_COORDS = (37, -96.65, 43.5, -84.81)\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_illinois(ILLINOIS):\n",
    "    fig, ax1 = plt.subplots()\n",
    "    m = Basemap(projection='merc', llcrnrlat=ILLINOIS[0], llcrnrlon=ILLINOIS[1],\\\n",
    "                urcrnrlat=ILLINOIS[2], urcrnrlon=ILLINOIS[3], resolution='c', ax=ax1)\n",
    "    m.drawstates()\n",
    "    m.drawcoastlines()\n",
    "\n",
    "    m.fillcontinents(color='beige', lake_color='aqua')\n",
    "    m.drawparallels(np.arange(36, 43, 1), labels=[0,1,1,0])\n",
    "    m.drawmeridians(np.arange(-92, -87, 1), labels=[1,0,0,1])\n",
    "    ax1.set_title('poop')\n",
    "\n",
    "    return ax1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plot_illinois(ILLINOIS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>A whole process to obtain the GNOME data</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir('E:\\\\SIF data\\\\GOME_2016\\\\GOME_2016\\\\')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import netCDF4 as nc\n",
    "gome_0711 = getNC_GOME('08\\\\ret_f_nr5_nsvd12_v26_waves734_nolog.20160801_v27_all.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(ILLINOIS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def inIllinois(ILLINOIS):\n",
    "    if "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cur_data = NetCDFFile('08\\\\ret_f_nr5_nsvd12_v26_waves734_nolog.20160801_v27_all.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# cur_data.variables['Latitude_corners'][:]\n",
    "# plt.hist(cur_data.variables['Quality_flag'][:])\n",
    "import math\n",
    "\n",
    "lat_corners = cur_data.variables['Latitude_corners'][:]\n",
    "lon_corners = cur_data.variables['Longitude_corners'][:]\n",
    "arg_min = None\n",
    "\n",
    "for i in range(len(lat_corners)):\n",
    "    if math.fabs(lat_corners[i][0]-40) <= 2 and math.fabs(lon_corners[i][0] + 88) <= 2:\n",
    "        arg_min = i\n",
    "        \n",
    "\n",
    "print(lat_corners[arg_min], lon_corners[arg_min])\n",
    "print(len(lat_corners), len(lon_corners))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "averaged_sif = cur_data.variables['Daily_averaged_SIF'][:]\n",
    "lons = cur_data.variables['longitude'][:]\n",
    "lats = cur_data.variables['latitude'][:]\n",
    "clouds = cur_data.variables['cloud_fraction'][:]\n",
    "qa = cur_data.variables['Quality_flag'][:]\n",
    "print(len(lons))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(lat_corners[arg_min], lon_corners[arg_min])\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.scatter(lat_corners[arg_min], lon_corners[arg_min])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "arr = get_Illinois(lats, lons, averaged_sif, ILLINOIS, clouds, qa)\n",
    "lat, lon = arr[:,0], arr[:,1]\n",
    "sifs = arr[:,2]\n",
    "print(arr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from numpy import array\n",
    "\n",
    "plt.figure(figsize = (20,10))\n",
    "m = Basemap(projection='merc', llcrnrlat=0, llcrnrlon=49,\\\n",
    "                urcrnrlat=-, urcrnrlon=ILLINOIS[3], resolution='i')\n",
    "m.drawstates()\n",
    "m.drawcoastlines()\n",
    "# m.fillcontinents(color='beige', lake_color='aqua')\n",
    "m.drawparallels(np.arange(36, 43, 1), labels=[0,1,1,0])\n",
    "m.drawmeridians(np.arange(-92, -87, 1), labels=[0,1,1,0])\n",
    "\n",
    "\n",
    "l = [-87]\n",
    "la = [40]\n",
    "s = [1]\n",
    "x,y = m(lon,lat)\n",
    "m.hexbin(x,y,C = sifs, gridsize = 100)\n",
    "\n",
    "\n",
    "m.colorbar(location='bottom')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# from mpl_toolkits.basemap import Basemap\n",
    "# import matplotlib.pyplot as plt\n",
    "# import matplotlib.colors as colors\n",
    "# from numpy import array\n",
    "# from numpy import max\n",
    "\n",
    "# # map = Basemap(llcrnrlon=110,llcrnrlat=0,urcrnrlon=140.,urcrnrlat=40.,\n",
    "# #              resolution='i', projection='merc')\n",
    "# map = Basemap(projection='merc', llcrnrlat=ILLINOIS[0], llcrnrlon=ILLINOIS[1],\\\n",
    "#                 urcrnrlat=ILLINOIS[2], urcrnrlon=ILLINOIS[3], resolution='i')\n",
    "\n",
    "\n",
    "# lats = [5,6,7,39,9,10,11,12,13,14,31.17]\n",
    "# lons = [121,121,121,-88,121,121,121,121,121,121,121.5]\n",
    "# c = [1,1.25,1.5,1.75,2,2.25,2.50,2.75,3,3,2.5]\n",
    "\n",
    "\n",
    "# la = [39, 40, 41]\n",
    "# lo = [-88, -87, -88]\n",
    "# c = [2,3,3]\n",
    "# x, y = map(lo, la)\n",
    "\n",
    "\n",
    "# map.drawcoastlines()\n",
    "# map.drawstates()\n",
    "\n",
    "# map.hexbin(array(x), array(y), C =c, gridsize = 8)\n",
    "\n",
    "# map.colorbar(location='bottom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "ILLINOIS = (36.7281, -91.8524, 42.493703, -87.020001)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(lats)\n",
    "plt.plot(lons)\n",
    "plt.plot(x)\n",
    "print(lats[118000], lons[118000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "def get_Illinois(lats, lons, sifs, ILLINOIS, clouds, qa):\n",
    "    illinois_ = []\n",
    "    for i in range(len(lons)):\n",
    "        if lons[i] >= ILLINOIS[1] and qa[i] == 2 and lons[i] <= ILLINOIS[3] and lats[i] >= ILLINOIS[0] and lats[i] <= ILLINOIS[2] and clouds[i] <= 0.2:\n",
    "            illinois_.append([lats[i], lons[i], sifs[i]])\n",
    "            \n",
    "    return np.array(illinois_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# print(illinois_gome_lons)\n",
    "# print(illinois_gome_lats)\n",
    "\n",
    "arr = get_Illinois(lats, lons, averaged_sif, bounding_box,clouds)\n",
    "print(arr.shape)\n",
    "# print(arr[:,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x,y = clean_data(fpar_series, x_values, qc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(x,y)\n",
    "plt.figure()\n",
    "plt.plot(x_values, fpar_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a_lot_of_sif = []\n",
    "x_ = []\n",
    "for i in range(1,12):\n",
    "    for j in range(1,31):\n",
    "#         print(i,j)\n",
    "        if i < 10:\n",
    "            zero = '0'\n",
    "        else:\n",
    "            zero = ''\n",
    "            \n",
    "        if j < 10:\n",
    "            zero2 = '0'\n",
    "        else:\n",
    "            zero2 = ''\n",
    "            \n",
    "        \n",
    "        date = '2016' + zero + str(i) + zero2 + str(j)\n",
    "        dir_ = zero + str(i) + '\\\\\\\\'\n",
    "        file = dir_ + 'ret_f_nr5_nsvd12_v26_waves734_nolog.' + date + '_v27_all.nc'\n",
    "        \n",
    "        try: \n",
    "            read_in = NetCDFFile(file)\n",
    "            data = getNC_GOME(file)\n",
    "            lons = data[0]\n",
    "            lats = data[1]\n",
    "            sifs = data[2]\n",
    "            clouds =  read_in['cloud_fraction'][:]\n",
    "            qa = read_in['Quality_flag'][:]\n",
    "            processed_data = get_Illinois(lats, lons, sifs, bounding_box, clouds, qa)\n",
    "            print('good')\n",
    "            if processed_data != []:\n",
    "                a_lot_of_sif.append(np.mean(processed_data[:,2]))\n",
    "                print(i,j)\n",
    "                x_.append((i-1)*30+j)\n",
    "        \n",
    "            \n",
    "        except:\n",
    "            continue\n",
    "#         gome_0711 = getNC_GOME('05\\\\ret_f_nr5_nsvd12_v26_waves734_nolog.20160501_v27_all.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(len(a_lot_of_sif))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(x_, a_lot_of_sif)\n",
    "plt.title('filtered sif plot with at most 20% cloud coverage and the best qa')\n",
    "plt.xlabel('days')\n",
    "plt.ylabel('sif value')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "a = np.array([[1,1,2],[2,2,1],[3,1,1]])\n",
    "plt.imshow(a, label = ['a','b','c'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "colors = [ im.cmap(im.norm(value)) for value in values]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "map = Basemap(projection='ortho', \n",
    "              lat_0=0, lon_0=0)\n",
    "\n",
    "map.drawmapboundary(fill_color='aqua')\n",
    "map.fillcontinents(color='coral',lake_color='aqua')\n",
    "map.drawcoastlines()\n",
    "\n",
    "lons = [0, 10, -30, -20]\n",
    "lats = [0, -10, 40, -20]\n",
    "\n",
    "x, y = map(lons, lats)\n",
    "\n",
    "map.scatter(x, y, marker='D',color='m')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import animation\n",
    "\n",
    "# First set up the figure, the axis, and the plot element we want to animate\n",
    "fig = plt.figure()\n",
    "ax = plt.axes(xlim=(0, 2), ylim=(-2, 2))\n",
    "line, = ax.plot([], [], lw=2)\n",
    "\n",
    "# initialization function: plot the background of each frame\n",
    "def init():\n",
    "    line.set_data([], [])\n",
    "    return line,\n",
    "\n",
    "# animation function.  This is called sequentially\n",
    "def animate(i):\n",
    "    x = np.linspace(0, 2, 1000)\n",
    "    y = np.sin(2 * np.pi * (x - 0.01 * i))\n",
    "    line.set_data(x, y)\n",
    "    return line,\n",
    "\n",
    "# call the animator.  blit=True means only re-draw the parts that have changed.\n",
    "anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                               frames=200, interval=20, blit=True)\n",
    "\n",
    "# save the animation as an mp4.  This requires ffmpeg or mencoder to be\n",
    "# installed.  The extra_args ensure that the x264 codec is used, so that\n",
    "# the video can be embedded in html5.  You may need to adjust this for\n",
    "# your system: for more information, see\n",
    "# http://matplotlib.sourceforge.net/api/animation_api.html\n",
    "anim.save('basic_animation.html', fps=30, extra_args=['-vcodec', 'libx264'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filtered_data[:,1][0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.basemap import Basemap\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "from numpy import array\n",
    "# from numpy import max\n",
    "\n",
    "plt.figure(figsize = (20,10))\n",
    "# map = Basemap(llcrnrlon=110,llcrnrlat=0,urcrnrlon=140.,urcrnrlat=40.,\n",
    "#              resolution='i', projection='merc')\n",
    "map = Basemap(projection='aea', llcrnrlat=10, llcrnrlon=-130,\\\n",
    "                urcrnrlat=49, urcrnrlon = -80, resolution='i')\n",
    "\n",
    "\n",
    "lats = [5,6,7,39,9,10,11,12,13,14,31.17]\n",
    "lons = [121,121,121,-88,121,121,121,121,121,121,121.5]\n",
    "c = [1,1.25,1.5,1.75,2,2.25,2.50,2.75,3,3,2.5]\n",
    "\n",
    "\n",
    "la = [39, 40, 41]\n",
    "lo = [-88, -87, -88]\n",
    "c = [2,3,3]\n",
    "x, y = map(lo, la)\n",
    "x1, y1 = map(filtered_data[:,1], filtered_data[:,0])\n",
    "\n",
    "\n",
    "\n",
    "map.drawcoastlines()\n",
    "map.drawstates()\n",
    "\n",
    "map.hexbin(array(x1), array(y1), C = filtered_data[:,2], gridsize = 150)\n",
    "\n",
    "map.colorbar(location='bottom')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy.random import uniform\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "\n",
    "# create north polar stereographic basemap\n",
    "m = Basemap(lon_0=270, boundinglat=20, projection='npstere',round=True)\n",
    "#m = Basemap(lon_0=-105,lat_0=40,projection='ortho')\n",
    "\n",
    "# number of points, bins to plot.\n",
    "npts = 10000\n",
    "bins = 40\n",
    "# generate random points on a sphere,\n",
    "# so that every small area on the sphere is expected\n",
    "# to have the same number of points.\n",
    "# http://mathworld.wolfram.com/SpherePointPicking.html\n",
    "u = uniform(0.,1.,size=npts)\n",
    "v = uniform(0.,1.,size=npts)\n",
    "lons = 360.*u\n",
    "lats = (180./np.pi)*np.arccos(2*v-1) - 90.\n",
    "# toss points outside of map region.\n",
    "lats = np.compress(lats > 20, lats)\n",
    "lons = np.compress(lats > 20, lons)\n",
    "# convert to map projection coordinates.\n",
    "x1, y1 = m(lons, lats)\n",
    "# remove points outside projection limb.\n",
    "x = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), x1)\n",
    "y = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), y1)\n",
    "# function to plot at those points.\n",
    "xscaled = 4.*(x-0.5*(m.xmax-m.xmin))/m.xmax\n",
    "yscaled = 4.*(y-0.5*(m.ymax-m.ymin))/m.ymax\n",
    "z = xscaled*np.exp(-xscaled**2-yscaled**2)\n",
    "\n",
    "# make plot using hexbin\n",
    "fig = plt.figure(figsize=(12,5))\n",
    "ax = fig.add_subplot(121)\n",
    "CS = m.hexbin(x,y,C=z,gridsize=bins,cmap=plt.cm.jet)\n",
    "# draw coastlines, lat/lon lines.\n",
    "m.drawcoastlines()\n",
    "m.drawparallels(np.arange(0,81,20))\n",
    "m.drawmeridians(np.arange(-180,181,60))\n",
    "m.colorbar() # draw colorbar\n",
    "plt.title('hexbin demo')\n",
    "\n",
    "# use histogram2d instead of hexbin.\n",
    "ax = fig.add_subplot(122)\n",
    "# remove points outside projection limb.\n",
    "bincount, xedges, yedges = np.histogram2d(x, y, bins=bins)\n",
    "mask = bincount == 0\n",
    "# reset zero values to one to avoid divide-by-zero\n",
    "bincount = np.where(bincount == 0, 1, bincount)\n",
    "H, xedges, yedges = np.histogram2d(x, y, bins=bins, weights=z)\n",
    "H = np.ma.masked_where(mask, H/bincount)\n",
    "# set color of masked values to axes background (hexbin does this by default)\n",
    "palette = plt.cm.jet\n",
    "palette.set_bad(ax.get_axis_bgcolor(), 1.0)\n",
    "CS = m.pcolormesh(xedges,yedges,H.T,shading='flat',cmap=palette)\n",
    "# draw coastlines, lat/lon lines.\n",
    "m.drawcoastlines()\n",
    "m.drawparallels(np.arange(0,81,20))\n",
    "m.drawmeridians(np.arange(-180,181,60))\n",
    "m.colorbar() # draw colorbar\n",
    "plt.title('histogram2d demo')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = plt.hexbin(np.array([1,2,3]), np.array([2,3,4]), C = np.array([1,2,1]), gridsize = 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sifutil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "'''Module to download CDL data'''\n",
    "from io import StringIO\n",
    "import re\n",
    "import requests\n",
    "from bs4 import BeautifulSoup\n",
    "import pandas as pd\n",
    "import sifutil\n",
    "import warnings\n",
    "BASE_CDL_URL = 'https://nassgeodata.gmu.edu/axis2/services/CDLService/GetCDLStat'\n",
    "CHAMPAIGN = 17019\n",
    "import numpy as np\n",
    "import math\n",
    "\n",
    "def from_csv(filepath):\n",
    "    df = pd.read_csv(filepath)\n",
    "    return df\n",
    "\n",
    "\n",
    "def get_by_fips(year, fips):\n",
    "    '''Grab CDL data for a county by FIPS code'''\n",
    "    url = BASE_CDL_URL + '?year=' + str(year) + '&fips=' + str(fips) + \"&format=csv\"\n",
    "    res = requests.get(url, verify = False)\n",
    "    returnurl = BeautifulSoup(res.text, 'lxml').find('returnurl').text\n",
    "    #print(returnurl)\n",
    "    rawdata = requests.get(returnurl, verify = False).text\n",
    "    raw_iter = StringIO(rawdata)\n",
    "    df = pd.read_csv(raw_iter, sep=\" *, * \", engine='python')\\\n",
    "           .apply(pd.to_numeric, errors='ignore')\\\n",
    "           .set_index(\"Category\")\n",
    "    return df\n",
    "\n",
    "\n",
    "def get_by_box(year, llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat):\n",
    "    '''Grab CDL data by bounding box'''\n",
    "    x1, y1 = sifutil.convertProjection(llcrnrlon, llcrnrlat, sifutil.WGS84, sifutil.CONUS_ALBERS)\n",
    "    x2, y2 = sifutil.convertProjection(urcrnrlon, urcrnrlat, sifutil.WGS84, sifutil.CONUS_ALBERS)\n",
    "    url = BASE_CDL_URL + '?year=' + str(year) + '&bbox=' + str(min(x1,x2)) + \",\" +\\\n",
    "              str(min(y1, y2)) + \",\" + str(max([x1, x2])) + \",\" + str(max([y1, y2])) + \"&format=csv\"\n",
    "    print(url)\n",
    "\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter(\"ignore\")\n",
    "        res = requests.get(url, verify = False)\n",
    "        returnurl = BeautifulSoup(res.text, 'lxml').find('returnurl').text\n",
    "    #print(returnurl)\n",
    "    with warnings.catch_warnings():\n",
    "        rawdata = requests.get(returnurl, verify = False).text\n",
    "    raw_iter = StringIO(rawdata)\n",
    "    df = pd.read_csv(raw_iter, sep=\" *, * \")\\\n",
    "           .apply(pd.to_numeric, errors='ignore')\\\n",
    "           .set_index(\"Category\")\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "# year = 2016\n",
    "\n",
    "# url = BASE_CDL_URL + '?year=' + str(year) + '&bbox=' + str(min(x1,x2)) + \",\" +\\\n",
    "#           str(min(y1, y2)) + \",\" + str(max(x1, x2)) + \",\" + str(max(y1, y2)) + \"&format=csv\"\n",
    "    \n",
    "# print(url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import lsq_linear\n",
    "fpars = []\n",
    "def get_fpars(query):\n",
    "    \n",
    "    from my_functions import coords_to_ind, get_cloud\n",
    "    base_lat = 38.4\n",
    "    base_lon = -89.59\n",
    "    \n",
    "    \n",
    "    \n",
    "    base_lat = 40.7\n",
    "    base_lon = -88.2\n",
    "    \n",
    "#     arg_max = None\n",
    "#     max_qc = None\n",
    "#     for i in range(105, 361,8):\n",
    "#         query = 'FPAR_A2016' + str(i) + '.hdf'\n",
    "#         raw_data = SD(query, SDC.READ)\n",
    "#         QC_reg_data = raw_data.select('FparExtra_QC')[:]\n",
    "#         FPAR_reg_data = raw_data.select('Fpar_500m')[:]\n",
    "#         row = 0\n",
    "#         left, right = base_lon, base_lon - 2\n",
    "#         bottom, up = base_lat, base_lat + 2\n",
    "#         l,b = coords_to_ind(left, bottom)[0], coords_to_ind(left, bottom)[1]\n",
    "#         r,u = coords_to_ind(right, up)[0], coords_to_ind(right, up)[1]\n",
    "\n",
    "#         Quality_flag = good_qc(QC_reg_data[u:b, r:l].flatten())\n",
    "#         if not arg_max:\n",
    "#             arg_max = query\n",
    "#             max_qc = Quality_flag\n",
    "#         elif Quality_flag > max_qc:\n",
    "#             arg_max = query\n",
    "#             max_qc = Quality_flag\n",
    "\n",
    "    #init\n",
    "    \n",
    "    \n",
    "    row = 0\n",
    "    rhs = np.zeros(16)\n",
    "#     mat = np.zeros((16,4))\n",
    "    raw_data = SD(query, SDC.READ)\n",
    "    QC_reg_data = raw_data.select('FparExtra_QC')[:]\n",
    "    FPAR_reg_data = raw_data.select('Fpar_500m')[:]\n",
    "    \n",
    "    for i in range(4):\n",
    "        cur_lon = base_lon\n",
    "        for j in range(4):\n",
    "\n",
    "            start_lon, start_lat = coords_to_ind(cur_lon, base_lat)[0], coords_to_ind(cur_lon, base_lat)[1]\n",
    "    \n",
    "            qc = QC_reg_data[start_lat, start_lon]\n",
    "        \n",
    "            \n",
    "            if get_cloud(qc)[0] == '0' and get_cloud(qc)[1] == '0':\n",
    "\n",
    "                end_lon, end_lat = coords_to_ind(cur_lon-0.01, base_lat+0.01)[0], coords_to_ind(cur_lon-0.01, base_lat+0.01)[1]\n",
    "\n",
    "#                 print(start_lat, end_lat)\n",
    "#                 print(start_lon, end_lon)\n",
    "                rhs[row] = FPAR_reg_data[(start_lat+end_lat)//2, (start_lon + end_lon)//2]\n",
    "                cur_lon -= 0.01\n",
    "            row += 1\n",
    "        base_lat += 0.01\n",
    "        \n",
    "        \n",
    "#     result = lsq_linear(mat, rhs, bounds = (0, 100))\n",
    "    return rhs\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 86.   0.   0.   0.   0.   0.   0.   0.  78.  78.  72.  72.  87.  89.  75.\n",
      "  80.  81.   0.   0.   0.   0.   0.   0.   0.  93.  76.   0.   0.   0.   0.\n",
      "   0.   0.  76.  89.  78.  88.  79.  83.  69.  73.   0.   0.   0.   0.   0.\n",
      "   0.   0.   0.  74.  74.  77.  75.  78.  78.  86.  77.  87.  87.  91.  75.\n",
      "  73.  77.  79.  72.]\n"
     ]
    }
   ],
   "source": [
    "#testing\n",
    "rhs = get_fpars('FPAR_A2016233.hdf')\n",
    "print(rhs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# fparreg_workflow()\n",
    "fparreg_workflow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fparreg_workflow():\n",
    "    big_mat = get_proportion_matrix()\n",
    "#     rhs = get_fpars('FPAR_A2016361.hdf')\n",
    "    print(big_mat)\n",
    "    save_matrix(big_mat)\n",
    "    mat = load_matrix('dajuzhen2.npy')\n",
    "    print(mat)\n",
    "    \n",
    "    \n",
    "def get_proportion_matrix():\n",
    "    from my_functions import get_fractions, get_by_box\n",
    "    mat2 = np.zeros((16,4))\n",
    "    base_lat, base_lon = 38.3, -89.59\n",
    "    \n",
    "    \n",
    "    base_lat = 40.7\n",
    "    base_lon = -88.2\n",
    "    row = 0 \n",
    "    \n",
    "    for i in range(4):\n",
    "        cur_lon = base_lon\n",
    "        for j in range(4):\n",
    "            print(base_lat, cur_lon)\n",
    "            mat2[row,:] = get_fractions(get_by_box(2016, cur_lon - 0.01, base_lat, cur_lon, base_lat + 0.01))\n",
    "            cur_lon -= 0.01\n",
    "            print(row)\n",
    "            row += 1\n",
    "        base_lat += 0.01\n",
    "    return mat2    \n",
    "    \n",
    "    \n",
    "def get_processed_matrix_and_rhs(mat, rhs):\n",
    "    indices = []\n",
    "    for i in range(len(rhs)):\n",
    "        if rhs[i] != 0:\n",
    "            indices.append(i)\n",
    "    indices = np.array(indices)\n",
    "#     print(indices)\n",
    "    return mat[indices, :], rhs[indices]\n",
    "\n",
    "\n",
    "def save_matrix(mat):\n",
    "    from tempfile import TemporaryFile\n",
    "    outfile = TemporaryFile()\n",
    "    np.save('dajuzhen.npy', mat)\n",
    "    \n",
    "    \n",
    "    \n",
    "def load_matrix(file):\n",
    "    return np.load(file)\n",
    "\n",
    "\n",
    "\n",
    "def run_regression():\n",
    "    from my_functions import get_cloud, coords_to_ind\n",
    "    from scipy.optimize import lsq_linear\n",
    "\n",
    "    time_series = np.zeros((4, 45))\n",
    "    ct = 0\n",
    "    x_values = []    \n",
    "    \n",
    "    prefix ='FPAR_A2016'\n",
    "    suffix = '.hdf'\n",
    "    ct = 0\n",
    "\n",
    "#     print(prefix+suffix)\n",
    "\n",
    "    for i in range(1,361,8):\n",
    "        a = str(int(i))\n",
    "        if i < 10:\n",
    "            a = '00'+ a\n",
    "        elif i < 100:\n",
    "            a = '0' + a\n",
    "\n",
    "        query = prefix + a + suffix\n",
    "#         print(query)\n",
    "        try:\n",
    "            data = SD(query, SDC.READ)\n",
    "            m2 = load_matrix('dajuzhen.npy')\n",
    "            rhs = get_fpars(query)\n",
    "#             print(rhs)\n",
    "            mat2, rhs2 = get_processed_matrix_and_rhs(m2,rhs)\n",
    "            \n",
    "#             print(mat2)\n",
    "            \n",
    "#             result = np.linalg.lstsq(mat2,rhs2)\n",
    "#             print(result[0])\n",
    "            \n",
    "            result = lsq_linear(mat2, rhs2, bounds = (0, 100))\n",
    "#             print(result.x)\n",
    "#             print('result', result[0])\n",
    "            ct += 1\n",
    "# #             print('result', result[0])\n",
    "            time_series[:,ct-1] = np.array(result.x)\n",
    "            x.append(i)\n",
    "                 \n",
    "        except Exception as e:\n",
    "            print(e)\n",
    "            continue\n",
    "                 \n",
    "    return x_values, time_series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "large_time = []\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "arrays used as indices must be of integer (or boolean) type\n",
      "arrays used as indices must be of integer (or boolean) type\n",
      "SD: no such file\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-179-749022387364>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[1;31m# print(rhs)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[1;31m# fparreg_workflow()\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtime_series\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrun_regression\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-178-919efdc9a651>\u001b[0m in \u001b[0;36mrun_regression\u001b[0;34m()\u001b[0m\n\u001b[1;32m     78\u001b[0m             \u001b[0mdata\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mSD\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mquery\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mSDC\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mREAD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m     79\u001b[0m             \u001b[0mm2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mload_matrix\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'dajuzhen.npy'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m---> 80\u001b[0;31m             \u001b[0mrhs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mget_fpars\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mquery\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     81\u001b[0m \u001b[1;31m#             print(rhs)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m     82\u001b[0m             \u001b[0mmat2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrhs2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mget_processed_matrix_and_rhs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mrhs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-174-38d1e5313b86>\u001b[0m in \u001b[0;36mget_fpars\u001b[0;34m(query)\u001b[0m\n\u001b[1;32m     39\u001b[0m     \u001b[0mraw_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mSD\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mquery\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mSDC\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mREAD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m     \u001b[0mQC_reg_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mraw_data\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mselect\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'FparExtra_QC'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m---> 41\u001b[0;31m     \u001b[0mFPAR_reg_data\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mraw_data\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mselect\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Fpar_500m'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     42\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mc:\\users\\bowen\\anaconda3\\lib\\site-packages\\pyhdf\\SD.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, elem)\u001b[0m\n\u001b[1;32m   1812\u001b[0m         \u001b[0mstart\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstride\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__buildStartCountStride\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0melem\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m   1813\u001b[0m         \u001b[1;31m# Get elements.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1814\u001b[0;31m         \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstart\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcount\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstride\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1815\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m   1816\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__setitem__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0melem\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "#testing1\n",
    "# m2 = load_matrix('dajuzhen.npy')\n",
    "# print(rhs)\n",
    "# fparreg_workflow()\n",
    "\n",
    "\n",
    "# July 16 days \n",
    "\n",
    "\n",
    "\n",
    "x, time_series = run_regression()\n",
    "\n",
    "\n",
    "large_time.append(time_series)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# # print(mat_3)\n",
    "# mat2, rhs2 = get_processed_matrix_and_rhs(m2,rhs)\n",
    "\n",
    "# result = np.linalg.lstsq(mat2, rhs2)\n",
    "# print(result[0])\n",
    "\n",
    "#test whole function\n",
    "# x, time_series = run_regression()\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 45)"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "large_time[0].shape\n",
    "# x, time_series = run_regression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x83baa28390>"
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzsnXd4m+XV/z+3ZFvy3ttJnL2c7YQV\nCCUNIWwoKZSy29LSlg46oO3bt0DppLx9344fdEHZu5RRZgYzIYvsxImd4W3L25a3ref3x61HlrWs\nZVsx9+e6csl6lm470tF5zvmec4SmaSgUCoVi4mIY7wUoFAqFYnRRhl6hUCgmOMrQKxQKxQRHGXqF\nQqGY4ChDr1AoFBMcZegVCoVigqMMvUKhUExwlKFXKBSKCY4y9AqFQjHBiRrvBQBkZGRohYWF470M\nhUKhOKXYtWtXo6ZpmSMdFxGGvrCwkJ07d473MhQKheKUQghR7s9xKnSjUCgUExxl6BUKhWKCowy9\nQqFQTHAiIkavUCg+PfT391NVVUVPT894L+WUwWw2U1BQQHR0dFDnK0OvUCjGlKqqKhITEyksLEQI\nMd7LiXg0TaOpqYmqqiqmTp0a1DVGDN0IIR4WQliEEAectqUJId4RQpTaH1Pt24UQ4g9CiDIhxD4h\nxNKgVqVQKCYsPT09pKenKyPvJ0II0tPTQ7oD8idG/0/gApdtdwEbNU2bCWy0PwdYB8y0/7sVeDDo\nlSkUigmLMvKBEerfa8TQjaZp7wshCl02Xwaca//5UeBd4E779sc0OZ/wYyFEihAiV9O02pBWqVB8\nytlb2crGw/Ve958zK5PiwrQxXJHiVCLYGH22brw1TasVQmTZt+cDlU7HVdm3uRl6IcStSK+fyZMn\nB7kMheLTwa/fKGHr8SY8OXaaBq/uq2XT91YpT1nhkXAnYz29yzxOH9c07a/AXwGKi4vVhHKFwgel\nFiufLy7gt1ctctv31LYKfvzSfo7UdzAnJ2kcVvfpYGBggKioU1O/Euyq6/WQjBAiF7DYt1cBk5yO\nKwBqQlmgQvFpp6Wzj0ZrLzOzEj3uP39+Nv/17/28vr9OGXo/eeyxx/jd736HEIKFCxdy3333ccst\nt9DQ0EBmZiaPPPIIkydP5qabbiItLY3du3ezdOlSEhMTqaio4Pjx41RUVPCd73yHb33rW+P964xI\nsIb+FeBG4Nf2x5edtn9TCPEMcBrQpuLzCkVolFqsAMzITvC4PyPBxGlT03ljfy13rJk1lksLmXte\nPcihmvawXnNeXhI/u2S+1/0HDx7kF7/4BR999BEZGRk0Nzdz4403csMNN3DjjTfy8MMP861vfYt/\n//vfABw9epQNGzZgNBq5++67KSkpYfPmzXR0dDB79mxuu+22oPXtY4U/8sqnga3AbCFElRDiS0gD\nv0YIUQqssT8HeB04DpQBfwO+PiqrVig+RZRaOgCYmeXZ0ANcuCCHUouV0vqOsVrWKcumTZu46qqr\nyMjIACAtLY2tW7dy7bXXAnD99dfz4YcfOo5fv349RqPR8fyiiy7CZDKRkZFBVlYW9fXek+SRgj+q\nmy942bXaw7Ea8I1QF6VQKIYorbcSF2MkLznW6zFr5+fw368c5I0DdczM9hziiUR8ed6jhaZpIyat\nnffHx8cP22cymRw/G41GBgYGwrvAUUD1ulEoIpwyi5UZWQkYDN6NU1aSmeIpqby+X0VKR2L16tU8\n99xzNDU1AdDc3MyZZ57JM888A8CTTz7JypUrx3OJYefUTCErFJ8iSi0dnDUjY8Tj1hXlcu9rhzje\nYGVapvcwz6ed+fPn85Of/IRVq1ZhNBpZsmQJf/jDH7jlllu4//77HcnYiYSQ0Zbxpbi4WFODRxQK\nd9q6+1l0z9vcecEcbjt3us9ja1q7OfPXm/jB2tl84zMzxmiFgXP48GHmzp073ss45fD0dxNC7NI0\nrXikc1XoRqGIYMrsihtfiVidvJRYlkxO4Y0DKnyjGI4y9ApFBFOmK268SCtdubAolwPV7VQ0dY3m\nshSnGMrQKxRjyCt7a7j2bx/jb8i0tN6KKcpAQWqcX8dfUJQDoLx6xTCUoVcoxpCNh+vZcqyJE42d\nfh1falfcGH0obpyZlBbHwoJk3jhQF8oyFRMMZegVijGktF7G3HeVt/h1fJnF6ld83pkLinLYU9lK\ndWt3wOtTTEyUoVcoxohBm8axBmnoP6kY2dBbeweobu0OuABqXVEuAG8qr15hRxl6hWKMqG7ppnfA\nhkH459Ef03vcBOjRT82IZ25uEm+o4qlRIyHh1KpTUIZeoRgj9J41q2ZlcrTeSlt3/wjH+y+tdOXC\nohx2lrdQ16YGcCuUoVcoxgxdE3/1ctnJe/cI4ZtSSwcxRgOT0/xT3DizboEM37x1UIVvXOns7OSi\niy5i0aJFFBUV8eyzz7Jx40aWLFnCggULuOWWW+jt7WXjxo1cccUVjvPeeecdrrzySsfz733veyxd\nupTVq1fT0NAAwLFjx7jgggtYtmwZZ599NiUlJQC8+uqrnHbaaSxZsoTPfvazjkZod999N7fccgvn\nnnsu06ZN4w9/+MOo/M6qBYJCMUaUWqxkJZo4e2YmBgGflLdw7uwsr8eX1VuZlhlPlDFwf2xGVgKz\nshN4fX8tN55ZGMKqR5k37oK6/eG9Zs4CWPdrr7vffPNN8vLy+M9//gNAW1sbRUVFbNy4kVmzZnHD\nDTfw4IMP8u1vf5tvfOMbw3rU33zzzYD8sli6dCkPPPAA9957L/fccw9/+tOfuPXWW3nooYeYOXMm\n27Zt4+tf/zqbNm1i5cqVfPzxxwgh+Pvf/85vf/tbHnjgAYAxaXusPHqFYozQm5PFm6KYm5vErhE9\nemvA8Xln1hXlsv1kMw0dvUFfYyKyYMECNmzYwJ133skHH3zAyZMnmTp1KrNmyV7+N954I++//z5C\nCK6//nqeeOIJWltb2bp1K+vWrQPAYDBw9dVXA3Ddddfx4YcfYrVa2bJlC+vXr2fx4sV89atfpbZW\n5kmqqqpYu3YtCxYs4P777+fgwYOO9YxF22Pl0SsUY4CmaZRZrHxuaT4Ay6ak8uKuKgYGbR499u6+\nQSpbuvjc0oKgX/PCBbn838ZS3jpYx3WnTwn6OqOKD897tJg1axa7du3i9ddf50c/+hHnn3++12Nv\nvvlmLrnkEsxmM+vXr/c6SlAIgc1mIyUlhT179rjtv/3227njjju49NJLeffdd7n77rsd+8ai7bHy\n6BWKcHPoZfj9Augf0rHXt/di7R1weOjLpqTS2TfIES+DQo41WNE0/1sfeGJWdgIFqbFsPdYU9DUm\nIjU1NcTFxXHdddfx/e9/ny1btnDy5EnKysoAePzxx1m1ahUAeXl55OXlcd9993HTTTc5rmGz2Xjh\nhRcAeOqpp1i5ciVJSUlMnTqV559/HpBf7nv37gVkeCg/X37JP/roo2P1qzpQHr1CEW4OvwptFdBy\nErJkt0FdcTPDPvd16eRUQMbp5+clu10ikGZm3hBCUJgerwqnXNi/fz8/+MEPMBgMREdH8+CDD9LW\n1sb69esZGBhg+fLlfO1rX3Mc/8UvfpGGhgbmzZvn2BYfH8/BgwdZtmwZycnJPPvss4DsZX/bbbdx\n33330d/fzzXXXMOiRYu4++67Wb9+Pfn5+Zx++umcOHFiTH9nZegVinBTsU0+tpQ7DH2Ziya+IDWW\nrEQTu8pbuP6MQrdLlFo6iDIIpqTHu+0LhNxkM++XNoR0jYnG2rVrWbt2rdv23bt3ezz+ww8/5Ctf\n+cqwbVar/P/8+c9/Pmz71KlTefPNN92ucdlll3HZZZe5bXcO4QAcOHDA59qDRYVuFIpw0l4jvXmA\n1nLH5lKLlZS4aDISYgDpbS+bkuo1IVtab6UwI56YqNA+orkpsVg6eukftLnvfOdnULYxpOtPdJYt\nW8a+ffu47rrrxnspIaE8eoUinFR8PPRzy0nHj2UWKzMyE4bNIl02JZU3DtRhae8hK8k87DJlFiuz\nc0Kf/ZqXbEbTwNLRS36K08zZzib46H+h8SjMcBv/rLCza9eu8V5CWFAevUIRTiq3Q3QcpM+QoRs7\nZRarW2J16RR7nN7Fq+8dGORkU2dI8XmdnGT5BVLrGqev2m5/3AnjMGUuEibbnUqE+vdShl6hCCeV\nH0P+Mkib7gjdNFl7ae7sY7rLHNf5eUnERBnc+t4cb+jEpsGMAJuZeSLP7sXXuLZCqLQb+k4LtFWF\n/DqBYDabaWpqUsbeTzRNo6mpCbPZPPLBXlChG4UiXPR1Qu0+WPkd6O2A8i1g188Dbl0oTVFGFuYn\nuxl6vcfNrBCklTq53jz6yu1gSoLedqjeCSmTQn4tfykoKKCqqsrRNkAxMmazmYKC4GsqlKFXKMJF\n9S7QBmHS6TL23dcB3S0Ow+2pynXZlFQe+egkPf2DmKONAJTVd2AQsgtlqCSao0k0RVHr7NEP9su1\nLr4Wdj8hwzfzr/B+kTATHR3N1KlTx+z1FCp0o1CED11WOWk5pBbKn1tOUmaxEh9jJC/Z/dZ76ZRU\n+gZtHKxpc2wrtVgpTI/HFGUMy7JyU8zUOHv0dfthoBsKV0LuQmn0FRMaZegVinBR+TFkzoXYVEi1\ntxxoLafMYmV61nDFjY5eOOUcvgm1x40rucmxwz16PT4/6TTIL4aaPTAY/rJ7ReSgDL1CEQ5sNqjc\nAZNPk89T7Ibe7tF7M9yZiSampMdJQ9/XRV97AycbO0NqfeBKXoqZ2jYnj75yGyQVQHK+TBwPdIPl\nUNheTxF5KEOvUISDhhLobZPxeQBzEsSm0td4grr2Hp8e+rLJqewqb0XbcDfaI+sYsGnMzApdcaOT\nmxxLo7WP3oFBuaFyO0xaIX8uWCYfq3eG7fUUkYcy9ApFOKi0F0rpBhQgZQo9DccBfBrupVNSabT2\n0le+g6j2SiDw8YG+0JU3dW090FYN7VUybAOQOhXi0qFKxeknMsrQKxThoGIbxGdC2rShbalTEK2y\nHYJPj35KKqBhaDqCcbAHk+h309yHgkNL39ozVCg1abl8FEKGb5RHP6EJydALIb4rhDgohDgghHha\nCGEWQkwVQmwTQpQKIZ4VQsSEa7EKRcRS+bH0kp0TrilTiO2qwRQFk1JjvZ46KzuRGaY2ogc6AZid\nohEbEx7FDThp6du6ZdgmKhZyFg4dkF8MDUegpz1sr6mILII29EKIfOBbQLGmaUWAEbgG+A3we03T\nZgItwJfCsVCFImLpqJd9bSafPnx7aiFRWj/L03p9jgM0GgTnZ7Y6nhelhbdiNDdZfsnUtvXIRGz+\nUjA6jaorWAZoUOO5e6Pi1CfU0E0UECuEiALigFrgPOAF+/5HgctDfA2FIrKp1PXzroZeKm+WJY3s\nKZ+WYHH8PDtlMGxLA4iNMZIaF42luRVq9w7PIwDkLZWPKnwzYQna0GuaVg38DqhAGvg2YBfQqmma\nLsqtAvJDXaRCEdFUbgOjCXIXDdvcmyDbCsyLbR7xErOM1Y6fpyWEX9OemxyLuWEf2AaGErE6cWmy\nN49KyE5YQgndpAKXAVOBPCAeWOfhUI/3oUKIW4UQO4UQO1XPC8UpTcXHMhwSNTwddaw/FZsmmGps\nHPESmd0nsGgpAEyO7Q/7EvNSzGS12meZFqxwP6CgWHr0qtHYhCSU0M1ngROapjVomtYP/As4E0ix\nh3IACoAaTydrmvZXTdOKNU0rzszMDGEZCsU40t9tD4ec5rartKmfelLJttX7voamEdV0lCMxclRd\ntqnH9/FBkJscy9SeQ9Jzj093PyC/GKz1Y97JUjE2hGLoK4DThRBxQtZ2rwYOAZuBq+zH3Ai8HNoS\nFYoIpvoTsPW7J2KRPeirtEwSu6s9nOhEew30ttOfJyWP5gHPA8NDITfZxELtCAP5yz0f4CicUuGb\niUgoMfptyKTrJ8B++7X+CtwJ3CGEKAPSgX+EYZ0KRWTiKJTy4NHXW2mOycOgjxb0RsNhAM5dtQbN\nGAM9bb6PD4KZUQ1kiHaa05Z4PiC7CIwxKiE7QQmpTbGmaT8Dfuay+TjgIQioUExAKrdDxiyZ0HSh\n1NJBX8IkaH0XBnohyuT5Gg1HADBkzwNzCnS3ej4uBKb1HJTLTVhAlqcDokxSW68SshMSVRmrUASL\nzSYVN65yRaBvwEZ5UxeG9EJA8x37thyWVbXx6WBOHhWPPrttL+1aLMc0H8MrCoqhVnWynIgoQ69Q\nBEtTKXS3uOvngfKmTgZsGok50+WGlhPer9NwBDLnyJ9jU6An/B59nGUXu7WZ1LT3ej8ovxj6u1Qn\nywmIMvQKRbBU2OPzXhKxAFmTZssNToPCh6FpsvOlbuhHw6PvacNgOUxJ1FxqW30oelRCdsKiDL1C\nESyV22Tnx/QZbrv08YGTC6eBIdoxKNyNjlo5tzXT/oUwGjH66l2ARnXCAmraur0flzoVYtNUQnYC\nogy9QhEsFR4amdkps1jJT4klzmySg7e9efQWqbgha658HA2PvnI7IGhPXzR80pQreidLlZCdcChD\nr1AEg7UBmo95lFWC9OgdU6JSC7179HbFjVvoJpwVqpXbIHs+qWnp1LZ2o/m6dkGxDCWpTpYTCmXo\nFYpgqNohHz0Y+kGbxvEGKzP0nvIpU7x79A2HIS4D4jPk89gU0AahzxqeddoGoWonTFpBXnIsnX2D\ndPT6UNXkF6M6WU48lKFXKIKh0e6JZ89321XV0kXvgM3Jo58C3c3Q66Hi1VlxA9Kjh/CFbxpKZA5g\n0mnkptj70vtKyOarTpYTEWXoFYpgaD4ute/mJLdduuJmhj4+0DEo3MWr1zSwlECWs6GXjc3ClpB1\ntFBe4RhA4jMhq3eyrP4kPK+viAiUoVcogqH5xPCxgU5sPmLBaBBD4wNTdUN/cviBHbVyoPhoevSV\n22VoKHXq0AASXx492BOyqpPlREIZekXkcfxd6elGMs0npBzRhTKLlae3V/KFFZNIjrVPcUoplI+u\nCdkG++/obOhj7R59OIqm+nvg2Gap8xeCrEQTBmEfKeiLgmKw1kH7CM3YFKcMytArIo8Xvwzv/nK8\nV+Gd/m5or/Lo0f/6jcPERhv5zmdnDW2MS4OYRPfQjcWDoQ+nR7/9r9Jgr7gVgCijgewksxwS7ov8\nYvlYpeL0EwVl6BWRRVczdDaMaV90n3JDT+gG28XQbznWyIbDFr7+melkJDg1MBNChm88efRx6ZDg\nNI9Bj9GHaui7muGD38GMNTBtlWNzbrJ5ZI8+R+9kqfT0EwVl6BWRReNR+ThGhr6urYflv9jAG/tr\n/T+p+bh8dDL0NpvGL/5zmPyUWG45yz2k41Fi6dz6QMdkT+6Gmoz94AGphV9zz7DNuSmxvoumQHay\nzC5SEssJhDL0ishCLyCy1svWvqNMmcVKo7WP7z2/lyN1fg780BuUpQ0Z9Jd2V3Owpp0fXjAbc7TR\n/Rzdo9fvHlx73OgYo2SYJxSPvqVchm0Wf9FN/pmXbKZmpKIpoC2tCK1mt+zQqTjlUYZeERJ1bT0c\nqA5jyb7u0YOcvDTKNFrll4mmwVcf30lbtx/zWpuPy1h6bCoA3X2D3P/WERYVJHPJwjzP56RMkZ0h\nO+3zYzvqpDF3NfQQegfLzb8AYYDP/NhtV25yLL0DNlq6vP+eTdZefrUnDtFn5eFX3qG5sy/4tSgi\nAmXoFSFx/1tHuOmRHYHHub3hbOgDDd9UbofewCpKdUP/p2uXUNXSzXee2Y3NNsLv0nxchm3sPW7+\n/sFx6tp7+MlF8zAY3PveAO4SS11xk+XB0IfS76Z2L+x7Fk6/DZLz3Xbn2Yumalq9x+k3HrbwyaC8\nWzmwYzNn/XoT9712CEt7+GfZKsYGZegVIXG80UqjtZemcHl9jUchd5H8ORB5X28HPLIONtwd2MtZ\n+4g2Cs6bk8XPLpnH5iMN/O/GUt8n6YYesHT08OB7x1g7P5sVU92nTDnQi6b0hKwnaaVOsB0sNQ3e\n/qnsQLnyux4PcWjpfcTp3z5UR3fSdLToOP57aS/rinJ4ZMtJVv52Mz/99wGqWroCX5tiXFGGXhES\nlc3yQ3+0PgwDrfu7ZXx52mfk87ZK/89tKQfbgPRm+zr9Pq3J2kt6vAkhBNedPoWrlhXwh42lvHOo\n3vMJg/3QWukw9L9/5yh9AzbuWjfX9wt58uhj02R1rSvBevTHNsKJ92DVD4dkmi442iB4Ud509g7w\nfmkjq+fnIXIXkdJygP+5ejGbv3cun1tawDM7Kjj3/nfZVd4S+PoU44Yy9Iqg6ewdoNEqPfnS+jA0\n4WoqAzTIXSirOdsC8Oh1T7m3HQ78y+/TGq29pCfEACCE4L7Li1iQn8x3n93DsQYPv1NrhWw6ljqV\nI3UdPLujkuvPmMLUjHjfLxQTL426vk6LPRHrocVxUIbeNgjv/EzeORTf4vWwjHgT0UbhVUv//tEG\n+gZsrJ2fA3lLoW4fDPYzOT2OX125gE3fO5cBm8b2E82BrU8xrihDrwiaiuahW/hSSxg8ej0+nzFb\nxpcDidHr0sWkAtj1iN+nNXX2DdO8m6ONPHT9MmKiDHz18V1YXTs9NuuKm2n84vXDJJii+Pbqmf69\nmC6x1BU3nuLzEFwydt9zUH8AVv+39yHkgMEgyPGhpX/7UD0pcdEsL0yVDc4GeoZ65gOT0uJIMEXR\n0DH6iihF+FCGXhE0uqFPMEVxNBwefcNRQED6dEieFJihby2HmAQ485uy0Kd2n1+nNVn7HB69Tn5K\nLH+6dgknGjv51tO7eWpbBX/cWMrdrxzkubffA+Dzz9fz/tEGvrV6JilxMZ4u7Y4usbTWS0PuKT4P\n0qPvs/o/pLu/GzbdB3lLYP6VIx6emxzrsd9N/6CNjYfrWT0nmyijQV4P3PT0mYkmLB0qMXsqoQy9\nImj0+Pw5szIore8IXXnTeFQaw+hYSC4ILBnbUi495kXXQJQZdv1zxFM0TaPB2ju8itXOmdMz+NG6\nOWwqsfDjl/bzwDtHefGTKgytJ+gWZlKz8rn9vBnccEah/2tMmSK/vOoPyOdeDX2A1bHb/iJbMqy5\nFwwjf6Tzks0eO1huO95Me88Aa+dnyw1p0+SXTs3wTpaZCSbl0Z9iRI33AhSnLuVNXSTHRrNsShqv\n769zC4METONRyLD3iEnKl/H2njavicVhtNoNfWwqzL9ChjLW3AumBK+nWHsH6BuwkZHg2SP/8tnT\nWDs/hyijIC0+BlOUEZ78B7TP4C83LA/890udIhPGZZvkc18ePUivPz595Ose+jdMOh2mnuPXMnJT\nYqnfX4vNpg2Tg759qA5ztIGzZ9oTxEJIr96lZXFmoonDdWoC1amE8ugVQVPR3MXktDhm2QdshKS8\nsQ1CY+mQoU8ukI/+JGQ1TXr0urJl2c3Q1wEHXvR5mp5ITo/3/uU0KS2O3ORYaeRBVsWmeWhx4A+p\nhfKx9C35hZSQ5fm4QDtYdtRDhvuAcm/kJZvpH9QcNQQgWzi8fbCec2ZmEhvjVNmbtxQsh2QnTDuZ\nicqjP9VQhl4RNJUOQy8HbISkvGmtgMFeyJwtnzsMvR9x+q4m6O8c0qpPWgFZ80ZMyjbZDV1Gop93\nIbZBKY/00od+RPT1NZV5V9xAYB0sbTbotEBCtt/L0LX0NU5a+v3VbdS190i1jTP5S+VdiB5uQhr6\njp4BevoH/X5NxfiiDL0iKAZtGpUtXUxOjyMr0USiOSo0j96huHHx6Nv9MPS6ZFH36IWAZTfJJGLN\nHu8vaTf06fF+JlPbq2GwL3hDn1wgWxOA97ANBGbou5ulIQ7E0DtGCg7F6d8+VIfRIAvHhqEnZJ3C\nN5n2L0bl1Z86KEOvCIq69h76BzUmp8UhhGBWdiKllhA8eldDn5ANhij/PHpdWql7zAALr4aoWJ9e\nvR66yfTXo/fQtTIgjNFS/gkjGPoAxgla7YVdARj6PA8e/VsH61lRmEaq65deUj7EZw1LyOp/L4sn\nQz/oR68gxZijDL0CkNK68ib/K0ormqTiZnJaHACzshNCU940HJFFUnH2NgIGIyTm+WfoXT16kHHu\noith/wueh3IjpZUAaf569LqG3sNkKb/R1+hNQw+BefRBGPqUuGjM0QaHR3+swUqZxTqktnFGCBm+\ncZJYZiZ48egth+FXBVCxze+1KMYGZegVADy65SRrfv8+bT66GjpT0Sy/FHRDPyMrkZaufoeXHDCN\npUPxeZ3kfP+SsS3lsp2AKXH49mU3Sz36/hc8v6S1l5S4aKKNfn4Mmo+D0SS93GDR7zp8efQx8fJu\nxp9krNUiH70ldj0ghCAveagvvd7uYY1rfF4nb6n8IrZ/YWbpoRuri6Hf/jdZYKX62EccytArAHi/\ntJG+ARslfsrmKpq7iDIIcpNlvFdX3gRVIatp0HgEMlwqTJML/Ot301o+3JvXKSiWAzS8hG+aOnv9\nj8+DNPSphX5p1b0yex3Mvsi3By6E/20QOurkYwAePcg4va6lf+tgHUX5SeSnxHo+OG8JoMnOmEB6\ngpw9O8yj77VKSSu4D0FXjDshGXohRIoQ4gUhRIkQ4rAQ4gwhRJoQ4h0hRKn9MTVci1WMDgODNnad\nlL1LSvwcvlHe1EV+aqysoITQlDedjdDdIlsfOJOUL3vSjzT8Qi+WckVPytbuddOCAzR29JEeiO6/\n+UTw8XmduRfDF57yrrjR8beDpdUC0fE+6wU8oVfHWtp72F3Rytp5Xrx5kKEbcPwNjQZBWryLxPLA\ni1LSGh3nPjJRMe6E6tH/H/CmpmlzgEXAYeAuYKOmaTOBjfbnigjmUG07nX1SKne41j+PXpdW6oSk\nvNETsZmzhm9PLgBbv5QPesNmk16/J48eYOHnpfHx4NU3dvY64s0joml2DX2Iht5f/PXorfWQGJg3\nD1JLb+no4c2D8o7gfG9hG4D4DEiePDxOn2iiwbkNwq5/QuZcmLrKfWSiYtwJ2tALIZKAc4B/AGia\n1qdpWitwGfCo/bBHgctDXaRidNE7Ec7ISuCwnx59hYuhD0l502gfH5jhwdCD7zh9R62UPHry6EEa\nzKIrYf+Lcoaq88t29Lr1ufGKtV5OiAq2WCpQAjH0AYZtQFbH2jR44uNypqQPFb15JX+Jm/LG4dHX\n7pX7im+Woa2Wk0MjEyOR9hpiF857AAAgAElEQVRoOjbeqxhTQvHopwENwCNCiN1CiL8LIeKBbE3T\nagHsjx6zREKIW4UQO4UQOxsaGkJYhiJUtp1oZmpGPGfPzOBoXQeDI0xYau/pp6Wrf5ihhxCUN42l\n0uvWpYc6DkPvI07vSXHjytKbZEHV0Tcdm/oGbLT3DPjfssEhrRwjQ+9vB0trfUCJWB09t3K03sra\n+TmIkUJJeUukAe+STsGwfjc7H5FS1oVXy/+H/k5ZxBaJDA7A41fAo5d+qubhhmLoo4ClwIOapi0B\nOgkgTKNp2l81TSvWNK04M9PD8AXFmGCzaew42cyKwjTm5iTR3T84rP2wJ3Rp5ZT04YY+aOVNwxFI\nn+Ge5PSnOtahoS/0fkz2PLfr6HNQ/fboQ9XQB8poe/TJQ4nX8+f5cX6ePU5v9+qzkkw0WHvRetph\n//Pyrik2ZejOKlLDN7sekS2i26ug8tMjAw3F0FcBVZqm6X+tF5CGv14IkQtgf/QRYFWEm9auPl7c\nVTWiV65z1NJBa1c/K6amMSdXJlRLRojT610rJ3nw6AFKA43Te5JWgkxIRsf77mLZWg4ISJnk/ZiY\neHkdJy9Tr4oNyKM3RMlY9VigJ2N93R3198gvg6BCN9Kjz0iIYclkP/QSeYvloz1On5lgon9Qo/uT\n56SEddnNcr/ez6f1ZMBrGnW6W2HzL6FghexwetD/ATWnOkEbek3T6oBKIYT+CV0NHAJeAW60b7sR\neDmkFSr85kB1Gxf/8UO+9/xeXt9f69c5enx+xdQ0ZmUnYhC4x+lbK6HkdcfT8ubhxVI6DuVNIHH6\nvk5oq3BX3IBUpowksWwph8Rcn8M2AJlQ7BwKEQ4Z+gCKpZIngXGMGr6ak2Uiut/7EG9HkjoIQ59k\njiYr0cS6olyM3gaau64nfSZU2w29XUtv+OSfkDVfSlkBUuxfhJEosXz/fqnuuugBmLkGDr0s+xd9\nCghVdXM78KQQYh+wGPgl8GtgjRCiFFhjf64YZV7YVcXnHtzCoE0jLT7Gb0O/7UQzeclmClJjMUcb\nmZoR7+7Rb/o5PHeDI6ZZ0dxFWnwMieboYYdlJZpIClR502gfxO2qodcZqWiqtXzIuPgiPtPF0MvQ\nTUAe/ViFbcC/DpYdgVfFOvPyN8/ixxeOMOvWmbyhhGxmookF4jjmxv0yCavH+E0JssI50kI3Tcdk\n3/4lX5SjKudfKcNe5VvGe2VjQkiGXtO0PfY4+0JN0y7XNK1F07QmTdNWa5o20/6ohkuOIr0Dg/zk\npf18//m9LJ2cyqu3r+TihblsPmKh03UMnguaJmd/rpia5kjGzclNGq6lt9mgbKP0Lvvk9srmLrew\nDUjlzczsxMC09Lqh9xS6AbtHP0KM3lciVsfF0OudK/3S0WtaeDT0geBPGwRH+4PAk7Eg4/TDWhKP\nRP5SqXJqryUz0cS1xo0MGGOlhNUZfZJWJPGOfcTieT+Vz2etlQKAT0n4RlXGnsLUtnVz9V8+5slt\nFXx11TQe/9IKMhJMXLQgl55+G5tKfKdHTjR20tDRy2nThoZbzM1JpKK5a2hWau1u6GqUP9sLeMqb\nutzCNjqzshM4aglAedN4RHZ09GZEkwpkiGLAQwOtgT7oqPEurXQmPkMWZukva+3FFGUg3h9D19UM\nvW2Ra+gTfWjgw4kjIbubrJg+LjVu4VjW+e6DYXSJZaRw4n0oeQ1WfnfobxUTD7MugEOv+D+y8RRG\nGfpTlC1ljVz8hw8pre/gwS8u5Ufr5jqqVIsL08hMNI0YvnGOz+vMyUkC4Iju1ZduGDqhu4WBQRvV\nrd1M8WLoZ2Yl0hqI8qbxqGwS5i3G7mhX7CF8014Fms1Pj94eo7d/ATVZ5TSsEWWFMPaKGwCzPUHq\nqzrWagGEDJWMBTkLQBih5hMSjvyLeNHLltRL3Y/TRyZGQvzbNghv/ljmV874xvB9RVdKJ+bk++Oz\ntjFEGfpTkP1VbVz/8HZS42N4+ZsrWbcgd9h+o0FwYVEOm0p8h2+2n2gmIyGGaRnxjm0O5Y3e86bs\nHdnIC6C7hdq2HgZtmlePfmagypuGo+6FUs4k2xuIeYrTe2pP7I34TNm33R7zbuzs8z8R22LvWjlW\nGnrw06Ovk19gY5UgjomDrLlQ/Qli1z85Kqayd9DDl58+MjGQmb+jxZ4noX4/rLlHziJ2ZsZn5UD5\nAxM/fKMM/SnI/20sJcEUxYu3ncmMLM8VjRcuyKV3wHf4ZptLfB4gPyWWRFMUJbUdMmRRtRNmnS93\n9rRS3uRZWqkTkPJmcACaj7m3PnAm2S6b9BSn96dYSifeXqthD980dngeCu6R5uNICacfrxMu/EnG\nWgObLBUW8pbIUEj9fjbHX0hDp4c7N11iOd7hm94O2PhzmHSaTL66Eh0Lsy+Ew6/KMOAERhn6SGWw\nH964U36onDhU086Gw/XcfFYhybHRXk4eCt/8Z59T+EbTZPk3UNXSRXVrNysK04adJ4RgTm6i9OiP\nbQI0WGBPtnW3OIqpXIuldAJS3rSWy/YFnqSVOkl58tGToW8pl9p2f9oGx9vDG3ZD39QZQPuD5uMy\nhBRt9u/4cOBvjH6sDX3+UpmYj47nUOZaz1OmIqVo6oP/kfmdtb/y3kSu6Er5ZXr83ZGv11Hn1kbj\nVEEZ+kjFchi2PQSPXS5lYfbY8p82S2/+5jN9hxH08M0w9c2eJ+H3RVB/yCk+n+527pycJEpqO9BK\n35Z93qefJ3d0t1LR3EWM0UB2kmejF5DypsFLjxtnomNlDNrTSMHWcmmADX4kVB0efQM2m+aI0ftF\n8/GxDduAnEYVHT+CoR8Pj96ekF3wOZKS0z0ben1k4ngqb1rKYeufZVuGgmXej5t+HpiS4eBLI1zv\nJPxpBbz4pbAuc6xQhj5Saa2QjzlF8MYP4ZVvUlbTyBsH6rjxzCkkx3n35nX08M1GPXyz82HQBmHn\nw2w/0UySOYrZOYlu583JTcTa24etdKP9g5Ag4/TdLVQ0d1KQGuuzyMZv5Y1jfKAXDb1Ocr53j97f\ncIqToW/v6WfApvnforj5eGhTpYLFnOw9GatpQfe5CYmchfCZn8CqO8lMNNHS1U/fgEvPGGO0NPbj\n5dEffQseu1R+2az+me9jo0yydXTJfzwru0CGdV64RSqvSt/xb+pZhKEMfaSiG/rr/gXn/BB2P4Hx\n8UuZFN3Ol1b6p/4oLkwjK9HE6/tq5R1C9S5pPPY+w77j1ayYmubRYM/JSWK+OImxu1FWEALEpjpC\nN5O9hG10/FbeNB6VHqkej/ZG8iTPyVhvA0c8EWe/c+lsDKwqtqdNtk4YS8WNjjnZe4y+u0WGvcba\nozcYYNUPIbnAUR3b6DppCuQX8FjH6JtPwFPXwFOfB0M0XPfCUDLfF/OvlEa8bKPn/RvvkZ+dNfcC\nGux9OqzLHguUoY9UWiukIiAuHc77CfUX/IXsrjJeM/2UtNb9fl3CaBCss4dv+nc+LuPZlz8IfR0s\nbN0wTFbpzOycRM41yGlCTF8tH2NToaeVCh8aeh2/lTeNIyhudDwVTfV1Srmkvx69MVr+Dp0NgVXF\n6nNix8PQx6Z4D93oIwSD6EUfLrzOjoWxLZrq74bNv4I/nyZzWmvuhdu2QOFK/86ftkqGKD0VTx19\nC7b+CZZ/Gc76NhSeDXueiuw2zB5Qhj5Saa2Qpf32JNL9lXO5ZvBe4mPN8PA62PusX5e5aGEegwN9\n2PY+Yx9jdyHtSbO4zrjBLRGrk2CKYq1pP+Wm2ZBgD3nEptBvbaa9Z2BEQ68rb3wmZDVNSiu9VcQ6\nk5Qvq3KdjZ5+x6MrPPwhTmrp9aHgfiVjx0NDr+PLow9iKHi40T16zwnZQrnGPt+dUENC02QPpj+f\nBu/9WoZgbt8pDXJUACMijdEw9xI48sbw3kJt1fDS1yB7AZz/C7lt8Rfle6Jiq//Xrz8oRy2OI8rQ\nRyq6oUe2HHhpdzVLV5yN8avvQsFyeOnWYY3GvFE8JZXL4w9i6m2CxdeBELyfdClFhpMUCS/DF7qa\nmW87ynvaoqFtsakMdMoE7kiGXlfe+JRYWuvl7bK/Hj0M9+oD0dDrxGe6hG788ejHuA+9M+YU6Pbm\n0UeQofcUunF0sawYvQXsfBie+YJM2N/4Klz18JBKK1DmXyG7cJa+LZ8PDsCLX5Zx+/WPDCmu5l0K\nMYmw+0n/rttwBB5aCTv+Fty6woQy9JGKk6H/f+8ewygEX1s1XcoEb/i3NJAb7h6xfNtgEHw54SMs\nWgrWyecC8HDHcnqEmSgvQ7M5vhkDNl62zqPbPmJQxuildzlSjN4xbcqX8saRiA3S0AeiodeJ1z36\nXgwCUuP88PpaTkhjGhM/8rHhxldP+hD73ISDjJFCNzC64ZuDL8nxhV/7EKaeE9q1Cs+Wd3x68dR7\nv4GKLXDx74eLBWLiYf7l8rX98dI3/0JWb+sN6MYJZegjke5W6e2mTKamtZsXdlWyvriAHPtUIIzR\nsjlT4xHY94zva3XUM7t9Ky8Ons3GI020dvXxSf0gZdnr5EDn7hb3c0o30BedzG7bDEot9vCLOYWo\nXmnoJ6X6NvQg4/Q+lTf+SCt1vHn0UbFDahp/sDc2a7D2kRYf41973rFuZuZMbAr0tnuehGStl7+/\nKWns12UnJspAalz0+Gjp+7uhcjvMWC0/D6FijIJ5l8mYfMnrsqXx4utg0dXuxy65Tk7ROjRCB/aa\nPUPHePqcjSHK0Eci+u1uymT+8t4xNA1uO3f68GPmXiI1zZt/JQdQeGPfswhtkI3mNby+v5YdJ+1v\nuOIvwUC3e6zfZoOyDfQXrsKGQVbIAsSmEm3rJifeQLxp5JL7EZU3jaUy2ezPrXZCtkwku3r0TjkM\nv4jPhO4WWjq6SI+P0PbEzpiTAU1+6btitUhvPpDffxTITDRh6fDw/kvIkl9Eo+XRV26DwV45jDxc\nFF0pPxPPXS8dkAt/6/m4SafJiWh7RgjfbLpPht9Sp/o3FnIUUYY+ErEb+uboHJ7eUcmVS/MpcPWi\nhYDP3i0LiXb+w/N1NA12PwGTTqNoYTGbjzSwqcRCTJSBGYvOhPxiu7beyeuu2wedFmLnXUBstJHD\nes8buwRybqp/jaocrRBcE7IDvfJO4ugb8pbYH0NlMEJi3vDeKYFIK3XiMwCN3o4GMhL9CNv0dcq2\nvOMRnwff1bHjURXrgWFDwp0RQv7/jJbE8sT7ssHalDPCd83JZ0BCjnQq1j/iPVwnBCy+Fso/Gsrh\nuFK+VfaJWvkdOf3MV3O6MUAZ+kjEbugfOWhjYNDG18+d4fm4aatg2mfg/d95Ls2u2inDO0uu46KF\nufQN2Hh+ZyWLJ6VgjjZC8S1yf/lHQ+eUvQOAYeYaZuckDvPoAWYl+tfSVZdYOpQ39Yfgjbvggdmy\n+EQDVt3p17UA96KplorAe884FU355dHrRmrcPHp7fYEnI9ExDsVSHshMMHlOxoJdSz9KHv2J9yF/\nGZjcC/6CxmCEK/8CX3gGsuf7PnbRF2RB1p6n3PdpmhzWk5ANK26V/4/Ko1e40VqBFh3P33e1ctni\nfAozfCQCV/83dDdLra8re56QwxXmX8GyyalkJ5kYsGmcpuvn518hvcadDw+dU7YRchdBQhZz7T1v\nNE2jP0Z6l4UJ/X79ClmJJrLNAyQffhr+thoePEPeeUw7F65/Cb69V8o9/cVZS9/dIsMZAXv00tCL\nzkb/FDe6kQpEwhlOTgGPPivJTENHr+dcTGqhvPMKt+a8px2qP5GOTriZdi5M/8zIxyXlyarxPU+7\nt2M+tkk6T2d/X94VxKaoGL3CA60VNMXk0t1vc4/Nu5K/FOZdDlv+BNahCUr0dcH+F+U+UyIGg2Bd\nkWxn7CiUiomDRdfK4QvWBuk5Vm6HGbIadk5OEi1d/Vg6erEMyBavk2O9eG8uCCH4nflhrqj6jQyB\nrP0V3FEC6/8pPyCGAN96SfmyIZttMDhpJTgamyUMtPqnode/WMZqILgr3jpYDvTJL/exGjjig8wE\nEz39tqFBNc6kTpHJ5HAbufItspVHqEqbUFn8RRk6PfHe0Dbdm0+eBMvso7P1Qe/jiDL0EYjWWs6h\nrhTOnpnhiHX75Lz/goEe+OB3Q9sOvyKLjJZc59h0y1lT+cKKySx3LpQqvkV2I9zzBBzfLD9AM3VD\nL1/7cG07VT3SA86J8TGs2oXpUQ18bJtH843vwRlfh3j3Bmp+k1wg12m1BCetBIdHny7a/Gt/0FYp\ne/zEj9FgD1e8efT6SMRICN3YtfQWn8qbk+F90RPvQ5QZClaE97qBMvtCacSdNfUl/4Ga3TIsqQ/T\niU2ViWNfg95HGWXoI5CBpnKO9adzy1l+JgEzZkqDvuMfQ97u7idktn/KmY7DJqfH8asrF8j4vE7m\nLKkh3vkIHH1bGpf8YmBo2lRJXQcnO+WbNjPK/zdrhminVkvlncNh0BA7T5oK1qM3p6AJI+mi3b/Q\nTVuVzA2Ml7JFj9G7GnprnXyMgNCNz+rY0dLSn3hPKl/Gsm20J6LNsGC9HFPY3SrvNjfdJxU5i74w\ndJx+ZzaO4Rtl6CMMrauF6IEOuuPyWTUrAI34qjtlMundX0nt98kP5MR7f4xU8c3yw7j/ORlWsU8s\nSo6LJi/ZTEltO8fa5XUSbX5OjgKie5rpN6Xx+v46/38Pbzi09JVyrebkkZuhuWIw0GdKI512/zpX\ntlUNve54EJMgE36ut/16n5sI8uh9a+lPhu8FOxuh/sD4h210lnxR3k0feFH+azgMn/nx8KlfvpLq\nY8QYzSBT+MvhkgPMA+bNm4/Bn4IeneR8meHf8kcZE0fI+Ls/zLlEVgV2NTri845duUmU1HXQnR5H\nB/Ek+qse6OtC9HeSVVDAR2WNtHX1+9Va2SsOQ18dWHtiF7qjU+0evZ8xen8Sc6OFwSALotw8er0q\ndvxj9Fm+DL05STYLC6fyRh/EE079fCjkLoas+fDJY/L/KXsBzLti+DF2xdp4Km+URx9hfLjjEwBO\nW7Ik8JNXflcahsOvSM/cnxatIBtALbtJ6odnrB62a05OImUWK8caOumOSvLfK+mSk5ymFxYyYNNC\nD9+YU+Qgjraq4DT0djqMqf6Fbgb7pYZ+PD16sHew9OLRB1IVPEokx0YTbRTeJZbh7mJ54n3ZayYv\niM/HaCCE9Opr98h2Gef9xF1oEDv+Hr0y9KONpsG/bpVzKUegprUbS2UpALGZQRTpxKXBWd+SPy/5\nYmDnrrpTtnZ1UXLMyU1iwKZRZrEyEJPsf5zRnjAsKJhMfkosb+yvHeGEERDCLrGstPcBCs7QtxiS\nyRTtw/MUnmivAbTxN/Se+t101ElPOZAOjaOEEILMBBOWdm+GvjD8Hn3hWWM3EN0fFl4tnaSC5TDr\nAvf95vGP0UfQX2uC0ngU9j0L/V2ybYEPHttaTj4WbDEJGPTbvUA581syGTT30sDOi4rx2DJ4Xu6Q\n6kcLpPDDPptVJGSxriiOx7aW097TT5I5lPBNvtRPD/QErW1v0pKYLvyY++mQVo63ofcgzYsQDb1O\nZuIIRVOHX5OJSn9GPvqirUoOk1/+5dCuE27iM+Da52RhnaecmD+D3keZU9vQH/y3rExLzHH6lys/\nBIm5MlkV6psrVI6+JR9r9/k8rKtvgKe3V/BYcgeG+CnBKz2iYmR3vTBRmB5PTJSBvgEbUfFp0HXU\nvxN1CWB8BusWJPP3D0+w6bCFy5f4GU7yRHKBfWA5QXv0lsEE4umWUrfoWO8HOgz9pKBeJ2yYk6UH\n74zVMq4DR1zJTDRR3eql31LqFCmLDUcYzBGfj5BErDMuIc9hmJIBoZKxQTPQI99ANbvthsWlAi9j\nNnxj2/g2ftL7W7eWy/9oL0qRf31STVt3PzNTmiFlnHqreCDKaGBWdgIHqtsxJaZDU2ChG+IzWZIc\nR06Smdf314Zm6JOcDEWQMfqa/gT7+hplDxJv6MPIk0JYbzjwNGXKWg+TTx+f9XggM9HMnkov7ZT1\nO6+W8vAY+rh0yJoX2nXGGoPBPv9XhW6CY9E18h/I5FlngzT8HXWyX/T+5+UfN87zJKVRp6dNTqLJ\nXgD1+6UszMN4M5tN45GPTrAgL4nYzhpIiRBFgZ05OUkcqG4nISVTfllp2shfnp2NsnthTDwG4IKi\nHJ7aXoG1d4AEP7pfesTZUKQEV61a2WtvJ9HZ4NvQt1VJoxIzckvmUcV1ytR4DQX3QWaiiebOXgZt\nmnvrZ2eJZeFZwb+IpsHx96Q3H2hVdSTgKak+hpyCf7Ehtp9o5mcvH5B9NozRsv9E/jKYc5Ec+AtD\nMz/DwB3P7uGXrx/2/4Tj74JtAM75nnzuJXzzQVkjxxo6+eppaYje9qCN2GhxwxlTuGvdHKLiU2Xl\nbK8fWvrOxmGqkAsXyKZqm0sswS9EVxHFZ/kOu3hh0KZxsiduaH2+aKsaf28epKEf6BlqRd3bLp9H\nWIzepkGTpzh98iRAhK68aToGHTWRGbbxB6fBPePBKW3oyyxWHt1azt4qD7eNemtZb21EA+TdIxb+\ntbuaRz464Vkz7Am90nTOJfKDWed5qPfDH54gM9HE2nx7w7AIM/QLC1LkdCs9QezPLWhnw7DWAcum\npJKZaOKNAyGob/R4eZBhm5auPho0+6COLj8M/XjH58G9OtZRLDX+GnodfUi4xzYIUTHyCzNU5Y3e\nTyZS9POBMs4dLEM29EIIoxBitxDiNfvzqUKIbUKIUiHEs0KIUdOAXbwoF3O0ged3VrrvdMQGQ/fo\nBwZt/OI/h8lOMtE/qPHMdj/mYNpssuXv9NVSCpazQPZ6d6HMYuW9ow3ccPoUojvsv0eEGXoHgRR+\ndDYM8+iNBsEF83PYXNJAV59/rY7d0IeUBJmIbbT20qQlD63PF+NdFavjZujHf4SgKz5nx4JdYnky\ntBc58Z7M0YxXy+hQGecOluHw6L8NOMczfgP8XtO0mUAL8KUwvIZHkszRXDA/h1f21tDT79IqNDpW\nDqsI1KOv2Q0nPhi26dmdlZRarNxzaRFnz8zgyW0V9A96GO/mTN1e+aGceb58nrMQGkrk4A0n/rnl\nBDFRBq49bfKwyVIRSSA9O1xCNwDrFuTQ3T/Ie0dGMLLeiI6VrZUDaW/sRJO1jy5MDBrNvg19T5sM\nkUSCoXeV5nVETp8bHZ/VsRB60ZTNJj+TU88Z94laQTPOHSxDMvRCiALgIuDv9ucCOA94wX7Io0D4\ntH4eWF88iY6eAd466KGfStrUwGP0b/0XPHOtw5h19PTzP28fZUVhGmvnZ3PjGYXUtffwzqERKj1L\n3wEEzPisfJ67UMbrLYcch7R19/PirmouW5Qne6+0Vsiqv2A19KONI3QzwhtW09xCNwArCtNIj4/h\n9QMh9L5Z/09YcFVQpzZaewGBLS7Dd4w+UjT04N7BMoL63Oj47HcD8g6so9b3yEtfWA7Ktsynanwe\n5GenpzX8vfn9JFSP/n+BHwK6e5sOtGqapt+bVwGjmtE6Y1o6+SmxvLCryn1n2tTAQzeNR6Q39/FD\nADz47jGaOvv4r4vnIoTgM3OyKEiN5dEtJ31f5+hbsld8gt2rzVkoH53i9C/vqaa7f5AbziiUG1or\nAp+DOpb4G6PvaZPaaRePPspo4Pz5OWw6XO9+BzYG6PNrRXyGb48+UjT04Dl0Y4yJKGfAHG0k0Rzl\nw6MvlI9tHkKs/nBcj8+fyoY+RTp6fZ3j8vJBG3ohxMWARdO0Xc6bPRzq8StMCHGrEGKnEGJnQ0OQ\nt/KAwSD43LICPixrpLrVpYVu6lT5wei1+nexrmZpAIwx8PGDVNfV8vcPT3DFknwWFsgPnNEguP70\nKWw70UxJnZcKy85GqN4FM9cOX0tMokN5o2kaT22roCg/iQUFdq+ttSLoROOY4G8pd1eTfPTQx/3C\nBTl09g3y/tHg/8+DpcnaS5RBYEzMGsHQ2w1SJHn0+t/capFhmwhzBrzOjoWh93SwcfoT78tqb397\nN0Ui49wGIRSP/izgUiHESeAZZMjmf4EUIYQulC4AajydrGnaXzVNK9Y0rTgzM7TmTOuXFciWMq5e\nvZ648fcN1lQmH8/5AfS2seeF3yKAH6wd3hrg88WTMEUZeGyrl7hj2QZAcwzwAKT2N6fIkZDdW9VG\nSV0H1yy3x+M1bcijj1SiY+UgjpGSsU5Vsa6cPi2dlLho3gglfBMkjdZe0hNiEPGZI4duDNGREQd3\nC93URVTYRiczwYehD6Vd8WC/HMt3KnvzMO4dLIM29Jqm/UjTtAJN0wqBa4BNmqZ9EdgM6EHUG4GX\nQ17lCExKi+P0aWm88EnV8NmVgUosG+3l/UWfo23yGs5qeJavn5FJXspwzXZqfAyXLc7jJXs1qxul\nb0utd+7i4dtzFkDdAbDZeGZ7BbHRRi5bbFeS9LTKkFEkG3oh7HrgEbwSp6pYV6KNBs6fl82GQ/X0\nDoxt+KbJ2ieHguuhG2/x0rZqqfCJhMKcaLOcpqQbCN2jjzCykszeVTcJ2fJ3CCYhW7UT+qynrqxS\nZ5w7WI7GO/lO4A4hRBkyZv+PUXgNN9Yvm0R5UxfbTzQPbUy1G3p/4/SNpWCMQUuZzC87LyFFdPK1\nuE0eD73hjEK6+wfdcwODA9Kjn7nG3VDkLIT+TjrrjvLK3houWZRLot7kK9IVNzr+yMR8GHqAdQty\n6egd4KOyEbTsYabR2ktGokmua7DPe+FXpGjodZw7WEZYVayOT4/eYJDv60A8+rr98PI34LHLZHvq\nwrMDWo+lo4eDNV7aMowHp3DoxoGmae9qmnax/efjmqat0DRthqZp6zVN87O6KDTWLcghwRTF886G\nNzZFtnN18ei3HW/ijxtL2VPZis3m5NU1lkLadN481Miz1RlUZ56Dafv/82gQivKTWTYllce3nhx+\njart8kOpyyqdyZUJ2d3b36erb5AvrHAy6qeMofejwk8Pi8R5nrV61vQMEs1RvLLHY1Rv1Gi09pER\nHzP0BeQtTq+PEIwUdEwDgxsAAB6HSURBVGne4ID820ZQsZROZqIJa++A9xqJlClei6bae/rZU9mK\nNjgg23k/chE8tBIO/Eu227713YDnDd/14n4+/9DW4Gs2ws04d7CMgHvT8BAXE8VFC3J5fX8tnc4T\n6dOmDZNYHq5t55Z/7uCBd45y+Z8/YvkvNnDHs3t4eU81gw1HGEyfzq/fLGFWdgLZl/xUfgPv+LvH\n17zhjCmcbOri/VIng1H6tuxN7WkyUeYcMERRW7KdOTmJLJ7k1ODsVDH0/uiBOxukF+qlX3pMlIGr\nlhXwyt4aSuv9H00YCpqm0dSpe/T2LyBPcXrboJxLGwmJWB3do9cb90WiRx+glr5/0MaGQ/Xc/uQO\nrr3vEV5/6C7afjMfnr1OfhbW/BzuOAQX/17ONQ6A6tZu3j1iobNvcGQZ9FjhrzR5lJgwhh5gfXEB\nXX2D/Md5yIWTxLKho5cvP7qTBHMUb37nbP7vmsWcMyuTzUcsfO+ZndiaTvDUMTPlTV385KJ5RE1e\nIStbt/zRo3JnXVEuGQmm4UnZo2/D5DOGkmjORJnoTplFVucRvrBiMsJZOdFSLqdDmVPcz4sk/IrR\nuxdLuXL7eTOJj4ni12+U+P3SvQODfPnRHby024OUdgQ6+wbp6beRPpJH31En+/lEkqHXO1g6qmIj\nL0Y/sqEvhJ42juzazAuPPMBz911H8jMX89vSi3gt6gf8OPppSnpSeaTg5/R/c5ccoBOkhPS5HZVo\nQFp8DC/trg7uFwo3MQkgjMqjDwfLpqQyLSOeF3Y6GYLUqdBWRU9PN197YhdNnb387YZi5uQkcdni\nfH5/9WJ2/tcaXrk2j2gxSKNpClcudRrMfe5dUi6482G319MrWjcfsVDe1Clv+S0Hh6ttXDikFTLf\nUM7li/KG74h0Db2OXvjhi84Gr2EbnbT4GG77zHQ2lljYeqzJr5f+8+ZjbDhs4W/vB97WQm+4lZFg\nGlqbJ0MfSRp6Hb2DpaNYKvIM/UjVsXVGGW6a/erlXFV+L+vZwMzMBGJW3AxX/BXt9k/Yvupx7imb\nzm1P7Q26zmLQpvHczkrOnpnJ1csn8UFpo71QbpwRYlzbIEwoQy+E1NRvP9nMyUZ7YULaNNBs/P65\nDewqb+GB9Ysdmngdo0EwL1p6S9/9wkX8z+ed1DKTVsC0c2HLH6Cvy+01v3jaZIxC8MTH5UO95531\n80509Q3wdnMWGaKNZFvz8J2RLq3UiU2RKohBD2ojnc5Gj9JKV245ayp5yWZ++frh4XkODxyubef/\nbS4jIyGGQ7XtHG/wszbCjv5hT0+I8R26iSQNvY4eutE9+ggaOqLjq99N/6CN27en8heuYsfCe7He\n8h4x/1VDyu2bMV74G1h0NSJ9Ot9aPZOfX17ExhILN/xjO+09Pt5jXnjvqIXath6uXTGJK5bkM2jT\neG3v2OaCvDKObRAmlKEH+NzSAgyCITWMXWJ5pGQvd6yZxUULcz2f2CRntZI+033fqruk97frEbdd\n2Ulm1hbl8OyOSgaPvA3Jkz2O5AN4bW8tu/vsnqJzy+JTQUOv40+s0aWhmTfM0Ua+d/5s9le38eo+\n7x/GgUEbP3xhHylx0Tz+pdMQAl7bF1gXTL0qNiPBBFEmOfXHl0cfCS2Kdcz20I3e5yY+8mL0qXEx\nGA3Co0f/581l7KjuZfLn7mP5ld8mYfJirzNfrz99Cn+4Zgm7K1u4+i8fY+kIrG3CU9sqyUgwsXpu\nNrOyE5mXm8RLY5z094o/d8OjxIQz9DnJZs6emcmLn1QxaNPYZJEThS4u6OX282Z4P7GxVI4fNCe5\n75tyhpR3ffR/cgSdCzeeUUhvTxe2Y5th1vlewy9P76igJ2O+fFK3d2hHdwv0dZxiht7LLahtUIa6\n/DD0AFcsyWdebhK/ffOI19v1v394gv3VbdxzaRFzc5NYPiWN13x8MXii0Tl0A0NaelfaqqQH7el9\nMF6Yk0GzyXmp5mSprY8wjAZBenyM25DwfVWt/HFTGVcsyWfdAi9OlguXLMrjHzcu52RjJ+sf2kpF\nk/udtCfq2nrYfMTC+uICoo3StF2xJJ+9la2caByf1gPDUKGb8LK+uIDath7+9sFxvvFyJd2YuXxK\nz/DkpyuNR2WZtTfOvUveOr/9U9newKkL5fLCVK7PrSba1sPTLXM8GqySunZ2V7Ry6YrZMm/g7NGf\nKoobGFkP3NUMaH4beoNB8OML51Ld2s1jW0+67T/eYOX37xxl7fxsLlwg47wXLczlaL2VowEodprs\nHn1avF0JFJ/p2dC3V0dWfB6GpHkNRyIyPq/jOiS8p3+Q7z67h8wEE3dfOj+ga50zK5Mnv3Iabd39\nXPXQFr88++d3VjJo07hm+dD/3yWL8hAC/h1iUlbTNG5/ejd//yCE+RYqdBNe1szLJiUuml+/UUJq\nXAxRmdOJ8lWVp2nS0Gf4kHEVroQ5F8OOv8HfzoNfFcBfz4XX7kDsfoIf5e2mX8Rwz4F0Lv/zR25G\n6JntlcQYDVy5tMBeIes0hORUMvQjlXL7aH/gjZUzM1g1K5M/bSqjtavPsd1m07jzxX2Yogz8/LIi\nxxf1ugU5GAQBxV6brL0kx0YTE2UYWl+XhyRwW2VkxedhSMHVWBrRhj7Lpd/Nb94s4VhDJ/evX0hy\nbHTA11s6OZWnv3I6bd393PXi/uFV7y7YbBrP7KjkrBnpTEmPd2zPSTZz5vR0/r2n2uf5I/HOoXpe\n3VvDq6HE+1XoJryYooxcXTyJ+Bgjf7uxmOiMab6rYzsbZAzUl6EHuPoJ+M5+WP8onH6blEztfx5e\n+SbGgy8QPX0VD960kkZrL5f88UMe3XISTdPo6R/kX59UcUFRjvQocxfK9ejVjqeUoR/Jo7cnOAMw\n9AA/unAO1t4B/ripzLHtiW3l7DjZwk8vnkdW0lC4IivRzGlT03ltX63fH95Ga59MxOp48+gjZeCI\nM/pdVF9HRBt658ZmW8oaeeSjk9x4xhTOnhl8L6u5uUncecEcNpVYeGaH9+6XH9ibGg4rQrRz+eJ8\nypu62F0ZnJEdGLTxmzelDLikroOBkWZReEOXydqCPD8EJqShB/jhBXPY8qPVzM9LlsqblpMyfuyJ\nRnsiNsNDItYZIaQxnn85rLkXbnoN7iyHb+6Ez/0DLvodn5mTxRvfPoczp6fzs1cOcss/d/D41nLa\newa4ZoX9ljJnkXysOyAfWytODQ09jByjH6H9gTfm5CRx1bICHtt6koqmLqpauvjNGyWcPTODq5a5\nG96LF+VyvLGTQ7VeOoi60GjtHYrP6+vrahr+nui1yt8rkhKxMLwmI8INfaO1l7aufr7//F6mZcRz\n17q5IV/3pjMLOWtGOj9/7ZCUMXvg6W0VpMXHsGae+9/ngqIcTFGGoMM3z++q4lhDJxfMz6F3wBZ8\nvN+cInMtvf69Z8PJhDX0RoMYul1MnSp7m7R7ue3Sm5mN5NF7wmCQXxALrnL03c5MNPHwTcu597L5\nbDnWxC9eP0xhehxnTLOXcecskI96+OZU0dCDU9tcb6Eb3aMP3Iu7Y81sjAbBb98q4Uf/2o8G/PKK\nBR5zK+uKcjEahN/qG2nonT36DPmhc/7CarcbgkiL0Q8z9JGnuNHJTDAxYNO447k91Hf08j9XLyY2\nxhjydQ0Gwf1XLcJoENzx3F4GXaS4lo4eNhyu56plBZii3F8v0RzNZ+dl89q+2pEnw7nQ3TfI7985\nytLJKXxnjXQE/XUu3BjHNggT1tAPY6Qulo2lEB0XVk9OCMENZxTy6u0rOWtGOt9dM2vIYCXmSEOo\nz5BtrQh6DuqYYzBKaaIvj14YgqpqzEk285Wzp/Havlo+KG3kzgvmMCktzuOxafExnDk9ndf21fgV\nvmnq7HPx6D0UTUWihh6GDATI906Ekpkow2sbSyx849zpw1t8hEheSiw/v6yIXeUtPPTesWH7XthV\nxYBLEtaVKxbn09zZxwelgc1BePijE1g6evnRhXOZnplATJSBQzXBGvrxa4PwKTH0el96L3H6xqOQ\nPn1U2tLOyk7kyS+fzmWLnb5EhJCdLGv3nVoaep1YHxPtOxsgLl1+IQTBV1dNJzvJxIqpaVx/uu8v\nv0sW5v3/9u49OO7qOuD490i70j6k3ZUsWd6VsWUbgy2IH9g8WhICGBISPDXTCSm0pbRlmnYKDWTS\nScljptN2mEkmLYR0OtNxQ1qmYUiYQgohJIbhMQlhAhbGYIyKbYz8kmxLliXZej9u//j9VlpJu6v1\nrnZ/+9s9nxmPtOuVf5eLdHT3/M49l2O9w7x3PH2XwvHJKfqGxq0WxXHJ2iAU0xGCiaoTSj2LeEW/\nNGTN7+XNIe67cYE0aBZ2bIpx64Yoj7x0gPdPWP/Pp6YMP37rGFevqmd1Y03Kr73ukkbqAl5++k7m\nN1N7B8f499c+4qb1TVzZUo+3soJLm2qzX9E72MGyPAJ9qNk6SCLlin6Bipt8WPYJ67Dw86fdU0Mf\nl67fTQZ9btKpqfaw64Hr+NE9V1NRkT6V9dnLluGtlAVr6nd3WLuQm0IZBHqpsPZTFJP4uygo6hz9\nZbEQOzbFePSOzTPVTYtIRHjotsupD1bxlZ/sZWR8kjc+OsPR3iH+8Or0Pz9Vngq2b4jx4v6TnMtw\nx+2/vnKQwbEJ/u6WmQ2QrdEQH3QOZFwEYIzhkZcO0N41oKmbvKuotPLnyQ4KHx+2VtSFDvTRDda5\nqgd3WY9dFejTbPzIsP1BOpFAVUaBIhzw8qm1jfz8va6ULRS6+of58pN7WdUQnL1hZzrQJ7RB6D8O\ntbGUuzYd5Sv+QB+o8vDoHZtZk2ZlnatIoIrv3r6Rg6fP891dH/Lk7qNEAl4+e9nCKa3bNjczOjHF\nrv0Ld7Q81jvEj357hNu3XMTaptrp51tjIc4MjnE6VfO2OY72DvHoywet5mqauimAVAeF9x4GzMIV\nN4stXnnT/rz10VWBPk1P+gwami2m7RuidPaP8M6x+b94RsYn+av/fpvhsQl23rVldi23v85avc8N\n9MWWtonzh6321/56p0fiuE9f0shd16zksdc/Ztf7J/n9zcvxeRdOFV6xIsKK+gDP7l24+uafX/yQ\nygrhKzfPXgC2xqw0WqZ5+rYO6/vyRN9wwkHvGujzJ96Xfu5brlwqbnIdjzcIh1+zHrsu0OcndXOh\nbm5tospTwc/enV19Y4zhmz99n3eP9/PIH2yatSoDrHd5gSXzUzfFGuh9EavHTTEcb1gEvv75daxq\nCDIxZbjzqsyqpESE2zbF+M2hHk4NpN5p+/6Jfp7d28mfX7uKZeHZ7SbWLbO+jzLN07cdsdKGXX3D\n9pnLVZqjz6u6VVbXxbkdC3sOAgL1awo7nvhh4ZOjVv7V74Ia+jifnbqZ+0tzYhRG+wsa6Gt9Xq6/\npJEX9nXNKrt7/I0Ont5znAduWstnUr2tT9w0NTVVfAeOJGq6HJqvcHoURSNQ5eE///RKvn/n5vm/\nxNPYsbmZKUPa082+88v/IxLw8pefnh8Tan1eVi4JZLyi322v6Dv7RqwiDIfaIBRhMjJPEkssaxIC\nUc8Bq266KnkZX14t2wDH3nTXah6sFb2ZtH5xVif8kA1mtys2V9s3xnjxg1Ps7ujlmtVLeOOjHv7p\n5+3c3NrEl9NVfwSWzIx5sNvaa1Gsgf5z33Z6BEWnpSFIS0Nw4RcmWNNYw8blYR56oZ0fvH6YFfUB\nVtQHWVEfYOWSAENjk/z6YA/funV9yrYNrdFQRiv6s4NjHDp9nppqD6fPjTA+OYXXoTYIZRToE0os\nV1w983zPgcLn5+PsM2TdF+gTysRmBfrsdsXmatu6pfi8FTz/XifNET/3PrGHVQ1BHv7ixvSVO8FG\n6LK7iBZraaVadP/yxY38Yt9JjvQOcfTMEL851MPTCamc5oifu34ndWlvazTEL/ef5PzoBDXVqUPo\n20es1fxnLmvimT0nODUwwnK/rujzK7ICkNkllsZAzyG44nedGVN8h6zrAn1CG4TEseewKzYXwWoP\n29Y18Yt9J9lzpI+JKcPOu7ZQ61ugkVawcWbMxbpZSi26i5fW8jfbZqd7RsYnOdY7xNHeIdY01iTd\nYRvXGgthDHx4coAtK1PfHN99pBdvpXDLZct4Zs8JOvtGWO6LwLkLO0thMZRPjt5TbaVoEkssBzph\nfNC5Ff3SVmi41Op37yapysSybGi2GLZviHJmcIz2kwN8/47NaTfPTAs2WvcUJkZ1RV/mfN5K1jbV\nsm1904LpoEwrb97uOMsnmsOsbrT+va7+Ycc6WJbPih6gvmV2iaVTFTdxnmq47y1nrp2LVDv8smhR\nvFhuWLeUDcvD7NjUzA3rMtw9mnikYP9xqxupGxrLKUctC/moC3jT5ulHxid573g/f3ZtC9GwH7BL\nLP0RGE6/kzsfyizQr4b2n808zrRrpZotVU/6wW6rfKy68Kcz+byVPHffJy/si+IppqGemT70bmgs\npxwlIrTGQmlX9PtO9DM2OcXWlnqC1R7Cfi9dfSMQjljvIqcms24Tko3ySd2AVWI5dGamD3zPASso\nFfFuw6KUqlVxvIbeLcEysQ1CMdfQq6LTGg2l7U0f3yi1ZaX1sxKL+OmMr+hhJgYVSHkF+njlTTxP\nH6+4cUtgKhapNn4MdjuStslaYuqmmGvoVdFpjYXS9qZv6+hlTWNw+ujKWNhnp24WOM8hT8os0Nu1\n9PE8/ZlDzuXn3UwkeRuEwe6CV9zkJB7o+45ZY9dArzLUGrV6DyXL009NGdqOnOXKlpmKnFjET1f/\niGNtEMor0NsHg9B7GEbPWas4zc9nJ1kbhMEz7gr01SHrnUnXXutxsR04oorW6sZgyt70H3Wfp394\nfDptAxCN+OgfHmfYY5d16oo+j6prrX4hvR9bq3nQFX22fHM6WBrjvtSNiPWLqTMe6HVFrzKTrjd9\nvO1B4oq+OWJV3pwetz4WetNUeQV6sLtYdsxU3CzRFX1W5tYDjw3CxHBBO1cuimADDNg19MV2Vqwq\naql607d19NJQU83KJTNtVWJ2oD8xajdJ09RNntWvtlI3PQdAKmfy9urCzN3K7VD7g5xNj1cgFHN0\nKMpdUvWmbztylq0r62addRy1u2AeG7LPLXZL6kZELhKRV0WkXUT2i8j99vP1IvKSiBy0P1744aH5\nVLfKys2f3Gfl7D3VC36JSmLuzViH2h/kLD7emib9XlAXJNkO2VMDIxztHWJry+yw1xTyUSFw4twk\nePyuSt1MAF81xqwHrgHuFZFW4EHgZWPMWuBl+3HxiJdYfvwrzc/nwl9nHYE4aR/L5uCu2JzEx6v5\neXWBkvWmb0uSnwcrp7+01seJvhFH2iBkHeiNMV3GmD325+eAdqAZ2AE8br/sceC2XAe5qOKpmvEh\nrbjJxXQbBPsb1u2pGw306gIl603fdqQXv7dyerWfKBbx2f1uCt/BclFy9CLSAmwG3gSajDFdYP0y\nAJI2HhGRL4lIm4i0dXd3J3tJfsRX9KAr+lzMbYPgYEOznAR0Ra+yN7c3fVvHWTZdFMFbOT+0Tu+O\ndeDwkZwDvYjUAE8DDxhjMjt2BTDG7DTGbDXGbG1sLOAq0F9nnegEuqLPxdwdfoM9VlMwr9+5MWVj\nekWvNfTqwrVGQ3ScGeT86ASDoxN80DXAlS3Jb0vGIn46+0cw/rB7UjcAIuLFCvJPGGOesZ8+JSJR\n+++jwOnchrjIRKwulqAr+lz453SwdFsNfVydfcBEw8XOjkO5UmJv+r3H+picMmxpSd6jPhb2MTYx\nxagn7J4VvVi1Q48B7caYhxP+6jngbvvzu4Fnsx9eniy52HrLHkh9aIBawNye9G5rfxDXeCn89W9h\nzTanR6JcKLHyZndHLxUCV6xI3uo6atfSn5Ngwcsrc2lTfC1wF7BPROythXwD+DbwlIjcAxwFbs9t\niHlw47fgyr9wehTulix147aTsuKWrnd6BMqlEnvTH+sdZt2yUMqTzeK7Y/umgjSOD1oVa5ULnIK2\nSLIO9MaY14FUbR+Le3lUv3r2TVl14Xz2fY6RhBV98xXOjUcpB8R70793vJ+OnkG+sCX1Tf34pqme\nyQBrwXo3XFOYd8HltzNWLY6KSuum9vBZmJqy+vy7MXWjVI5aoyH2dw4wODaZMj8PUB+sotpTwanx\nwrdB0ECvsue3A/1IH0xNaKBXZSmxZj5VxQ1Yq//miJ8TI3agL2CeXgO9yl68DUK8/YHbGpoptQji\nvembI/7p82FTiUZ8HJnud6MreuUG8Z70bm1/oNQiWN0YpNpTkXY1HxcL+zk8aN+ALWDqprwOB1eL\nyxeB/hPubX+g1CLwVlbwH3+ylVUNwQVfG434eeW8F6op6IpeA73KXrw5kwZ6VeauuySz7/3miI9+\nY/ep1xy9coXp1E08R7/E2fEoVeRiET8TeJj01mjVjXIJf8SqtjnbAf56qNQ3iEqlE79ZO+qp1Zux\nyiXiu2PPHNS0jVIZiEWs0sqhilpN3SiXiAf6noNacaNUBgJVHiIBLwMENXWjXCJ++MjogAZ6pTIU\nC/vpNUFN3SiX8CfUDWvqRqmMxCJ+eib8uqJXLuFPaMeqgV6pjMQiPrrGfJqjVy4xa0WvqRulMhGL\n+OkeD8DECIyPFOSaGuhV9rwBqLT7duiKXqmMRMM++rF30RYofaOBXmVPZOaGrDY0UyojzRE//cYO\n9AVK32igV7mJp290Ra9URqIR/8yKvkCVNxroVW6mA72u6JXKRFNtNQPUWA80daNcwR+BCs9MCkcp\nlZansoKqGvskKl3RK1eoWQq1MajQbyWlMhUI2++ANUevXOGGb8KdTzo9CqVcJVxnd3otUOpG2w2q\n3NQus/4opTIWjQTpN0FCw2eRAlxPV/RKKVVgsYiffhNg9NyZglxPA71SShVYzC6xHDvXW5DraaBX\nSqkCi4Z99JkaJof0ZqxSSpWkZntFLyP9BbmeBnqllCqwSMDLoNTgGdM6eqWUKkkiwpQvjG/iHBiT\n9+tpoFdKKSf46/EwAeNDeb+UBnqllHKAN2j3iSrA7ti8BHoRuUVEPhSRQyLyYD6uoZRSbuYLWbtj\nx87nv8Ry0QO9iFQC/wZ8DmgF7hSR1sW+jlJKuVlN2GrtffbM6bxfKx8r+quAQ8aYw8aYMeDHwI48\nXEcppVwrXG8F+r7e7rxfKx+Bvhk4lvD4uP2cUkopW33DUgDO97kz0Cfr0TOvfkhEviQibSLS1t2d\n//9QpZQqJkubYrwd+BTeuuV5v1Y+ulceBy5KeLwc6Jz7ImPMTmAnwNatW/NfSKqUUkXEXxNmy9ee\nL8i18rGi3w2sFZFVIlIF3AE8l4frKKWUysCir+iNMRMich+wC6gEfmiM2b/Y11FKKZWZvBw8Yox5\nAXghH/+2UkqpC6M7Y5VSqsRpoFdKqRKngV4ppUqcBnqllCpxGuiVUqrEiSlA0/sFByHSDRzJ8ssb\ngJ5FHE4p0DlJTudlPp2T+dw0JyuNMY0LvagoAn0uRKTNGLPV6XEUE52T5HRe5tM5ma8U50RTN0op\nVeI00CulVIkrhUC/0+kBFCGdk+R0XubTOZmv5ObE9Tl6pZRS6ZXCil4ppVQarg70egg5iMgPReS0\niLyf8Fy9iLwkIgftj3VOjrHQROQiEXlVRNpFZL+I3G8/X7bzIiI+EXlLRN615+Qf7OdXicib9pz8\nxG4tXlZEpFJE3hGR5+3HJTcnrg30egj5tP8Cbpnz3IPAy8aYtcDL9uNyMgF81RizHrgGuNf+3ijn\neRkFbjTGbAQ2AbeIyDXAd4BH7Dk5C9zj4Bidcj/QnvC45ObEtYEePYQcAGPMr4DeOU/vAB63P38c\nuK2gg3KYMabLGLPH/vwc1g9xM2U8L8Zy3n7otf8Y4Ebgf+zny2pOAERkOXAr8AP7sVCCc+LmQK+H\nkKfWZIzpAivoAUsdHo9jRKQF2Ay8SZnPi52i2AucBl4CPgL6jDET9kvK8Wfoe8DXgCn78RJKcE7c\nHOgzOoRclS8RqQGeBh4wxgw4PR6nGWMmjTGbsM5xvgpYn+xlhR2Vc0RkO3DaGPN24tNJXur6OcnL\nCVMFktEh5GXqlIhEjTFdIhLFWsGVFRHxYgX5J4wxz9hPl/28ABhj+kTkNaz7FxER8dgr2HL7GboW\n+D0R+TzgA0JYK/ySmxM3r+j1EPLUngPutj+/G3jWwbEUnJ1nfQxoN8Y8nPBXZTsvItIoIhH7cz9w\nE9a9i1eBL9gvK6s5McZ83Riz3BjTghU/XjHG/BElOCeu3jBl/yb+HjOHkD/k8JAKTkSeBK7H6rh3\nCvh74H+Bp4AVwFHgdmPM3Bu2JUtEPgn8GtjHTO71G1h5+rKcFxHZgHVjsRJrgfeUMeYfRWQ1ViFD\nPfAO8MfGmFHnRuoMEbke+FtjzPZSnBNXB3qllFILc3PqRimlVAY00CulVInTQK+UUiVOA71SSpU4\nDfRKKVXiNNArpVSJ00CvlFIlTgO9UkqVuP8H/wgA/1S/TWcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x83ba996470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "time_series = large_time[0]\n",
    "plt.plot(time_series[0], label = 'corn')\n",
    "plt.plot(time_series[1], label = 'soybean')\n",
    "# plt.plot(time_series[2], label = 'forest')\n",
    "# plt.plot(time_series[3], label = 'grass')\n",
    "plt.legend()\n",
    "# time_series\n",
    "s = \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_fractions(cdl):\n",
    "    \n",
    "        total_acre = sum(cdl['Acreage'])\n",
    "        if total_acre == 0:\n",
    "            corn = 0\n",
    "            soy = 0\n",
    "            forest = 0\n",
    "            grass = 0\n",
    "            return\n",
    "        if \"Corn\" in cdl.index:\n",
    "            corn = cdl['Acreage']['Corn'] / total_acre\n",
    "        else:\n",
    "            corn = 0\n",
    "        if \"Soybeans\" in cdl.index:\n",
    "            soy = cdl['Acreage']['Soybeans'] / total_acre\n",
    "        else:\n",
    "            soy = 0\n",
    "        pattern = re.compile(r' Forest')\n",
    "        trees = [cdl.index[i] for i in range(len(cdl.index))\\\n",
    "                 if re.search(pattern, cdl.index[i]) != None]\n",
    "        frst = 0\n",
    "        for tree in trees:\n",
    "            frst += cdl['Acreage'][tree]\n",
    "        forest = frst /  total_acre\n",
    "        grass = 1 - (forest + corn + soy)\n",
    "        return  np.array([corn, soy, forest, grass])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import lsq_linear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "result = np.linalg.lstsq(mat, rhs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(result[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_fpar_reg_series():\n",
    "    \n",
    "    from my_functions import get_cloud, coords_to_ind\n",
    "\n",
    "    time_series = np.zeros((4, 45))\n",
    "    qc = []\n",
    "    x_values = []\n",
    "    fpp = []\n",
    "    \n",
    "    \n",
    "    prefix ='FPAR_A2016'\n",
    "    suffix = '.hdf'\n",
    "    ct = 0\n",
    "\n",
    "    print(prefix+suffix)\n",
    "\n",
    "    for i in range(1,361,8):\n",
    "        a = str(int(i))\n",
    "        if i < 10:\n",
    "            a = '00'+ a\n",
    "        elif i < 100:\n",
    "            a = '0' + a\n",
    "\n",
    "        query = prefix + a + suffix\n",
    "        print(query)\n",
    "        try:\n",
    "            data = SD(query, SDC.READ)\n",
    "#             print('get')\n",
    "            FPAR_data = data.select('Fpar_500m')[:]\n",
    "            QC_data = data.select('FparExtra_QC')[:]\n",
    "        \n",
    "            \n",
    "            base_lat = 38.3\n",
    "            base_lon = -89.99\n",
    "\n",
    "            Q = QC_data[coords_to_ind(base_lon, base_lat)[1],coords_to_ind(base_lon, base_lat)[0]]\n",
    "            \n",
    "            if get_cloud(Q)[0] == '0' and get_cloud(Q)[1] == '0':\n",
    "#                 print('inside')\n",
    "                \n",
    "                \n",
    "                \n",
    "                result = get_sifs(query)\n",
    "                print(result.shape)\n",
    "                fpp.append(result)\n",
    "#                 print(result)\n",
    "#                 time_series[:,ct] = np.array(result)\n",
    "                ct += 1\n",
    "\n",
    "        except Exception as e:\n",
    "            print(e)\n",
    "            continue\n",
    "\n",
    "    return fpp\n",
    "\n",
    "# data = SD('FPAR_A2016361.hdf', SDC.READ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#testing\n",
    "time_s = get_fpar_reg_series()\n",
    "print(time_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mat = get_proportion_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "time_s = np.array(time_s)\n",
    "print(time_s.shape)\n",
    "\n",
    "for i in range(16):\n",
    "    if np.mean(time_s[:,i]) < 150:\n",
    "        plt.plot(time_s[:,i])\n",
    "\n",
    "\n",
    "        \n",
    "        \n",
    "1330Haiyan\n",
    "# plt.figure()\n",
    "# plt.plot(time_s[0], label = 'corn')\n",
    "# plt.plot(time_s[1], label = 'soybeans')\n",
    "# plt.plot(time_s[2], label = 'forest')\n",
    "# plt.plot(time_s[3], label = 'grass')\n",
    "# plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "base_lon = 20\n",
    "base_lat = 20\n",
    "for i in range(3):\n",
    "    cur_lon = base_lon\n",
    "    for j in range(3):\n",
    "#         rhs[row] = FPAR_data[coords_to_ind(cur_lon, base_lat)[1],coords_to_ind(cur_lon, base_lat)[0]]\n",
    "        cur_lon -= 0.25\n",
    "        row += 1\n",
    "        print(base_lat, cur_lon)\n",
    "    base_lat += 0.25\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "arg_max = 'FPAR_A2016249.hdf'\n",
    "\n",
    "print(get_sifs().x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(time_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from my_functions import get_fractions, get_by_box,hasdf\n",
    "\n",
    "\n",
    "base_lat, base_lon = 38.3, -89.99\n",
    "\n",
    "mat2 = np.zeros((16,4))\n",
    "row = 0\n",
    "for i in range(4):\n",
    "        cur_lon = base_lon\n",
    "        \n",
    "        for j in range(4):\n",
    "\n",
    "            print(base_lat, cur_lon)\n",
    "            mat2[row,:] = get_fractions(get_by_box(2016, cur_lon, base_lat, cur_lon - 0.1, base_lat + 0.1))\n",
    "            cur_lon -= 0.1\n",
    "            print(row)\n",
    "            row += 1\n",
    "        base_lat += 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from io import StringIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_by_box(year, llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat):\n",
    "\n",
    "    x1, y1 = sifutil.convertProjection(llcrnrlon, llcrnrlat, sifutil.WGS84, sifutil.CONUS_ALBERS)\n",
    "    x2, y2 = sifutil.convertProjection(urcrnrlon, urcrnrlat, sifutil.WGS84, sifutil.CONUS_ALBERS)\n",
    "    url = BASE_CDL_URL + '?year=' + str(year) + '&bbox=' + str(min(x1,x2)) + \",\" +\\\n",
    "          str(min(y1, y2)) + \",\" + str(max(x1, x2)) + \",\" + str(max(y1, y2)) + \"&format=csv\"\n",
    "    print(url)\n",
    "\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter(\"ignore\")\n",
    "        res = requests.get(url, verify = False)\n",
    "        returnurl = BeautifulSoup(res.text, 'lxml').find('returnurl').text\n",
    "    #print(returnurl)\n",
    "    with warnings.catch_warnings():\n",
    "        rawdata = requests.get(returnurl, verify = False).text\n",
    "    raw_iter = StringIO(rawdata)\n",
    "    df = pd.read_csv(raw_iter, sep=\" *, * \")\\\n",
    "       .apply(pd.to_numeric, errors='ignore')\\\n",
    "       .set_index(\"Category\")\n",
    "    return df \n",
    "\n",
    "\n",
    "def get_fractions(cdl):\n",
    "    total_acre = sum(cdl['Acreage'])\n",
    "    if total_acre == 0:\n",
    "        corn = 0\n",
    "        soy = 0\n",
    "        forest = 0\n",
    "        grass = 0\n",
    "        return\n",
    "    if \"Corn\" in cdl.index:\n",
    "        corn = cdl['Acreage']['Corn'] / total_acre\n",
    "    else:\n",
    "        corn = 0\n",
    "    if \"Soybeans\" in cdl.index:\n",
    "        soy = cdl['Acreage']['Soybeans'] / total_acre\n",
    "    else:\n",
    "        soy = 0\n",
    "    pattern = re.compile(r' Forest')\n",
    "    trees = [cdl.index[i] for i in range(len(cdl.index))\\\n",
    "             if re.search(pattern, cdl.index[i]) != None]\n",
    "    frst = 0\n",
    "    for tree in trees:\n",
    "        frst += cdl['Acreage'][tree]\n",
    "    forest = frst /  total_acre\n",
    "    grass = 1 - (forest + corn + soy)\n",
    "    return  np.array([corn, soy, forest, grass])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from my_functions import hasdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(mat2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import my_functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.array([[2,2,2,2,2],[1,1,1,1,1]]).T\n",
    "y = np.array([1,3,5,7,9])\n",
    "result=  np.linalg.lstsq(x,y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import p"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}