Clustering Coastal Edges.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# # Clustering- Coastal Edges\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\h5py\\__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import pandas_ml as pd_ml\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "\n",
    "import scipy\n",
    "\n",
    "%matplotlib notebook\n",
    "\n",
    "import matplotlib.pylab as py\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "#from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA as sklearnPCA\n",
    "from sklearn.metrics import f1_score\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.datasets import make_classification\n",
    "from sklearn.metrics import classification_report\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn import neighbors\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.model_selection import train_test_split \n",
    "from sklearn import preprocessing\n",
    "\n",
    "import itertools\n",
    "from matplotlib.colors import ListedColormap\n",
    "#import plotly.plotly as ply\n",
    "\n",
    "import tensorflow as tf\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>red</th>\n",
       "      <th>green</th>\n",
       "      <th>blue</th>\n",
       "      <th>nir</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>1.559241e+07</td>\n",
       "      <td>1.559241e+07</td>\n",
       "      <td>1.559241e+07</td>\n",
       "      <td>1.559241e+07</td>\n",
       "      <td>1.559241e+07</td>\n",
       "      <td>1.559241e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>3.869886e+05</td>\n",
       "      <td>4.795378e+06</td>\n",
       "      <td>8.719028e+01</td>\n",
       "      <td>1.142661e+02</td>\n",
       "      <td>1.349220e+02</td>\n",
       "      <td>2.134724e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>6.635469e+02</td>\n",
       "      <td>2.139469e+03</td>\n",
       "      <td>5.407919e+01</td>\n",
       "      <td>4.078841e+01</td>\n",
       "      <td>2.623506e+01</td>\n",
       "      <td>2.971749e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>3.856388e+05</td>\n",
       "      <td>4.791368e+06</td>\n",
       "      <td>2.000000e+00</td>\n",
       "      <td>4.500000e+01</td>\n",
       "      <td>7.200000e+01</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>3.864288e+05</td>\n",
       "      <td>4.793580e+06</td>\n",
       "      <td>5.200000e+01</td>\n",
       "      <td>8.900000e+01</td>\n",
       "      <td>1.210000e+02</td>\n",
       "      <td>3.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>3.870852e+05</td>\n",
       "      <td>4.795872e+06</td>\n",
       "      <td>6.300000e+01</td>\n",
       "      <td>1.030000e+02</td>\n",
       "      <td>1.280000e+02</td>\n",
       "      <td>7.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>3.874528e+05</td>\n",
       "      <td>4.796879e+06</td>\n",
       "      <td>1.070000e+02</td>\n",
       "      <td>1.280000e+02</td>\n",
       "      <td>1.390000e+02</td>\n",
       "      <td>2.900000e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>3.887862e+05</td>\n",
       "      <td>4.798792e+06</td>\n",
       "      <td>2.550000e+02</td>\n",
       "      <td>2.550000e+02</td>\n",
       "      <td>2.550000e+02</td>\n",
       "      <td>2.530000e+02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  x             y           red         green          blue  \\\n",
       "count  1.559241e+07  1.559241e+07  1.559241e+07  1.559241e+07  1.559241e+07   \n",
       "mean   3.869886e+05  4.795378e+06  8.719028e+01  1.142661e+02  1.349220e+02   \n",
       "std    6.635469e+02  2.139469e+03  5.407919e+01  4.078841e+01  2.623506e+01   \n",
       "min    3.856388e+05  4.791368e+06  2.000000e+00  4.500000e+01  7.200000e+01   \n",
       "25%    3.864288e+05  4.793580e+06  5.200000e+01  8.900000e+01  1.210000e+02   \n",
       "50%    3.870852e+05  4.795872e+06  6.300000e+01  1.030000e+02  1.280000e+02   \n",
       "75%    3.874528e+05  4.796879e+06  1.070000e+02  1.280000e+02  1.390000e+02   \n",
       "max    3.887862e+05  4.798792e+06  2.550000e+02  2.550000e+02  2.550000e+02   \n",
       "\n",
       "                nir  \n",
       "count  1.559241e+07  \n",
       "mean   2.134724e+01  \n",
       "std    2.971749e+01  \n",
       "min    1.000000e+00  \n",
       "25%    3.000000e+00  \n",
       "50%    7.000000e+00  \n",
       "75%    2.900000e+01  \n",
       "max    2.530000e+02  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filename='NOAA2017_DEM-0000014848-0000000000_valid_combined.xyz'\n",
    "\n",
    "d=pd.read_csv(filename, sep=' ')\n",
    "\n",
    "colNames=['x', 'y','red','green','blue','nir']\n",
    "\n",
    "d.columns=colNames \n",
    "\n",
    "d.to_hdf('Coastal.h5','table')\n",
    "\n",
    "\n",
    "\n",
    "#df = pd.read_hdf('GoldenGate19withRGB.h5','table')\n",
    "\n",
    "#df = pd.read_hdf('NJ_Sandy05.h5','table')\n",
    "df = pd.read_hdf('Coastal.h5','table')\n",
    "\n",
    "\n",
    "ndpoints=df.shape[0]\n",
    "\n",
    "df.describe()\n",
    "\n",
    "##print(df[1:1000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Original data is mislabeled. We will perfom K-mean clustering with 4 classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:2: FutureWarning: Method .as_matrix will be removed in a future version. Use .values instead.\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "adf=df[['x','y','red','green','blue', 'nir']]\n",
    "adfmatrix=adf.as_matrix()\n",
    "\n",
    "kmeans = KMeans(n_clusters=4, random_state=0).fit(adfmatrix)\n",
    "kmeans.fit(adfmatrix)\n",
    "y_kmeans = kmeans.predict(adfmatrix)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEICAYAAACavRnhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFWVJREFUeJzt3X+w3XWd3/HnywQUKwKSgEjQ6JrpLDgFMYNx7Oy4YjGgNcystqCF6LCT1sJUR7ddtDuLonbYdqoOW3WHlozBuiLr6ppFWDai1t0ZRC6oILAuVxYlghIIIK4Vhb77x/kEjpeT3HPuJ+Tkkudj5sz5ft/fz/f7+Xw44b7y/XFPUlVIktTjadMegCRp8TNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTaReSvC/J/572OKS9nWGifV6SNyeZSfKzJHcnuTLJP5/2uPaUJJXkxdMehxY3w0T7tCTvAj4K/BfgcOD5wMeBddMcl7TYGCbaZyU5CDgfOLuqPl9V/1hVv6qqv6yq/7iTff4syY+TPJjk60mOGdp2SpJbkjyU5EdJfq/VlyW5PMkDSbYn+ZskI//fS3JMki2t3U+SvLfVn57ko0nuaq+PJnl62/bWJH875ziPnW0k+WSSjyX5UhvbtUl+o237etvlO+3M7F93/UfVPssw0b7sFcAzgC9MsM+VwCrgMOAG4NND2y4G/m1VHQi8BPhKq78b2AosZ3D2817gCd9jlORA4MvAXwHPA14MXN02/2dgDXAccCxwAvAHE4z7dOD9wCHALPAhgKr6rbb92Kp6VlV9doJjSo8xTLQvOxS4t6oeGXeHqtpYVQ9V1cPA+4Bj2xkOwK+Ao5M8u6rur6obhupHAC9oZz5/U6O/FO/1wI+r6r9X1S9aP9e2bW8Bzq+qe6pqG4NgOGOCuX6+qr7Z5vppBqEk7TaGifZl9wHLkiwdp3GSJUkuSPL9JD8F7miblrX33wFOAX6Q5P8keUWr/zcGZwN/neT2JOfupIujgO/vZNvzgB8Mrf+g1cb146HlnwPPmmBfaV6GifZl1wC/AE4ds/2bGdyYfw1wELCy1QNQVddV1ToGl8D+Aris1R+qqndX1YuAfwm8K8mJI45/J/AbO+n7LuAFQ+vPbzWAfwSeuWNDkueOOR9ptzFMtM+qqgeBPwQ+luTUJM9Msl+Sk5P81xG7HAg8zOCM5pkMngADIMn+Sd6S5KCq+hXwU+DRtu31SV6cJEP1R0cc/3LguUne2W64H5jk5W3bZ4A/SLI8ybI27h2///Id4JgkxyV5BoPLb5P4CfCiCfeRfo1hon1aVX0YeBeDm9nbGJwdnMPgzGKuSxhcXvoRcAvwjTnbzwDuaJfA/h3wb1p9FYMb6z9jcDb08ar62oixPAT8CwZnLz8GbgN+u23+IDAD3AjcxODm/wfbfn/P4Km0L7d9fu3JrjG8D9jUnjb7VxPuKwEQ/3EsSVIvz0wkSd0ME0lSN8NEktTNMJEkdRvrl7WeCpYtW1YrV66c9jAkaVG5/vrr762q5fO122fCZOXKlczMzEx7GJK0qCT5wfytvMwlSdoNDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd32md+A77Hy3C9Nre87Lnjd1PqWpHF5ZiJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKnbWGGS5I4kNyX5dpKZVntOki1Jbmvvh7R6klyYZDbJjUmOHzrO+tb+tiTrh+ova8efbftmoX1Ikva8Sc5Mfruqjquq1W39XODqqloFXN3WAU4GVrXXBuATMAgG4Dzg5cAJwHk7wqG12TC039qF9CFJmo6ey1zrgE1teRNw6lD9khr4BnBwkiOA1wJbqmp7Vd0PbAHWtm3PrqprqqqAS+Yca5I+JElTMG6YFPDXSa5PsqHVDq+quwHa+2GtfiRw59C+W1ttV/WtI+oL6ePXJNmQZCbJzLZt28acqiRpUuP+41ivrKq7khwGbEnyd7tomxG1WkB9V8bap6ouAi4CWL169XzHlCQt0FhnJlV1V3u/B/gCg3seP9lxaam939OabwWOGtp9BXDXPPUVI+osoA9J0hTMGyZJ/kmSA3csAycB3wU2AzueyFoPfLEtbwbObE9crQEebJeorgJOSnJIu/F+EnBV2/ZQkjXtKa4z5xxrkj4kSVMwzmWuw4EvtKd1lwJ/WlV/leQ64LIkZwE/BN7U2l8BnALMAj8H3gZQVduTfAC4rrU7v6q2t+W3A58EDgCubC+ACybpQ5I0HfOGSVXdDhw7on4fcOKIegFn7+RYG4GNI+ozwEt2Rx+SpD3P34CXJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1GztMkixJ8q0kl7f1Fya5NsltST6bZP9Wf3pbn23bVw4d4z2t/r0krx2qr2212STnDtUn7kOStOdNcmbyDuDWofU/Aj5SVauA+4GzWv0s4P6qejHwkdaOJEcDpwHHAGuBj7eAWgJ8DDgZOBo4vbWduA9J0nSMFSZJVgCvA/5XWw/wauBzrckm4NS2vK6t07af2NqvAy6tqoer6h+AWeCE9pqtqtur6pfApcC6BfYhSZqCcc9MPgr8J+D/tfVDgQeq6pG2vhU4si0fCdwJ0LY/2No/Vp+zz87qC+nj1yTZkGQmycy2bdvGnKokaVLzhkmS1wP3VNX1w+URTWuebburPl//jxeqLqqq1VW1evny5SN2kSTtDkvHaPNK4A1JTgGeATybwZnKwUmWtjODFcBdrf1W4Chga5KlwEHA9qH6DsP7jKrfu4A+JElTMG+YVNV7gPcAJHkV8HtV9ZYkfwa8kcE9jvXAF9sum9v6NW37V6qqkmwG/jTJh4HnAauAbzI4y1iV5IXAjxjcpH9z2+erk/TR8x9Ce4+V535pKv3eccHrptKv9FQwzpnJzvw+cGmSDwLfAi5u9YuBTyWZZXC2cBpAVd2c5DLgFuAR4OyqehQgyTnAVcASYGNV3byQPiRJ0zFRmFTV14CvteXbGTyJNbfNL4A37WT/DwEfGlG/ArhiRH3iPiRJe56/AS9J6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6tbzFfSSpDFN69/pgT3zb/V4ZiJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG7zhkmSZyT5ZpLvJLk5yftb/YVJrk1yW5LPJtm/1Z/e1mfb9pVDx3pPq38vyWuH6mtbbTbJuUP1ifuQJO1545yZPAy8uqqOBY4D1iZZA/wR8JGqWgXcD5zV2p8F3F9VLwY+0tqR5GjgNOAYYC3w8SRLkiwBPgacDBwNnN7aMmkfkqTpmDdMauBnbXW/9irg1cDnWn0TcGpbXtfWadtPTJJWv7SqHq6qfwBmgRPaa7aqbq+qXwKXAuvaPpP2IUmagrHumbQziG8D9wBbgO8DD1TVI63JVuDItnwkcCdA2/4gcOhwfc4+O6sfuoA+5o57Q5KZJDPbtm0bZ6qSpAUYK0yq6tGqOg5YweBM4jdHNWvvo84QajfWd9XHrxeqLqqq1VW1evny5SN2kSTtDhM9zVVVDwBfA9YABydZ2jatAO5qy1uBowDa9oOA7cP1OfvsrH7vAvqQJE3BOE9zLU9ycFs+AHgNcCvwVeCNrdl64ItteXNbp23/SlVVq5/WnsR6IbAK+CZwHbCqPbm1P4Ob9JvbPpP2IUmagqXzN+EIYFN76uppwGVVdXmSW4BLk3wQ+BZwcWt/MfCpJLMMzhZOA6iqm5NcBtwCPAKcXVWPAiQ5B7gKWAJsrKqb27F+f5I+JEnTMW+YVNWNwEtH1G9ncP9kbv0XwJt2cqwPAR8aUb8CuGJ39CFJ2vP8DXhJUjfDRJLUzTCRJHUzTCRJ3cZ5mkvSU9TKc780lX7vuOB1U+lXTx7PTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHWbN0ySHJXkq0luTXJzkne0+nOSbElyW3s/pNWT5MIks0luTHL80LHWt/a3JVk/VH9ZkpvaPhcmyUL7kCTteeOcmTwCvLuqfhNYA5yd5GjgXODqqloFXN3WAU4GVrXXBuATMAgG4Dzg5cAJwHk7wqG12TC039pWn6gPSdJ0zBsmVXV3Vd3Qlh8CbgWOBNYBm1qzTcCpbXkdcEkNfAM4OMkRwGuBLVW1varuB7YAa9u2Z1fVNVVVwCVzjjVJH5KkKZjonkmSlcBLgWuBw6vqbhgEDnBYa3YkcOfQbltbbVf1rSPqLKCPuePdkGQmycy2bdsmmaokaQJjh0mSZwF/Dryzqn66q6YjarWA+i6HM84+VXVRVa2uqtXLly+f55CSpIUaK0yS7McgSD5dVZ9v5Z/suLTU3u9p9a3AUUO7rwDumqe+YkR9IX1IkqZgnKe5AlwM3FpVHx7atBnY8UTWeuCLQ/Uz2xNXa4AH2yWqq4CTkhzSbryfBFzVtj2UZE3r68w5x5qkD0nSFCwdo80rgTOAm5J8u9XeC1wAXJbkLOCHwJvatiuAU4BZ4OfA2wCqanuSDwDXtXbnV9X2tvx24JPAAcCV7cWkfUiSpmPeMKmqv2X0PQqAE0e0L+DsnRxrI7BxRH0GeMmI+n2T9iFJ2vP8DXhJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVI3w0SS1M0wkSR1M0wkSd0ME0lSN8NEktTNMJEkdTNMJEndDBNJUjfDRJLUzTCRJHUzTCRJ3QwTSVK3ecMkycYk9yT57lDtOUm2JLmtvR/S6klyYZLZJDcmOX5on/Wt/W1J1g/VX5bkprbPhUmy0D4kSdMxzpnJJ4G1c2rnAldX1Srg6rYOcDKwqr02AJ+AQTAA5wEvB04AztsRDq3NhqH91i6kD0nS9MwbJlX1dWD7nPI6YFNb3gScOlS/pAa+ARyc5AjgtcCWqtpeVfcDW4C1bduzq+qaqirgkjnHmqQPSdKULPSeyeFVdTdAez+s1Y8E7hxqt7XVdlXfOqK+kD4kSVOyu2/AZ0StFlBfSB9PbJhsSDKTZGbbtm3zHFaStFALDZOf7Li01N7vafWtwFFD7VYAd81TXzGivpA+nqCqLqqq1VW1evny5RNNUJI0voWGyWZgxxNZ64EvDtXPbE9crQEebJeorgJOSnJIu/F+EnBV2/ZQkjXtKa4z5xxrkj4kSVOydL4GST4DvApYlmQrg6eyLgAuS3IW8EPgTa35FcApwCzwc+BtAFW1PckHgOtau/OrasdN/bczeGLsAODK9mLSPiRJ0zNvmFTV6TvZdOKItgWcvZPjbAQ2jqjPAC8ZUb9v0j4kSdPhb8BLkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6GiSSpm2EiSepmmEiSuhkmkqRuhokkqZthIknqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG6LNkySrE3yvSSzSc6d9ngkaV+2KMMkyRLgY8DJwNHA6UmOnu6oJGnftSjDBDgBmK2q26vql8ClwLopj0mS9lmpqmmPYWJJ3gisrarfbetnAC+vqnPmtNsAbGir/xT43gK7XAbcu8B99zbOZe/zVJkHOJe9Vc9cXlBVy+drtHSBB5+2jKg9IRWr6iLgou7OkpmqWt17nL2Bc9n7PFXmAc5lb7Un5rJYL3NtBY4aWl8B3DWlsUjSPm+xhsl1wKokL0yyP3AasHnKY5KkfdaivMxVVY8kOQe4ClgCbKyqm5/ELrsvle1FnMve56kyD3Aue6snfS6L8ga8JGnvslgvc0mS9iKGiSSpm2EyZL6vaEny9CSfbduvTbJyz49yPGPM5a1JtiX5dnv97jTGOZ8kG5Pck+S7O9meJBe2ed6Y5Pg9PcZxjTGXVyV5cOgz+cM9PcZxJDkqyVeT3Jrk5iTvGNFmUXwuY85lsXwuz0jyzSTfaXN5/4g2T97PsKryNbhvtAT4PvAiYH/gO8DRc9r8e+BP2vJpwGenPe6OubwV+B/THusYc/kt4HjguzvZfgpwJYPfPVoDXDvtMXfM5VXA5dMe5xjzOAI4vi0fCPz9iD9fi+JzGXMui+VzCfCstrwfcC2wZk6bJ+1nmGcmjxvnK1rWAZva8ueAE5OM+gXKaXvKfN1MVX0d2L6LJuuAS2rgG8DBSY7YM6ObzBhzWRSq6u6quqEtPwTcChw5p9mi+FzGnMui0P5b/6yt7tdec5+wetJ+hhkmjzsSuHNofStP/EP1WJuqegR4EDh0j4xuMuPMBeB32iWIzyU5asT2xWDcuS4Wr2iXKa5Mcsy0BzOfdpnkpQz+Fjxs0X0uu5gLLJLPJcmSJN8G7gG2VNVOP5fd/TPMMHncOF/RMtbXuOwFxhnnXwIrq+qfAV/m8b+tLDaL5TMZxw0MvgfpWOCPgb+Y8nh2KcmzgD8H3llVP527ecQue+3nMs9cFs3nUlWPVtVxDL4V5IQkL5nT5En7XAyTx43zFS2PtUmyFDiIvfOyxbxzqar7qurhtvo/gZftobHtbk+Zr9apqp/uuExRVVcA+yVZNuVhjZRkPwY/fD9dVZ8f0WTRfC7zzWUxfS47VNUDwNeAtXM2PWk/wwyTx43zFS2bgfVt+Y3AV6rdydrLzDuXOdev38DgWvFitBk4sz09tAZ4sKrunvagFiLJc3dcv05yAoP/P++b7qieqI3xYuDWqvrwTpotis9lnLksos9leZKD2/IBwGuAv5vT7En7GbYov07lyVA7+YqWJOcDM1W1mcEfuk8lmWWQ5qdNb8Q7N+Zc/kOSNwCPMJjLW6c24F1I8hkGT9MsS7IVOI/BjUWq6k+AKxg8OTQL/Bx423RGOr8x5vJG4O1JHgH+L3DaXvqXlVcCZwA3tevzAO8Fng+L7nMZZy6L5XM5AtiUwT8e+DTgsqq6fE/9DPPrVCRJ3bzMJUnqZphIkroZJpKkboaJJKmbYSJJ6maYSJK6GSaSpG7/H5w50ttaAaGnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1=plt.figure(1)\n",
    "plt.hist(y_kmeans)\n",
    "plt.title('Class count')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15592413\n",
      "<class 'numpy.ndarray'>\n",
      "[0 1 2 3] [5423183 3935352 2735887 3497991]\n"
     ]
    }
   ],
   "source": [
    "print(y_kmeans.size)\n",
    "print(type(y_kmeans))\n",
    "unique,count=np.unique(y_kmeans,return_counts=True)\n",
    "print(unique,count)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndpoints=df.shape[0]\n",
    "nsubset=100000\n",
    "k=10  ## k nearest neighbors\n",
    "grid=np.floor(ndpoints/nsubset).astype(int)\n",
    "###\n",
    "edf=df.iloc[0:ndpoints:grid]\n",
    "cdf=df[['x','y']]\n",
    "ecdf=cdf.iloc[0:ndpoints:grid]\n",
    "sd=scipy.spatial.cKDTree(ecdf)\n",
    "nn=sd.query(ecdf.iloc[1:ndpoints:grid], k)\n",
    "adf=df[['x','y','red','green','blue','nir']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAADxCAYAAAA3MOvfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNX5+PHPMwlhh4QQFgkCKloQFWlELbUiVMTlB2LFolWQL0rrVrdWxS5UXEptK2BdqWjBDXFDimtEKViRTRZFUCioIPsOhiXJPL8/7hkYwiQzCTNzJ5nn7eu+knvuuec+40vn5NyziapijDHGxCLgdwDGGGOqD6s0jDHGxMwqDWOMMTGzSsMYY0zMrNIwxhgTM6s0jDHGxMwqDWOMqaFEpI6IzBGRRSKyRETuiZCntoi8JCIrRGS2iLStqEyrNIwxpubaB/RQ1VOAzkBvETmjTJ4hwDZVPQ4YBfylogKt0jDGmBpKPbvdaS13lJ3R3RcY735/BegpIlJemZlxj7KGadq0qbZt29bvMIwx1cD8+fM3q2rekZRx3jn1dcvW0tiet3jfEmBvWNJYVR0bnkdEMoD5wHHAo6o6u0wxrYDVAKpaIiI7gFxgc6RnWqURRdu2bZk3b57fYRhjqgER+eZIy9iytZQ57x4dU96Mlsv3qmpBRXlUtRToLCLZwOsi0klVPw/LEqlVUe76Ugl/PSUiGSKyQESmuvOZIrLQHWtFZLJLzxGR10Vkseu46RRWxs0i8rnryLklLL2JiBSKyHL3M8eli4g87Dp2FotIl7B7Brn8y0VkUKI/vzHGVIYCwRj/qVS5qtuB6UDvMpfWAK0BRCQTaAxsLa+cZPRp3AwsDZ2o6lmq2llVOwOzgNfcpbuBhap6MjAQGAPgKo9rga7AKcBFItLe3XMXME1V2wPT3DnA+UB7dwwFHndlNQGGA6e78oaHKhpjjEkFilKspTEd0YhInmthICJ1gZ8Cy8pkmwKE/oC+FPhAK1jJNqGVhojkAxcCT0W41hDoAUx2SR3xvvhR1WVAWxFpDnQAPlHVIlUtAf4D9HP3hHfgjAcuDkuf4DqBPgGyRaQlcB5QqKpbVXUbUMjhta4xxvgqji2NlsCHIrIYmIv3/TdVREaISB+XZxyQKyIrgNs4+Md3RInu0xgN3AE0jHCtH14rYac7XwRcAnwkIl2BNkA+8Dlwv4jkAnuAC4BQJ0NzVV0HoKrrRKSZSz/QseOscWnlpR9CRIbitVA4+ujY3i0aY0w8KEppnLasUNXFwKkR0v8Y9vteoH+sZSaspSEiFwEbVXV+OVkuB14MOx8J5IjIQuAmYAFQoqpL8cYNFwLv4FUuJdEeHyFNK0g/NEF1rKoWqGpBXt4RDYQwxphKC6IxHX5IZEujG9BHRC4A6gCNROQ5Vb3StRq6cvA1E67FMRi8jmxglTtQ1XF4TShE5AG8FgLABhFp6VoZLYGNLv1Ax46TD6x16d3LpE+P1wc2JpUENYggVDDk/ogt2vo1Dy55g+++34wilFLK/rAvswBQm0xKKKU4Qnp5+csvp4TisOdnAFlxTC+m5JC/SDs1bMWfOv+co+s3rfS/m6pSoNSnCiEWCas0VHUYMAxARLoDv1HVK93l/sBU1yzC5ckGilR1P3ANMCP06kpEmqnqRhE5Gu8V1pnutlAHzkj3842w9BtFZCJep/cOV7G8CzwQ1vndKxSjMdXZrv1FXDvrSVbu2VRunloIT55+LZ1y2h7x8zYV7eDCGRVOHAYgCOyJ8GIgXumlCU7/fNd3XDrzIf515nV0bNz6sOuJ4lcrIhZ+zdMYgPdFH64DMEFESoEv8Ka2h7zqWifFwA2uExtXxiQRGQJ8y8H3cm/h9X2sAIpwLRhV3Soi9+J1CAGMUNVyh5YZk6qKS0v4fPtqvtz+HQ8tfyu2e1D+b7Y37+u07GO5q9PFtG6QW6nnfrNzE/0/HlXpeKu72+c/y9s97k7KsxQoTuFtuMX2CK9YQUGB2uQ+kyreX/cZdy96MXrGSnqi4Bq6ND3mkLTS0lJ+N/9FPtj6RdyfV91kEuDj3vdFzSci86NNtovm5FNq6ZtvxfY67Oj89Uf8vMqyGeHGVAPri7bRZ8ZfE1b+r+YdNirehDm+QYvkPUyhNIX/lrdKw5gUN+Wbudy39HW/w0hrY8/8ZdKe5c0IT11WaRiTwko1aBWGjxoHalPYa3iSnyqURpwdkBqs0jAmhT259D2/Q0hLtclg0tm30bJu8lcZ8jrCrdIwxlTC17s3ctlHo/0OI62ckX0sd53Uj6PqN/E1Dm+ehlUaxpgY7C3dz08K/+R3GGkjE/hPrxHUCqTWV2EwhVsatnOfMSliwv/+YxVGDD7o+Ye4lfVx7wdSrsIItTRiOfxglYYxKWDcVx/wyPJ3/Q4j5T3V9Zc0qFWXOb0fILOGfn15y7EEYjr8UDP/rRtTjazctYEnV77vdxgp79Vut3FykzYHzid2u6WC3NFlpHC/QVAlpsMPqdUuMyYNDfjvGL9DSFmdGuQzrFNf2mcftoMBRzdsyrs97mbQfx9l/b4dlS67sOfv4xFi3CnCfs3wO4xyWaVhjI/WFW2LnimNPf3j6yu8npPVgCnn3HlImqpSuHYxH21cSlZGLVbs2sCq7zeQk9WAXi1Ppl/r02hZz98RUhXxJvel7ksgqzSM8dErX89KaPl/6nAJDy17k526L6HPSYS6UqtK94kIvVqdQq9Wp8Q5ouSxIbfGmIimrVucsLKzpS4n57bl/fMOndH8xfbVvPLtJ6zavYkd+4toWSeb/LpNeH1dai3MOeknt/kdgi9UhVK1loYxJoLNxd8nrOztuodLPnoIgF+0/jE3n3gBAB2zW/PH7MP3hrjquLO5ZObfj+iZtcjgwVN/wcqd63nu24/YWbyn0hsKXdSyM8M6XUKtjPT9egpaS8MYE1lyljN9fvVHbC3exT2df15unvz6uczp/QAAf/nsdV79bm65ecOdUL8lD/7wSlrWO7jkRrfmP+Cq9t0PnH+xfQ33ffYqK7/fcMhifA0DtWnXqDk/yj2eQceeTUYgdTuAk8XrCE/dr+bUjcyYNHB+i1N4Y/2CpDzr7fWLuIfyK41wd57UjztP8nZjLgmWsmDzKr7avZYu2e1oVi+bnKz6BAKxv0LpmJ3PC2fdXKW40411hBtjIpq7eUXSKoyQXfuKaFi7XqXuyQxkcFqz4zit2XEJisqUVZrOy4iISIaILBCRqe58pogsdMdaEZns0nNE5HURWSwic0SkU1gZt4rIEhH5XEReFJE6Lr2diMwWkeUi8pKIZLn02u58hbveNqysYS79SxE5L9Gf35hIPtv0NTfMezrpz90drH6jqNKNzQiHm4GloRNVPUtVO6tqZ2AW8Jq7dDewUFVPBgYCYwBEpBXwa6BAVTsBGXh7jAP8BRilqu2BbRzcV3wIsE1VjwNGuXyISEd374lAb+AxEbGXqCap9pbsY8j8sb48u0WdbF+eayonqIGYDj8k9Kkikg9cCBy2l6SINAR6AJNdUkdgGoCqLgPaikhzdy0TqCsimUA9YK2IiLv/FZdnPHCx+72vO8dd7+ny9wUmquo+VV0FrAC6xunjGhOTn7x/jy/PPbNJe7z/DUwq8xYsTN+WxmjgDiLvXtgPmKaqO935IuASABHpCrQB8lX1O+BvwLfAOmCHqr4H5ALbVbXE3b8GCK010ApYDeCu73D5D6RHuOcAERkqIvNEZN6mTZuq8rmNiWjIfx/z5bk/ym3PmK6DfXm2qRxFKNaMmA4/JKzSEJGLgI2qOr+cLJcDL4adjwRyRGQhcBOwACgRkRy8FkI74CigvohcCREHMofGL5Z3raJ7DiaojlXVAlUtyMvLKyd8YyrngsL7+GzXmqQ/94Uzf83o06zCqC5UoVQDMR1+SOToqW5AHxG5AKgDNBKR51T1ShHJxXst1C+U2bU4BgO4V0mr3HEesEpVN7lrrwE/Ap4HskUk07Um8oG1rrg1QGtgjXul1RjYGpYeEn6PMQmxcsd6Bsx6OCnPygAGtzuH7i1O5LhGLQhI6g7dNOWRuE3uE5HWwASgBd4bn7GqOqZMnu7AG3jftwCvqeqI8spMWKWhqsOAYWFB/UZVr3SX+wNTVXVvKL+IZANFqrofuAaYoao7ReRb4AwRqQfsAXoC81RVReRD4FJgIjAI74MDTHHns9z1D1z+KcALIvIQXqulPTAnUf8OTHqb9t1ihn02ManPnOUm55nqSyGerYgS4HZV/dT1I88XkUJV/aJMvpmqelEsBfo1T2MA3uuocB2ACSJSCnyBGwmlqrNF5BXgU7x/AQuA0NCTO4GJInKfSx/n0scBz4rICrwWxgBX1hIRmeTKLwFuUNXSxHxEk45UldFfTOXF1YldiDCSExq2TPozTWLEq5NbVdfh9QWjqrtEZCleP27ZSiNmSak0VHU6MD3svHuEPLPw/vKPdP9wYHiE9JVEGP3kWjD9yynrfuD+mAI3phKmr/2MOxa/GD1jgvyi3Y99e7aJHyUxGyy5+WqnArMjXD5TRBbhva7/jaouKa8cmxFuTBV9X7yXgR/9g9X7/N8T44ycY+l91Kl+h2HiQIHi2Neeaioi4csTj1XVwyYBiUgD4FXglrARqyGfAm1Udbfrg55MOX/Ag1UaxlTay6tm8dcv/+13GGQCxzU8ir90/gUt6+dEzW+qC6nMfhqbVbWgwtJEauFVGM+r6mtlr4dXIqr6log8JiJNVXVzpPKs0jAmRsXBEnq9N4LvKYmeOUFaZjXmjR53Rs9oqi2FuM32diNRxwFLVfWhcvK0ADa4wUJd8aZibCmvTKs0jKnA0u1ruHn202w/ONDPV+v276DrO3fzardbaN2wmd/hmASJ48593YCrgM/cHDjwlmw6GkBVn8AbYXqdiJTgjVAdoKrlrtlvlYYxYVSVq//7KEt3p/b0nZ/9dzSvn3U7rern+h2KiTNViVtLQ1U/IvKk5vA8jwCPxFqmVRrGAGu+38yVH/2DIi32O5SY9Zv59wObJpmaw+sIT911VK3SMGlp3779/K9oA+NXTufDTUuj32BM0tge4cb4Zvvm7dxz6d/5esm37N5SdCBdw46s9gH292gEPbKhYS2/QjUGCHWEp+5qxFZpmBppxisfc+/PR1W4Bbe4o87yIHWWb0ef3M7e1qDnNKbknEbQvE6ywjXmEH4tex4LqzRMjfP7vn9m9r8/rTBP2b/jQksg110NTNiBTthBENjbtTal1+ZBq7qJCfYIPHzqIL9DMAmQqBnh8WKVRgpSVdZ/vZHPZizli1lfsnzBSvbv38f+PaU0zKnP/xvai3N+cRZZWfYqpay1K9dHrTAiCf9fNFSBZAD15+xD56whKFB8bCb7B+ZCQaM4RVt1bWo34YzmJ/gdhkmQoLU0TDTBYJCHr/8nb459P2reL2c/zt+GPA7Ayed05G/v/8l2ZHOWfPTlEZcRsQJRyFhRQp0/bgA28P0fWlB6ZsMqld+7+Um8s+GzKsfXvWlHHiy4MnpGUy2pQnHQKg1TgRf/Opmn73y+Svcu/vALemVcxuvb/kWDxvXjHFn1k90svq2A8nbtqn/venY9mYW2rl3pMkecejkjuJyS0hIe/Hwyk9dFbxkJ8NdTr+InzTtU+nmmevFeT1mlYcrxwcSZVa4wwv08fyhv7jrycqq7jj8+PuHPELyKo/YLW9l7Z+zLkTer3Yip59x14DwzI5O7T7mUu0+5FIDZ67/iloXjKXW99/Uki35Hn8Yvjz+XOhlZ8fwIJsXFcUZ43Fml4bM/XxGfHd32f7+fb5aupk2H1tEz11A7tuzk0rwhSXxi5KFZWWRwVbuzyK3TkGMaNKNL7rExlXZ6i+OZ1dtW7U93NuTWRLRv7z6u++EdcS3z9h7DeWXd03Ets7pY/81Grmp3Q1KeFaoq9g08uIRHAJh69l00ret/J7mp7uz1VNoq3l/M2DufZfKYt5PyvB0bdvGfV2dx9s/OjCl/SUkJGlRqVeNRWB//ey7D+z6YtOeFKow997Sg+8mn8tuOfayiMHEXrz3CE8EqjQS55uRb+ebzNUl/7n39H2Laxafxm6euo1ETb3TP3MIF3H1e9DWKsurV4tG5f6FtNXjFVVJSwvm1L69w8l4iCHDrP4dywZBzk/tgkza80VNpvPaUiGQA84DvVPUiEZkJhMYqNgPmqOrFIpIDPA0cC+wF/k9VPxeRE4CXwoo8Bvijqo4WkSbuWlvga+AyVd3m1pAfA1wAFAFXq+qnLp5BwO9dWfep6vh4f+YL6l1O8V7/9lyYNXkuP5s8l1p1MisVx/6iYq498TYAzrmiG3c/d0uF+Yt2FfGHvg+yeHq5O0OSm5/Dkwv+SuPcxjHHUZ6SkhIeu+0Z3n7qA0p8+PfboEl9nl35GA0a1Uv6s036sMl9cDOwFGgEoKpnhS6IyKvAG+70bmChqvYTkR8AjwI9VfVLoLPLnwF8B7zu7rkLmKaqI0XkLnd+J3A+3naF7YHTgceB010lMxwowPsbdb6ITFHVuO3X+cYT7/haYYQ7kjg+fOG/TJ/4Me+VTIp4vW/OlRTt2Be1nC1rtnFp3jU0a5vH8ysfq3Qcqspz973M8/e/Sun+YKXvj4datTN5a49/e3+b9JPKr6cS2tsiIvnAhcBTEa41BHrg7UcL0BGYBqCqy4C2ItK8zG09gf+p6jfuvC8QaimMBy4OS5+gnk+AbBFpCZwHFKrqVldRFAK9j/yTHvTs8MhfstWRBpW+OVcdln5e1mUxVRjhNn69iet++NuY8hbt3sPPWgzh3EB/emVcxoThL/tWYQC8vOGw/3yNSZjQ6KlYDj8kuot+NHAHEOn/+H54rYTQ/rSLgEsA3JaDbYD8MvcMAML/5GuuqusA3M/QVmatgNVh+da4tPLSDyEiQ0VknojM27RpU7TPeIhgMMkv2ROsaMdermj3qwPnj9/2DMGSqn3GFQu+rvD6to07GHPjWPo2GsjOjTsrzJtMI6/6h98hmDQT1EBMhx8S9npKRC4CNqrqfBHpHiHL5RzaAhkJjHFbEn4GLICDmzGLSBbQBxgWy+MjpGkF6YcmqI4FxgIUFBRU6hvy6hE/5x83jKvMLSlv0zdbODfQPy5lnZvVH0ogo3aA0n3+tR4qY9ns5X6HYNKIqlCSwkNuExlZN6CPiHwNTAR6iMhzACKSC3QF3gxlVtWdqjpYVTsDA4E8YFVYeecDn6rqhrC0De61E+7nRpe+BggfApQPrK0gPW76XNeb47vGNpkrLbk/A6pLhQHQ9sSj/Q7BpJm0fD2lqsNUNV9V2+K9VvpAVUOrrPUHpqrq3lB+Ecl2rQmAa4AZYa+uwGuZlO2NnAKE1ocexMFO9SnAQPGcAexwr6/eBXqJSI4brdXLpcXVo5+MZPz/7JVGTfH7Sbf6HYJJI+nep1Gesn0TAB2AJSKyDK9VcXPogojUA84FXitzz0jgXBFZ7q6PdOlvASuBFcA/gesBVHUrcC8w1x0jXFrcHdWuBYXBl+l64amJKN4kyZiP76Nxrk3eM8mVypVGUib3qep0YHrYefcIeWbhDZGNdH8RkBshfQveiKqy6QpEXFNCVZ/Gmw+SENu37GBQ+xsp2r43emaT0m58ZAgdz7A9K0xy2TyNNKGqXJx7NUXbi6JnNinv0XkjOb6L9U0Zf6TyPA2rNOLkZ80HW4VRTQVqC+06tuH8a3rS66qzqdsg9bZ2NelDFUpsE6aabd/efeza/L3fYZhK6tzzRP781u/JrGX/G5jUYq+narhZU+b5HYKphBFT7uDMi07zOwxjIrI+jTRQu17lt/w0/rh13HVWYZiUp1Zp1Gw//OnJfodgYlAYfNnvEIyJSSp3hKdub0s1klUni2e+HON3GKYCQx68MnomY1KAavzmaYhIaxH5UESWisgSEbk5Qh4RkYdFZIWILBaRLhWVaZVGnOS3P4rC4Mtc8MvDpo0YnwUyhQG/6et3GMbESCgNBmI6YlAC3K6qHYAzgBtEpGOZPOFbSQzF20qiXFZpxNmtj/+KwuDLvLnneR54+27OG3y23yGltYzaGbyz76XoGY1JIaoS0xG9HF0X2oBOVXfh7W1UdmXv8raSiMgqjQTJqp3Faeedym/G3cgD7/7O73DSUp8bevHOnol4GzkaUz1Ucu2ppqFtHNwxtLxyRaQtcCowu8ylmLaMCLGO8CTYvm673yGkje4DfsStY39FPZugZ6or9fo1YrRZVQuiZRKRBsCrwC1lFoKFGLeMCLFKIwk6ndXB7xBqpKM7tuIfs0dSr34dv0MxJq7iOXpKRGrhVRjPq2rZRV+hkltGWKWRBC3bld211hyp51c/TrNWTf0Ow5i4U9cRHg/ivZsdByxV1YfKyTYFuFFEJgKnc3AriYisTyNJRky9y+8QahSrMExNphrbEYNuwFV4m+AtdMcFIvIrEQnt4xxxK4nyWEsjSc684Id+h2CMqSbiNSNcVT8icp9FeJ5yt5KIxFoaSZTX+rAtQY5I4+bpuTlQ7XpZ0TMZU015rYj4DLlNBKs0kmjUjBFxLW/MzPsoDL7MM1+NoWHTBnEtO5X1vam33yEYk1CpvHNfwisNEckQkQUiMtWdzwx7t7ZWRCa79BwRed1NY58jIp3CysgWkVdEZJmbDn+mS28iIoUistz9zHHp5U6LF5FBLv9yERlEEjVv04wXVz8Rt/IWfPgZAPnHHcVrG5+hMPgy9/77DjLr1uy3jr/43aV+h2BMQsWxTyPuktHSuBlvFiIAqnqWqnZW1c7ALA7u+303sFBVTwYGAuGLOY0B3lHVHwCnhJV3FzBNVdsD09w5lDMtXkSaAMPxRgh0BYaHKppkadoql8Lgy/zmmQr7mmLyxK0TDks748LTePv7FykMvsxDM/+EZNa8iW02B8PUZIoQDAZiOvyQ0KeKSD5wIfBUhGsNgR7AZJfUEe+LH1VdBrQVkeYi0gj4Cd6wMVR1v6qGZsv1Bca738cDF4elR5oWfx5QqKpbVXUbUAj48q7jvEHnHHEZ+/buq/D6Sd1OJCOzZr2BPKGrbcFqaj6N8fBDor9RRgN3AMEI1/rhtRJCsxMXAZcAiEhXoA3eJJNjgE3AM+4111MiUt/d0zw0ntj9bObSy5sWX6np8qnuqGMqnv8RDAYp2VuapGgSL+/oJjzyyUi/wzAmsdK1I1xELgI2qur8crJcDrwYdj4SyBGRhcBNwAK8FRozgS7A46p6KvA9B19Dlfv4CGlaQXrZ2IeG1nLZtGlTlEdV3T8//9sR3f/onIq/QCf8adIRlZ9KLruzDy98/aTfYRiTHCnc1EhkS6Mb0EdEvgYm4k0ueQ5ARHLx+hTeDGVW1Z2qOtj1dQwE8oBVeK2BNaoaWmTrFbxKBGBDaDVG93OjSy9vWnxM0+VVdayqFqhqQV5eXhU/fnRtO7bhXysrvw9H7fq1ea90Eg2yKx4xVfj89CpGllpGf3wf1/75Kr/DMCZp0rKloarDVDVfVdsCA4APVDW0E05/YKqq7g3ldyOkQgPwrwFmuIpkPbBaRE5w13oCX7jfpwChEVCDgDfC0ge6UVRncHBa/LtALzdSKwfo5dJ806rtUW4p9Re44nf9ImcKQOO8Rpw78Gye//Zxpu56LqaVWzev2RrnaJPrjD4/5N2SlzjxjBOiZzamhlAgGJSYDj/4NTZzAN7rqHAdgAkiUopXKQwJu3YT8LyrVFYCg136SGCSiAwBvsWrjMCbFn8B3rT4olB+Vd0qIvcCc12+EaqaEt+sWbVrMfjeKxh87xUx36OqiAilpaVkZGQcdj1Y7FdXWdU1bt6IRz55gBZtbL0uk6YUSPc9wlV1OjA97Lx7hDyz8IbIRrp/IXDY8r+qugWv5VE2vdxp8ar6NPB0TIGnEFXlk3c+5Y8XVtyP0fak1jwx/69kZB5eiaSiFsc244JretL3+t7Ua1jP73CMSQl+zcGIRc2eBVZDrP7qO/7vB7fElPfrz1bTO2sAl9x+IY3zGrJj064ER3e4+o3rIhkBGjZpQKMmDcnMyiTv6CYcd0o7TjnnRNp0aE1dW87cmPJZpWGqSlVjrjDCvfb3N+l1dXfe+9f0+AdVRq06JRScvRsC+2nWqphgSSannP97zu7fK+HPNqbm8a+TOxZWaaS46ZM+rvK9R1JhBAJBMrNKOOHU3Vx45WbadNhLbnOlTl3IzAQJQMV98Tey7vOzaNlpXJVjMCZtWUvDVNWy2cuP6P6MrAxK9x+c4Fe7bikdC3aS334vrY/ZS0GPnbQIG4QcqOR4uvIqDlVoljuT4L4lBGqfSDC4B3Y/CUUTgN1lcufi/V9Sy51nQb2hSINeIPU5OKiu7DOK0dJtIPUIZKTPgo2mhlNQn0ZGxcIqjRTXvku7Kt8bqBXg+W8eY0DLXwLQtcdWRjy7usJ7YhjJGxMR15m3rV/E5QAOteXwpKI/oEV/AGL7o8t7RjY0fY9AZnYlIjUmFaVupVGzFiaqgbr/vFuV7x23+O/kNm9CzlHZZGQq90xYjQgVHvEU7/Ki2w6buxJcf16yH2xMfKXpjHATB5m1Mpmyc0KlV6v967Q/kn9CPgAvrR5LTl6xD1/ifllFcP3xaCqPWzSmIlZpmCNRt0Fd3ts/iX999Q86nBFxKgsADXLq88C7v6cw+DKdzznpQLqIcErPnyQj1JSiG070OwRjKi80uS+WwwfWp1GNtDquBQ9//ECV7r1rwu0Ur3mSQIYfr438UkJw/1cEso73OxBjKiWVG8nW0kgjGc2nex3UfgeSTPtm+B2BMZUXlNgOH1hLo4YK7lsC2+4EvjrsWto0NABqd/U7AmMqTVL4LztradRAwfXdYVs/IlUYaWfrpQR3TEC14l0OjUkZsXaCW0e4iYfg+iFE2CIkve25D91wEsH1x3stMGNSWoyd4DVtPw2TfMH1XYCZfoeR2rb1I7jPWmAmxVlLwyRSsHgbwfU/5PDlOUxE2y4iuHOC31EYU75gjIcPrNKoxoLbJxBcfzxsOR1I/hLo1VrRfQTX/9jvKIw5XIrP07BKo5oKru8Pe+/zO4xqbiPB9ek36dGkPtHYjqjliDwtIhtF5PNyrncXkR0istAdf4xWplUa1VDlGSHpAAAbCklEQVQwGAQW+R1GDbGe4PouBEsqXsjRmKSKX5/Gv4DeUfLMVNXO7hgRrcColYaI3CgiOTGFF/n+DBFZICJT3fnMsFptrYhMduk5IvK6iCwWkTki0imsjK9F5DN3z7yw9CYiUigiy93PHJcuIvKwiKxw5XUJu2eQy79cRAZV9XP5qtg6cuNrN2zu6Y2uWn88wS2j/Q7ImLhQ1RnA1niWGUtLowUwV0QmiUhvkUovQnEzsDR0oqpnhWo1YBbwmrt0N7BQVU8GBgJjypRzjrsvfK/wu4BpqtoemObOAc7H22+8PTAUeBy8SgYYDpwOdAWGH0mF6JuMtn5HULMVP+ZVHjueRLXE72hMGqrE66mmIjIv7BhahcedKSKLRORtEYm6YFvUSkNVf4/35TsOuBpYLiIPiMix0e4VkXzgQuCpCNcaAj2AyS6pI94XP6q6DGgrIs2jPKIvMN79Ph64OCx9gno+AbJFpCVwHlCoqltVdRtQSPSmW8oJZNr+2kmx5+/oho4E99hSJCaJlMosI7JZVQvCjrGVfNqnQBtVPQX4Bwe/j8sVU5+GemtMr3dHCZADvCIiD0a5dTRwB5EHh/XDayXsdOeLgEsARKQr0AbID4UAvCci88vUpM1VdZ2LcR3QzKW3AsJfUq9xaeWlH0JEhoZq7k2bNkX5iD7JWxzHwmpBg5HQ7AsCLb6q8KDZUiDNNjnacQ3BYLHfUZh0kqR5Gqq6U1V3u9/fAmqJSNOK7om69pSI/BoYBGzGazH8VlWLRSQALMerFCLddxGwUVXni0j3CFku59AWyEhgjIgsBD4DFuBVUADdVHWtiDQDCkVkmXtXV27YEdK0gvRDE7zaeixAQUFBSq4CI4HMI/pvJtCiav0igUAGtJgDgGoQ3b8Ctv8Z9L9HEE01sPc/UO+nfkdh0kSy1p4SkRbABlVV98d6gIhbaR4Uy4KFTYFLVPWb8ERVDbqKoTzdgD4icgFQB2gkIs+p6pUikovXp9AvrLydwGD3QQRY5Q5Uda37uVFEXnf3zgA2iEhLVV3nXj9tdMWtAcJ2viYfb22NNUD3MunTY/h3kHJEjqTSiM86lSIBpPbx0PyZA2nBvcth1x+gdBWwA99mIMXbrret0jDJE6dKQ0RexPvOayoia/D6dGsBqOoTwKXAdSJSAuwBBmiU3cuifnuoarnjdlV1aQXXhgHDXODdgd+o6pXucn9gqqruDftw2UCRqu4HrgFmqOpOEakPBFR1l/u9FxAaFjYFrxU00v18Iyz9RhGZiNfpvcNVLO8CD4R1fvcKxVgt1R0Ce8ZV/r4mr0XPU0WBOu2hzsSI11RL0dItsHch7P49sD1hccSd/pvg+n9Do0kE6nX2OxpT08Wp0lDVy6NcfwR4pDJl+rU0+gC8L/pwHYAJIlIKfAEMcenNgdfdoK1M4AVVfcddGwlMEpEhwLd4lRHAW8AFwAqgCNeCUdWtInIvMNflG6GqcR2OlkyBxncSDOTA93+L/aZGkwhk/SBxQVVAJAPJbAYNenlHmOCeObBjEFDqS2wx23kZwZ3HEmjxtt+RmBoq1ol7fhHbR7liBQUFOm/evOgZfRYs/hq23gb6BYe/EsqDhn8jUP9MHyKrvOC+VbDtIiCFO5/rXUOgUcTuPJPGRGR+mWkBlVYnv7Xm33RbTHn/d9dtR/y8yrJNmGqIQK220Dxxr52SKVC7HbRYQjC4FzZfD8GP/A7pcEXPgFUaJkFSuaVhlYZJWYFAHWj2NADBknWwuQ9e53oqqCEd/CY1pXClYWtPmWohkNmSQIu5kPc50Clq/oQLJPWNgEknMc4G96s1YpWGqVYCGVkEWrwGzT4D8e+LW/LGR89kTFXZJkzGxFcgUJtA8xe8SYp1f5v052vR9KQ/06QPCcZ2+MEqDVPtBRpfC03eT+5Dd11PcP0VyX2mMSnAKg1TIwSyjnZLo7RP4lPnEdzypyQ+z6QNez1lTHIEWrwJjZ9L3gOLXyC49c/Je56p+awj3JjkCtTtijT/MnkP3P8MwfUnENzxePKeaWo2a2kYk1wiktyKA4U9owhuujJ6VmOisUrDmOQTEcj+ILkPLZ1DsPT75D7T1CiCjZ4yxjeBOvmQ8ZPkPnTTqdiabqbKrE/DGH8F8p4i2bsN6pZfJ/V5poax11PG+CvQYg5J/c+95F2C+1J0q2CT+qzSMMZ/gRbLgEbJe+C2bgR3vhE9nzFl2OspY1JEoMU8IIk77xX9luD6HxAssc5xUwnW0jAmdQRaTMLbtj5ZgrD5VILrjye469UkPtdUS2qjp4xJPc3mRs+TCN8PI7i+s42uMhVL55aGiGSIyAIRmerOZ4rIQnesFZHJLj1HRF4XkcUiMkdEOlVUjktrJyKzRWS5iLwkIlkuvbY7X+Gutw27Z5hL/1JEzkv05zepKRCo7ePTi9ANJxDcdKtVHiaidO/TuBlYGjpR1bNUtbOqdgZmAaE9Su8GFqrqycBAYExF5Th/AUapantgGzDEpQ8BtqnqccAolw8R6QgMAE4EegOPiUhGXD6lqX4alf1PLMlK3/Qqj/UnEgzu8zcWk1rStaUhIvnAhcBTEa41BHoAk11SR2AagKouA9qKSPPyyhERcfe/4pLGAxe73/u6c9z1ni5/X2Ciqu5T1VXACqBrXD6sqXYC9c4HsvwOAyiGjSdZq8N4Yq0wamKlAYwG7iDyhsr9gGmqutOdLwIuARCRrkAbIL+CcnKB7apa4s7XAK3c762A1QDu+g6X/0B6hHsOEJGhIjJPROZt2mRj7Wu0vNl+R3CA7pnpdwgmBQhp+npKRC4CNqrq/HKyXA68GHY+EsgRkYXATcACoKSCciRCmRrlWkX3HExQHauqBapakJeXV074piYIZNSH3Bl+h+HZeQ3Bfd/4HYVJAWlZaQDdgD4i8jUwEeghIs8BiEgu3muhN0OZVXWnqg52fR0DgTxgVQXlbAayRSTTFZEPrHW/rwFau2dlAo2BreHpEe4xaSpQqwXklve3TZJtO9fvCEwqSMfXU6o6TFXzVbUtXufzB6oaWje6PzBVVfeG8otIdmj0E3ANMMNVJBHLUe8F8IfApe6eQUBo+u0Ud467/oHLPwUY4EZXtcPb5m1O/D+9qW4CtRoCOX6HAUAwWOp3CMZv6VhpRDGAQ19NAXQAlojIMuB8vNFS0dwJ3CYiK/D6LMa59HFArku/DbgLQFWXAJOAL4B3gBtU1f4PNZ68j/2OwNnvdwDGTzG+morl9ZSIPC0iG0Xk83Kui4g87KYhLBaRLtHKzIyWIR5UdTowPey8e4Q8s4iywXOEclYSYfSTa8H0L6eM+4H7YwjbpJlARgZBGgK7/I0jUNfX55sUEL9WxL+AR4AJ5Vw/H+97tz1wOvC4+1mupFQaxlQb2X+H7UMT+oidezLZ9n0dFn7TnOXrc8lusJ/zT/4frXJ2QY4tM2Lit0SIqs4In9wcQV9ggnt9/4nrJmipquvKu8EqDWPCSO2zq/RH3u69tVi+IZvVW+oy48vjWLS6Jdv31EbJwHsLHGng3kFPTD/D/fY+bfLm8tTt/clt3LAKkZiaoBIjo5qKyLyw87GqOrYSjypvGoJVGsbEQkQOVBp7izP49JuWzFnRkpWbG/PFmmZs31eXg5VA6Ii59JhyfbNpB+fe5c1jvarnKdx6aY9KPMNUe5Xr5N6sqgVH8LSYpiGEs0rDmAi+Wp/DFU/+PIaclak0Ku/ZaYuYMusLPvjbDd6e5yY9JG9kVKWnIdgqt8ZEcMWTl7nfJMqReDuKirly5AtJeZbxX5JnhE8BBrpRVGcAOyrqzwBraRgTwdEks1KIxdJvN1IaDJIRsL/z0oEE41MjiMiLQHe8vo81wHCgFoCqPgG8BVyAtw5fETA4WplWaRhTVt4U4Am/ozjMjt17adKont9hmESL48Q9Vb08ynUFbqhMmVZpGFNGIKMeqdTKCGncIJm7DRo/+bWuVCys0jCmmrBXU2nEKg1jzJF47o4BfodgkiiVWxr2p4sxETSoU8vvEA546XdX0rFdS7/DMMlkCxYaU73cM7CX3yEA8OHffkX7fNvTJa2ot4xILIcf7PWUMRGcfUqFa2cmXP1aAT4cdROZGfZ3XboJzdNIVfZfpDERBAL+jp76vjhI1xvHcPm946NnNjWPamyHD6zSMCaFfbl2K+fc9ojfYZgkS9ftXo0xcbBjTzHfbtzqdxgmWWLtBLdKwxhTnknTF/kdgkmiVO4IT3ilISIZIrJARKa685kistAda0VkskvPEZHX3ZaDc0Skk0uv484XicgSEbknrOx2IjJbRJaLyEuhPcbdHuAvuS0MZ4dvQiIiw1z6lyJyXqI/v6m+6mVl+B3CAY3q1fY7BJNEaV1p4O31vTR0oqpnqWpnVe0MzAJec5fuBhaq6snAQGCMS98H9FDVU4DOQG+3GiPAX4BRqtoe2AYMcelDgG2qehwwyuVDRDri7U9+ItAbeExEUuebwaSU31x2jt8hHHBZ91P9DsEki5K+HeEikg9cCDwV4VpDoAcw2SV1BKYBqOoyoK2INFfPbpenljtUvM0FegCvuGvjgYvd733dOe56T5e/LzBRVfep6iq8lR0P22PcGIDG9VNjr+6Lz+xAdoPUiMUkRzp3hI8G7gAiNaT6AdNUdac7XwRcAiAiXYE2eBuChF5xLQQ2AoWqOhvIBbaraom7P7RNIYRtYeiu73D5y9va8BAiMlRE5onIvE2bNlXlc5saYNk36/0OgTHX9+WPA3v7HYZJtnTsCBeRi4CNqjq/nCyXAy+GnY8EclzlcBOwACgBUNVS9zorH+jq+jsq2qawvGsxbW2oqmNVtUBVC/LybDZuusrLaeB3CJQGfXpxbXyT5E2YKi2RLY1uQB8R+RqYCPQQkecARCQX77XQm6HMqrpTVQe7ymEgkAesCi9QVbcD0/H6IzYD2SISmtUevk3hgS0M3fXGwFaqsLWhSV8rvtvidwjc9sS//Q7BJJsqEozt8EPCKg1VHaaq+araFq/z+QNVvdJd7g9MVdW9ofwikh0a/QRcA8xQ1Z0ikici2S5PXeCnwDK3eciHwKXunkHAG+73Ke4cd/0Dl38KMMCNrmoHtAfmxP3Dmxph+sIVfodg0lUKv57ya+2pAXivo8J1ACaISCnwBQdHQrUExrtRTgFgkqpOddfuBCaKyH14r7PGufRxwLMisgKvhTEAQFWXiMgkV34JcIOqlibiA5rqb9fe/X6HYNJUKq89lZRKQ1Wn471WCp13j5BnFt5f/mXTFwMRxxuq6koijH5yLZj+5dxzP3B/TIGbtNa0cX1Wb9rhawwnHt3c1+cbHyjg06unWNiMcGPK8dv+3f0OgQl3VbjFs6mp7PWUMdXPj086xtfnt8xpgDe9yKSbtH89ZYypvNfuGex3CMYnfo2MioVVGsakqAyf9/QwPvHx1VMsrE/DmBR122NvRM9kahxvcp/GdPjBWhrGpKiPvvjG7xCMX1J4IQCrNIypQICU/v/X1FB+tSJiYa+njKnAyGsv8PX5Xa4bRZfrRjHny299jcMkke3cZ0z1dc4ph8039cWvRr/KE1Nn+R2GSYr4rT0lIr3dhnMrROSuCNevFpFNYRvjXROtTKs0jKlARkaAximya97YNz/xOwSTLHHYhMktvfQocD7efkWXu43oynoptDGeqh6291FZVmkYE8Vzd13hdwgmnWjctnvtCqxQ1ZWquh9vtfG+RxqeVRrGRNEyt7HfIRwQTOFJXyaO4rPda0ybzgE/E5HFIvKKiLSOcP0QNnrKmCgCAaH7KccwfdFKv0Oh4IbRBIDWeY3peerxXNfnR2Rk2N9+NU7sfxs0FZF5YedjVXWs+z2WTef+DbyoqvtE5Fd422T3qOiB9l+bMTH467UX+R3CAUHgm007ePq9uZx24xhOu24U23YV+R2WiSMJBmM6gM2hXUbdMTasmKibzqnqFlXd507/CfwwWmxWaRgTg4yMDJo0qOt3GBGVAj3veJK/TypEU3h8v4mR4v1lEMtRsblAexFp5za4G4C3Ed0BItIy7LQPsDRaofZ6ypgY3TOoFzc9mrpLezz/4ec8/+HnAPzk5HY8MPgC6tXJinKXSTVCfJYIUdUSEbkReBfIAJ52G9GNAOap6hTg1yLSB29Tuq3A1dHKtUrDmBh16+TvUumVMWPxKn5866OMveVSCk6I2rdpUk2cWoyq+hbwVpm0P4b9PgwYVpky7fWUMTXY0NGvUFpqC6FUO/EZPZUQCa80RCRDRBaIyFR3PjNs9uFaEZns0nNE5HU39GuOiHRy6a1F5EMRWSoiS0Tk5rCym4hIoYgsdz9zXLqIyMNuFuRiEekSds8gl3+5iAxK9Oc3NUvvguP9DqHSTrtxDF2uG0XB9aMY++YsSoNWiaS0+PVpJEQyWho3E9a5oqpnhWYfArOA19ylu4GFqnoyMBAY49JLgNtVtQNwBnBD2KzGu4BpqtoemObOwZsB2d4dQ4HHwatkgOHA6XgTX4aHKhpjYnH///m7FtWRCCo8MfUTTrthDCNfeN/vcEwFKjF6KukSWmmISD5wIXDY1HQRaYg3HniyS+qI98WPqi4D2opIc1Vdp6qfuvRdeBVQaIJKX7xxxbifF4elT1DPJ0C2GyVwHlCoqltVdRtQCPSO52c2NZuIcHoN6COYNPOzA4shnn3rP5jw7lyKS0r9DssAEOOrqRr6emo0cAeRG1L98FoJO935IuASABHpCrTBG1d8gIi0BU4FZruk5qq6DsD9bObSy5sJGdMMSREZKiLzRGTepk2bYvmcJo3846ZL/A4hrnbtLWH05I84/aaHuf/5QopLrfLwlZKelYaIXARsVNX55WS5HHgx7HwkkCMiC4GbgAV4r6ZC5TUAXgVuCatoyn18hDStIP3QBNWxockyeXl5UR5l0k1mRoAuxx3ldxgJ8epHn3P6jQ/T5bpRvD17id/hpK807dPoBvQRka/xFsrqISLPAYhILl6fwpuhzKq6U1UHu76OgUAesMrlr4VXYTyvqq+FPWNDaHKK+7nRpZc3EzLqDEljYnH3FT/1O4SE+92/3qPLdaO4/9n3/A4l7aTydq8JqzRUdZiq5qtqW7yZiB+o6pXucn9gqqruDeUXkWw3axHgGmCGqu4UEQHGAUtV9aEyj5kChEZADQLeCEsf6EZRnQHscK+v3gV6uZFaOUAvl2ZMpRzTMpeHrz/iBUOrhVc/XkKX60axq2hv9MwmPtLx9VQUAzj01RRAB2CJiCzDG/0UGlrbDbgKr6USGqobGsIyEjhXRJYD57pz8CazrARW4K2ncj2Aqm4F7sWbXj8XGOHSjKm0H59UfSb7xcPZtz/udwjpQRVKg7EdPkjKjHBVnQ5MDzvvHiHPLLwhsmXTPyJyXwSqugXoGSFdgRvKuedp4OmYAjcmigZ1M9m9pyR6xhpiz75i6tau5XcYNV8KryFmM8KNOQKP3FCzRlJFU2r7eSRHCr+esrWnjDkCJx1TM0dRladBXVsAMeEUbyZmirJKw5gj4I3TSA8Z6fNRfaagqbvUi1UaxhyBwvlf+R1C0kz/+3V+h5AeFN86uWNhlYYxR+B3T70ZPVMNUb9uHb9DSB8p3BFulYYxRyBdxk3Zm6kks0rDmJpnf3H6rNH00y6HjYY3CePfyKhYWKVhTBXdM/4dv0NImnuvtsWgk0aBFN7zxOZpGFMFqsrbadAJnteoLh+PuYmsWvb3ZVLZPA1japb7ni/0O4SEaJ1bn6OaNqHfjzvx0y4nEAhYb0byqY2eMqam+c/ilX6HEFdZGcJ/x/yajAx7+eA7BU3heRr2X4gxVdAsu4HfIcTN0U0b8ckjt1iFkUqCGtvhA/uvxJgqePDai/wOIS7q167F5HuH+B2GKSuF+zSs0jCmCvLzsnn7/mto1yLH71CqJBAQrrvoDGaOvtHvUExZqt7oqVgOH1ifhjFV1LxJQ14dfrXfYZiayOZpGGOMiY2ipak7cdQqDWOMSSW2NLoxxphKSechtyKSISILRGSqO58Zttf3WhGZ7NJzROR1EVksInNEpFNYGU+LyEYR+bxM2U1EpFBElrufOS5dRORhEVnhyusSds8gl3+5iAxK9Oc3xpjKUECDGtMRjYj0FpEv3XfhXRGu1xaRl9z12SLSNlqZyRg9dTOwNHSiqmepamdV7QzMAl5zl+4GFqrqycBAYExYGf8CIi1+cxcwTVXbA9PcOcD5ePuNtweGAo+DV8kAw4HTga7A8FBFY4wxKUHdJkyxHBUQkQzgUbzvw47A5SLSsUy2IcA2VT0OGAX8JVp4Ca00RCQfuBB4KsK1hkAPYLJL6oj3xY+qLgPaikhzdz4D2BrhEX2B8e738cDFYekT1PMJkC0iLYHzgEJV3aqq24BCIldGxhjjGy0tjemIoiuwQlVXqup+YCLed2O48O/QV4CeEmU7ykT3aYwG7gAaRrjWD6+VsNOdLwIuAT4Ska5AGyAf2FBB+c1VdR2Aqq4TkWYuvRWwOizfGpdWXvohRGQoXgsFYLeIfAkcDXxbQSx+SdW4wGKrKout8lIlrjZHWsAutr37vr7SNMbsdURkXtj5WFUd636P9H13epn7D+RR1RIR2QHkApvLe2DCKg0RuQjYqKrzRaR7hCyXc2gLZCQwRkQWAp8BC6j6HjeRakqtIP3QBO9f+tjwNBHZpKoFVYwnYVI1LrDYqspiq7xUjasqVDVebz9i+b6L6TsxXCJfT3UD+ojI13jNoh4i8hyAiOTiNZ0O7JWpqjtVdbDr6xgI5AGrojxjg3vthPu50aWvAVqH5csH1laQHovtMeZLtlSNCyy2qrLYKi9V4/JTLN93B/KISCbQmMhdAQckrNJQ1WGqmq+qbYEBwAeqeqW73B+Yqqp7Q/lFJFtEstzpNcCMsFdX5ZkChEZADQLeCEsf6EZRnQHscK+x3gV6uZFaOUAvlxaLHTHmS7ZUjQsstqqy2CovVePy01ygvYi0c9+tA/C+G8OFf4deivc9XWFLw695GgPwXkeF6wBMEJFS4Au8Xn0ARORFoDvQVETWAMNVdZwrY5KIDMF7n9nf3fIWcAGwAigCBgOo6lYRuRfvXybACFWtsFYNMzZ6Fl+kalxgsVWVxVZ5qRqXb1wfxY14fxhnAE+r6hIRGQHMU9UpwDjgWRFZgdfCGBCtXIlSqRhjjDEH2Cq3xhhjYmaVhjHGmJhZpWGMMSZmVmkYY4yJmVUaxhhjYmaVhjHGmJhZpWGMMSZm/x/ezRgbADxngQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#fig=plt.figure()\n",
    "#ax = Axes3D(fig)\n",
    "#ax.scatter(xs=adfmatrix[:, 0], ys=adfmatrix[:, 1], c=y_kmeans, s=1,cmap='viridis')\n",
    "#plt.scatter()\n",
    "#centers = kmeans.cluster_centers_\n",
    "#plt.scatter(centers[:, 0], centers[:, 1], c='black', s=200, alpha=0.5);\n",
    "\n",
    "###\n",
    "sub_adf= adf[::grid]\n",
    "###\n",
    "\n",
    "\n",
    "\n",
    "#ax2 = adf.plot.scatter(x='x', y='y', c=y_kmeans, colormap='viridis')\n",
    "ax2 = sub_adf.plot.scatter(x='x', y='y', c=y_kmeans[::grid], colormap='viridis')                       \n",
    "\n",
    "####\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}