SJER_sample_selection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import thetagrids\n",
    "import h5py\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from sklearn import metrics\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "from dask import delayed, compute\n",
    "from dask.diagnostics import ProgressBar\n",
    "\n",
    "pd.options.mode.chained_assignment = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "site = 'BARR'\n",
    "file_path = f'/media/data/NEON/{site}/filesToStack00200'\n",
    "files = [os.path.join(file_path, f) for f in os.listdir(file_path) if '.h5' in f]\n",
    "save_path = f'/media/data/NEON/{site}/observations/nsae.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook selects observations based on completeness of data then samples from them using two level stratification; first by wind direction sector, then by footprint calculation input parameters.\n",
    "\n",
    "\n",
    "Below we define some functions for later use.  The docstring describe what they do.  The function `find_sectors()` will be used to evenly sample from different wind directions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_valid_observations(f):\n",
    "    '''\n",
    "    Reads footprint statistics from f,\n",
    "    drops bad observations,\n",
    "    then returns dataframe.\n",
    "    '''\n",
    "    # open the hdf\n",
    "    hdf = pd.HDFStore(f)\n",
    "\n",
    "    # get the flux quality flags\n",
    "    qfqm_CO2 = hdf.get(f'{site}/dp04/qfqm/fluxCo2/nsae')\n",
    "    qfqm_H2O = hdf.get(f'{site}/dp04/qfqm/fluxH2o/nsae')\n",
    "    qfqm_T = hdf.get(f'{site}/dp04/qfqm/fluxTemp/nsae')\n",
    "    qfqm_foot = hdf.get(f'{site}/dp04/qfqm/foot/turb')\n",
    "\n",
    "    # Select observations with no bad flags\n",
    "    qfqm_CO2  = qfqm_CO2.loc[qfqm_CO2.qfFinl == 0]\n",
    "    qfqm_H2O  = qfqm_H2O.loc[qfqm_H2O.qfFinl == 0]\n",
    "    qfqm_T    = qfqm_T.loc[qfqm_T.qfFinl == 0]\n",
    "    qfqm_foot = qfqm_foot.loc[qfqm_foot.qfFinl == 0]\n",
    "\n",
    "    # get the footprint input stats\n",
    "    stat = hdf.get(f'{site}/dp04/data/foot/stat/')\n",
    "\n",
    "    # get indices of the dfs from above\n",
    "    istat  = stat.set_index('timeBgn').index\n",
    "    iqfqmC = qfqm_CO2.set_index('timeBgn').index\n",
    "    iqfqmH = qfqm_H2O.set_index('timeBgn').index\n",
    "    iqfqmT = qfqm_T.set_index('timeBgn').index\n",
    "    iqfqmf = qfqm_foot.set_index('timeBgn').index\n",
    "\n",
    "    # keep only entries in stat which correspond to good\n",
    "    # qfqm flags for all variables\n",
    "    good = stat[\n",
    "        (istat.isin(iqfqmC)) &\n",
    "        (istat.isin(iqfqmH)) &\n",
    "        (istat.isin(iqfqmT)) &\n",
    "        (istat.isin(iqfqmf))\n",
    "    ]\n",
    "\n",
    "    hdf.close()\n",
    "\n",
    "    return good\n",
    "\n",
    "\n",
    "def find_sectors(stat, theta=10):\n",
    "    '''\n",
    "    Returns a df of timestamps and sectors of the mean wind direction\n",
    "    '''\n",
    "\n",
    "    # make sure theta goes into 360 an even number of times\n",
    "    if 360 % theta != 0:\n",
    "        while 360 % theta != 0:\n",
    "            theta= theta + 1\n",
    "        print(f'theta has been forced to {theta} for even division of 360')\n",
    "\n",
    "    # set start angle, and empy list\n",
    "    stat['sector'] = theta * (stat.angZaxsErth // theta)\n",
    "\n",
    "    return stat\n",
    "\n",
    "\n",
    "def gmm_cluster(df,k_range=(2, 10), verbose=True):\n",
    "    '''\n",
    "    Performs Gaussian Mixture Model clustering.\n",
    "    df - dataframe containing only the features to\n",
    "         be used in clustering\n",
    "    k_range - tuple, representing range of numbers\n",
    "         of clusters to be tried.\n",
    "    '''\n",
    "    # empty list\n",
    "    silhouette_scores  = []\n",
    "\n",
    "    # make list of range\n",
    "    k_r = list(range(*k_range))\n",
    "\n",
    "    for k in k_r:\n",
    "        # make model\n",
    "        model = GaussianMixture(n_components=k,\n",
    "                                n_init=10,\n",
    "                                init_params='kmeans')\n",
    "        # fit\n",
    "        labels = model.fit_predict(df)\n",
    "\n",
    "        # calculate silhouette score\n",
    "        silhouette_scores.append(\n",
    "            metrics.silhouette_score(df, labels, metric='euclidean')\n",
    "            )\n",
    "\n",
    "        # choose the highest silhouette scoring k\n",
    "        k = k_r[np.argmax(silhouette_scores)]\n",
    "\n",
    "    # now create full model\n",
    "    modelk = GaussianMixture(\n",
    "        n_components=k,\n",
    "        covariance_type='full',\n",
    "        random_state=1\n",
    "        )\n",
    "\n",
    "    cluster = modelk.fit(df)\n",
    "    labels = modelk.predict(df)\n",
    "\n",
    "    if verbose:\n",
    "        print(f'*************** {k} Cluster Model ***************')\n",
    "        print('Weights: ', cluster.weights_)\n",
    "        print('Converged: ', cluster.converged_)\n",
    "        print(' No. of Iterations: ', cluster.n_iter_)\n",
    "        print('Lower Bound: ', cluster.lower_bound_)\n",
    "        print(f'*************************************************')\n",
    "\n",
    "    return labels\n",
    "\n",
    "\n",
    "def fetch_by_sector(files, sector, theta):\n",
    "    '''\n",
    "    Merges footprint stats and nsae observations, performs qc check,\n",
    "    returns df of observations for a given sector.\n",
    "     \n",
    "    sectors is a list or 1d array of sectors\n",
    "    '''\n",
    "\n",
    "    sect_dfs = []\n",
    "\n",
    "    for f in files:\n",
    "            # open the hdf\n",
    "        hdf = pd.HDFStore(f)\n",
    "\n",
    "        # get the flux quality flags\n",
    "        qfqm_CO2 = hdf.get(f'{site}/dp04/qfqm/fluxCo2/nsae')\n",
    "        qfqm_H2O = hdf.get(f'{site}/dp04/qfqm/fluxH2o/nsae')\n",
    "        qfqm_T = hdf.get(f'{site}/dp04/qfqm/fluxTemp/nsae')\n",
    "        qfqm_foot = hdf.get(f'{site}/dp04/qfqm/foot/turb')\n",
    "\n",
    "        # Select observations with no bad flags\n",
    "        qfqm_CO2  = qfqm_CO2.loc[qfqm_CO2.qfFinl == 0]\n",
    "        qfqm_H2O  = qfqm_H2O.loc[qfqm_H2O.qfFinl == 0]\n",
    "        qfqm_T    = qfqm_T.loc[qfqm_T.qfFinl == 0]\n",
    "        qfqm_foot = qfqm_foot.loc[qfqm_foot.qfFinl == 0]\n",
    "\n",
    "        # get the footprint input stats\n",
    "        stat = hdf.get(f'{site}/dp04/data/foot/stat/')\n",
    "\n",
    "        # get indices of the dfs from above\n",
    "        istat  = stat.set_index('timeBgn').index\n",
    "        iqfqmC = qfqm_CO2.set_index('timeBgn').index\n",
    "        iqfqmH = qfqm_H2O.set_index('timeBgn').index\n",
    "        iqfqmT = qfqm_T.set_index('timeBgn').index\n",
    "        iqfqmf = qfqm_foot.set_index('timeBgn').index\n",
    "\n",
    "        # keep only entries in stat which correspond to good\n",
    "        # qfqm flags for all variables\n",
    "        good = stat.loc[\n",
    "            (istat.isin(iqfqmC)) &\n",
    "            (istat.isin(iqfqmH)) &\n",
    "            (istat.isin(iqfqmT)) &\n",
    "            (istat.isin(iqfqmf))\n",
    "        ]\n",
    "\n",
    "        # find the sectors that exist in this df\n",
    "        good = find_sectors(good, theta=theta)\n",
    "\n",
    "        # get fluxes\n",
    "        flux_CO2 = hdf.get(f'{site}/dp04/data/fluxCo2/nsae')\n",
    "        flux_H2O = hdf.get(f'{site}/dp04/data/fluxH2o/nsae')\n",
    "        flux_T   = hdf.get(f'{site}/dp04/data/fluxTemp/nsae')\n",
    "\n",
    "        # change names of flux columns\n",
    "        flux_CO2.columns = ['timeBgn', 'timeEnd', 'flux_CO2']\n",
    "        flux_H2O.columns = ['timeBgn', 'timeEnd', 'flux_H2O']\n",
    "        flux_T.columns = ['timeBgn', 'timeEnd', 'flux_Temp']\n",
    "\n",
    "        # drop timeEnd, set index to timeBgn\n",
    "        flux_CO2 = flux_CO2.drop(columns=['timeEnd']).set_index('timeBgn')\n",
    "        flux_H2O = flux_H2O.drop(columns=['timeEnd']).set_index('timeBgn')\n",
    "        flux_T = flux_T.drop(columns=['timeEnd']).set_index('timeBgn')\n",
    "        good = good.drop(columns=['timeEnd']).set_index('timeBgn')\n",
    "\n",
    "        # merge dfs\n",
    "        dfs = [good, flux_CO2, flux_H2O, flux_T]\n",
    "        data = pd.concat(dfs, axis=1, join='inner')\n",
    "\n",
    "        if sector in data.sector.unique():\n",
    "            # group the data by sector\n",
    "            sect_group = data.groupby('sector')\n",
    "\n",
    "             # get data for sector\n",
    "            df = sect_group.get_group(sector)\n",
    "\n",
    "            # stick it on the list\n",
    "            sect_dfs.append(df)\n",
    "\n",
    "        else:\n",
    "            pass\n",
    "\n",
    "    # concat data from all files pertaining to sector\n",
    "    df = pd.concat(sect_dfs)\n",
    "\n",
    "    # drop rows with nulls\n",
    "    df = df.dropna()\n",
    "\n",
    "    hdf.close()\n",
    "\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# degrees per sector\n",
    "θ = 18\n",
    "\n",
    "# make empty df\n",
    "sectors = pd.DataFrame(columns=['timeBgn', 'sector'])\n",
    "\n",
    "# fill df with timestamps and sectors of valid observations\n",
    "for f in files:\n",
    "    # find footprint stats of the valid observations\n",
    "    stat = find_valid_observations(f)\n",
    "\n",
    "    # find sectors in which observations lie\n",
    "    sects = find_sectors(stat, theta=θ)\n",
    "\n",
    "    # scrunch the latest observations onto the df\n",
    "    sectors = pd.concat([sectors, sects], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "0it [00:00, ?it/s]\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "No objects to concatenate",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m/home/michael/Work/NEON/SJER_sample_selection.ipynb Cell 6'\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=11'>12</a>\u001b[0m     df_list\u001b[39m.\u001b[39mappend(fetch_by_sector(files, s , θ))\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=13'>14</a>\u001b[0m \u001b[39m# concat the list and cluster\u001b[39;00m\n\u001b[0;32m---> <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=14'>15</a>\u001b[0m df \u001b[39m=\u001b[39m pd\u001b[39m.\u001b[39;49mconcat(df_list, axis\u001b[39m=\u001b[39;49m\u001b[39m0\u001b[39;49m)\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=16'>17</a>\u001b[0m \u001b[39m# columns to use in clustering\u001b[39;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=17'>18</a>\u001b[0m cols \u001b[39m=\u001b[39m [\u001b[39m'\u001b[39m\u001b[39mveloYaxsHorSd\u001b[39m\u001b[39m'\u001b[39m,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=18'>19</a>\u001b[0m \u001b[39m'\u001b[39m\u001b[39mveloZaxsHorSd\u001b[39m\u001b[39m'\u001b[39m,\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=19'>20</a>\u001b[0m \u001b[39m'\u001b[39m\u001b[39mveloFric\u001b[39m\u001b[39m'\u001b[39m,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=24'>25</a>\u001b[0m \u001b[39m'\u001b[39m\u001b[39mdistXaxsMax\u001b[39m\u001b[39m'\u001b[39m, \n\u001b[1;32m     <a href='vscode-notebook-cell:/home/michael/Work/NEON/SJER_sample_selection.ipynb#ch0000008?line=25'>26</a>\u001b[0m \u001b[39m'\u001b[39m\u001b[39mdistYaxs90\u001b[39m\u001b[39m'\u001b[39m]\n",
      "File \u001b[0;32m~/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py:311\u001b[0m, in \u001b[0;36mdeprecate_nonkeyword_arguments.<locals>.decorate.<locals>.wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py?line=304'>305</a>\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(args) \u001b[39m>\u001b[39m num_allow_args:\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py?line=305'>306</a>\u001b[0m     warnings\u001b[39m.\u001b[39mwarn(\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py?line=306'>307</a>\u001b[0m         msg\u001b[39m.\u001b[39mformat(arguments\u001b[39m=\u001b[39marguments),\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py?line=307'>308</a>\u001b[0m         \u001b[39mFutureWarning\u001b[39;00m,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py?line=308'>309</a>\u001b[0m         stacklevel\u001b[39m=\u001b[39mstacklevel,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py?line=309'>310</a>\u001b[0m     )\n\u001b[0;32m--> <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/util/_decorators.py?line=310'>311</a>\u001b[0m \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n",
      "File \u001b[0;32m~/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py:346\u001b[0m, in \u001b[0;36mconcat\u001b[0;34m(objs, axis, join, ignore_index, keys, levels, names, verify_integrity, sort, copy)\u001b[0m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=141'>142</a>\u001b[0m \u001b[39m@deprecate_nonkeyword_arguments\u001b[39m(version\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m, allowed_args\u001b[39m=\u001b[39m[\u001b[39m\"\u001b[39m\u001b[39mobjs\u001b[39m\u001b[39m\"\u001b[39m])\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=142'>143</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mconcat\u001b[39m(\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=143'>144</a>\u001b[0m     objs: Iterable[NDFrame] \u001b[39m|\u001b[39m Mapping[Hashable, NDFrame],\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=152'>153</a>\u001b[0m     copy: \u001b[39mbool\u001b[39m \u001b[39m=\u001b[39m \u001b[39mTrue\u001b[39;00m,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=153'>154</a>\u001b[0m ) \u001b[39m-\u001b[39m\u001b[39m>\u001b[39m DataFrame \u001b[39m|\u001b[39m Series:\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=154'>155</a>\u001b[0m     \u001b[39m\"\"\"\u001b[39;00m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=155'>156</a>\u001b[0m \u001b[39m    Concatenate pandas objects along a particular axis with optional set logic\u001b[39;00m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=156'>157</a>\u001b[0m \u001b[39m    along the other axes.\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=343'>344</a>\u001b[0m \u001b[39m    ValueError: Indexes have overlapping values: ['a']\u001b[39;00m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=344'>345</a>\u001b[0m \u001b[39m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=345'>346</a>\u001b[0m     op \u001b[39m=\u001b[39m _Concatenator(\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=346'>347</a>\u001b[0m         objs,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=347'>348</a>\u001b[0m         axis\u001b[39m=\u001b[39;49maxis,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=348'>349</a>\u001b[0m         ignore_index\u001b[39m=\u001b[39;49mignore_index,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=349'>350</a>\u001b[0m         join\u001b[39m=\u001b[39;49mjoin,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=350'>351</a>\u001b[0m         keys\u001b[39m=\u001b[39;49mkeys,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=351'>352</a>\u001b[0m         levels\u001b[39m=\u001b[39;49mlevels,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=352'>353</a>\u001b[0m         names\u001b[39m=\u001b[39;49mnames,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=353'>354</a>\u001b[0m         verify_integrity\u001b[39m=\u001b[39;49mverify_integrity,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=354'>355</a>\u001b[0m         copy\u001b[39m=\u001b[39;49mcopy,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=355'>356</a>\u001b[0m         sort\u001b[39m=\u001b[39;49msort,\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=356'>357</a>\u001b[0m     )\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=358'>359</a>\u001b[0m     \u001b[39mreturn\u001b[39;00m op\u001b[39m.\u001b[39mget_result()\n",
      "File \u001b[0;32m~/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py:403\u001b[0m, in \u001b[0;36m_Concatenator.__init__\u001b[0;34m(self, objs, axis, join, keys, levels, names, ignore_index, verify_integrity, copy, sort)\u001b[0m\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=399'>400</a>\u001b[0m     objs \u001b[39m=\u001b[39m \u001b[39mlist\u001b[39m(objs)\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=401'>402</a>\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(objs) \u001b[39m==\u001b[39m \u001b[39m0\u001b[39m:\n\u001b[0;32m--> <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=402'>403</a>\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mNo objects to concatenate\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=404'>405</a>\u001b[0m \u001b[39mif\u001b[39;00m keys \u001b[39mis\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m    <a href='file:///home/michael/miniconda3/envs/geo3/lib/python3.10/site-packages/pandas/core/reshape/concat.py?line=405'>406</a>\u001b[0m     objs \u001b[39m=\u001b[39m \u001b[39mlist\u001b[39m(com\u001b[39m.\u001b[39mnot_none(\u001b[39m*\u001b[39mobjs))\n",
      "\u001b[0;31mValueError\u001b[0m: No objects to concatenate"
     ]
    }
   ],
   "source": [
    "# empty list\n",
    "df_list = []\n",
    "\n",
    "# get the sectors\n",
    "sc = sectors.sector.unique()\n",
    "sc.sort()\n",
    "\n",
    "# for each sector (this cannot be || because hdfs are open)\n",
    "for s in tqdm(sc):\n",
    "    \n",
    "    # append df for that sector to list \n",
    "    df_list.append(fetch_by_sector(files, s , θ))\n",
    "\n",
    "# concat the list and cluster\n",
    "df = pd.concat(df_list, axis=0)\n",
    "\n",
    "# columns to use in clustering\n",
    "cols = ['veloYaxsHorSd',\n",
    "'veloZaxsHorSd',\n",
    "'veloFric',\n",
    "'distZaxsMeasDisp',\n",
    "'distZaxsRgh',\n",
    "'distZaxsAbl',\n",
    "'distXaxs90',\n",
    "'distXaxsMax', \n",
    "'distYaxs90']\n",
    "\n",
    "# now cluster based on footprint inputs, and note n clusters\n",
    "df['label'] = gmm_cluster(df[cols])\n",
    "n_clust = df.label.unique().shape[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.0      9\n",
       "72.0     14\n",
       "108.0    16\n",
       "90.0     17\n",
       "0.0      21\n",
       "36.0     22\n",
       "54.0     23\n",
       "342.0    23\n",
       "126.0    32\n",
       "324.0    38\n",
       "Name: sector, dtype: int64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# look at some of smallest sector counts\n",
    "sectors.sector.value_counts().sort_values().head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the above output we can see that there is descent representation across sectors. Let's choose 20 as a sample size from each sector. There will be three underrepresented sectors, but they will not be hugely underrepresented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 20/20 [00:00<00:00, 15363.75it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[########################################] | 100% Completed |  0.1s\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# number of samples\n",
    "n_samples = 20\n",
    "\n",
    "# function to go in dask loop to get samples\n",
    "@delayed()\n",
    "def lil_func(s, sector_groups):\n",
    "    ddd = []\n",
    "\n",
    "    d = sector_groups.get_group(s)\n",
    "\n",
    "    # find the smallest cluseter count in d\n",
    "    x = np.inf\n",
    "    for i in range(n_clust):\n",
    "        x = min(x, len(d.loc[d.label == i]))\n",
    "        # now many to take from each cluster\n",
    "        n_s = min(n_samples // 2, x)\n",
    "    \n",
    "    # sample from each cluster\n",
    "    for i in range(n_clust):\n",
    "        ddd.append(d.loc[d.label == i].sample(n_s))\n",
    "    d = pd.concat(ddd, axis=0)\n",
    "    return d\n",
    "\n",
    "# now bust the df up by sector again\n",
    "sector_groups = df.groupby('sector')\n",
    "\n",
    "# for each sector\n",
    "lazy = []\n",
    "for s in tqdm(sc):\n",
    "    lazy.append(lil_func(s, sector_groups))\n",
    "\n",
    "with ProgressBar():\n",
    "\n",
    "    to_concat = compute(*lazy)\n",
    "\n",
    "    d = pd.concat(to_concat, axis=0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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>angZaxsErth</th>\n",
       "      <th>distReso</th>\n",
       "      <th>veloYaxsHorSd</th>\n",
       "      <th>veloZaxsHorSd</th>\n",
       "      <th>veloFric</th>\n",
       "      <th>distZaxsMeasDisp</th>\n",
       "      <th>distZaxsRgh</th>\n",
       "      <th>distZaxsAbl</th>\n",
       "      <th>distXaxs90</th>\n",
       "      <th>distXaxsMax</th>\n",
       "      <th>distYaxs90</th>\n",
       "      <th>sector</th>\n",
       "      <th>flux_CO2</th>\n",
       "      <th>flux_H2O</th>\n",
       "      <th>flux_Temp</th>\n",
       "      <th>label</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>timeBgn</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2019-09-04T01:30:00.000Z</th>\n",
       "      <td>33.300662</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.987935</td>\n",
       "      <td>0.547851</td>\n",
       "      <td>0.413994</td>\n",
       "      <td>42.41</td>\n",
       "      <td>3.882275</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>560.0</td>\n",
       "      <td>280.0</td>\n",
       "      <td>240.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>3.881349</td>\n",
       "      <td>57.971081</td>\n",
       "      <td>-38.694133</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-08-30T03:30:00.000Z</th>\n",
       "      <td>22.813947</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.315621</td>\n",
       "      <td>0.230000</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>42.41</td>\n",
       "      <td>4.241000</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>600.0</td>\n",
       "      <td>280.0</td>\n",
       "      <td>200.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>2.401831</td>\n",
       "      <td>-20.109489</td>\n",
       "      <td>-26.950858</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-08-23T15:30:00.000Z</th>\n",
       "      <td>38.385578</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.416736</td>\n",
       "      <td>0.321403</td>\n",
       "      <td>0.215207</td>\n",
       "      <td>42.41</td>\n",
       "      <td>4.241000</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>520.0</td>\n",
       "      <td>240.0</td>\n",
       "      <td>200.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>-13.989782</td>\n",
       "      <td>185.337514</td>\n",
       "      <td>137.965709</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-09-24T06:00:00.000Z</th>\n",
       "      <td>50.010893</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.540548</td>\n",
       "      <td>0.230000</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>42.41</td>\n",
       "      <td>2.503394</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>760.0</td>\n",
       "      <td>320.0</td>\n",
       "      <td>320.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>-5.419283</td>\n",
       "      <td>-0.765701</td>\n",
       "      <td>2.163597</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-09-05T20:00:00.000Z</th>\n",
       "      <td>85.135861</td>\n",
       "      <td>40.0</td>\n",
       "      <td>1.207982</td>\n",
       "      <td>0.680086</td>\n",
       "      <td>0.672338</td>\n",
       "      <td>42.41</td>\n",
       "      <td>4.241000</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>680.0</td>\n",
       "      <td>280.0</td>\n",
       "      <td>280.0</td>\n",
       "      <td>72.0</td>\n",
       "      <td>-8.367664</td>\n",
       "      <td>251.247036</td>\n",
       "      <td>179.738166</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-08-24T11:30:00.000Z</th>\n",
       "      <td>337.174517</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.359577</td>\n",
       "      <td>0.230000</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>42.41</td>\n",
       "      <td>0.965674</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>440.0</td>\n",
       "      <td>200.0</td>\n",
       "      <td>324.0</td>\n",
       "      <td>0.074285</td>\n",
       "      <td>2.722504</td>\n",
       "      <td>-9.942192</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-10-25T04:00:00.000Z</th>\n",
       "      <td>324.195835</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.373068</td>\n",
       "      <td>0.230000</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>42.41</td>\n",
       "      <td>0.088155</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>1640.0</td>\n",
       "      <td>640.0</td>\n",
       "      <td>200.0</td>\n",
       "      <td>324.0</td>\n",
       "      <td>3.449447</td>\n",
       "      <td>-21.971639</td>\n",
       "      <td>-11.122484</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-07-05T08:00:00.000Z</th>\n",
       "      <td>325.018644</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.465912</td>\n",
       "      <td>0.230000</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>42.41</td>\n",
       "      <td>4.241000</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>600.0</td>\n",
       "      <td>280.0</td>\n",
       "      <td>320.0</td>\n",
       "      <td>324.0</td>\n",
       "      <td>0.976421</td>\n",
       "      <td>-10.823095</td>\n",
       "      <td>-3.762834</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-09-04T02:00:00.000Z</th>\n",
       "      <td>347.345615</td>\n",
       "      <td>40.0</td>\n",
       "      <td>1.345506</td>\n",
       "      <td>0.381645</td>\n",
       "      <td>0.334950</td>\n",
       "      <td>42.41</td>\n",
       "      <td>4.241000</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>600.0</td>\n",
       "      <td>280.0</td>\n",
       "      <td>520.0</td>\n",
       "      <td>342.0</td>\n",
       "      <td>16.087157</td>\n",
       "      <td>106.646516</td>\n",
       "      <td>-88.122698</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-10-24T15:30:00.000Z</th>\n",
       "      <td>344.744146</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.557173</td>\n",
       "      <td>0.230000</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>42.41</td>\n",
       "      <td>1.871786</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>840.0</td>\n",
       "      <td>360.0</td>\n",
       "      <td>320.0</td>\n",
       "      <td>342.0</td>\n",
       "      <td>-10.371942</td>\n",
       "      <td>44.400944</td>\n",
       "      <td>89.807029</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>144 rows × 16 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                          angZaxsErth  distReso  veloYaxsHorSd  veloZaxsHorSd  \\\n",
       "timeBgn                                                                         \n",
       "2019-09-04T01:30:00.000Z    33.300662      40.0       0.987935       0.547851   \n",
       "2019-08-30T03:30:00.000Z    22.813947      40.0       0.315621       0.230000   \n",
       "2019-08-23T15:30:00.000Z    38.385578      40.0       0.416736       0.321403   \n",
       "2019-09-24T06:00:00.000Z    50.010893      40.0       0.540548       0.230000   \n",
       "2019-09-05T20:00:00.000Z    85.135861      40.0       1.207982       0.680086   \n",
       "...                               ...       ...            ...            ...   \n",
       "2019-08-24T11:30:00.000Z   337.174517      40.0       0.359577       0.230000   \n",
       "2019-10-25T04:00:00.000Z   324.195835      40.0       0.373068       0.230000   \n",
       "2019-07-05T08:00:00.000Z   325.018644      40.0       0.465912       0.230000   \n",
       "2019-09-04T02:00:00.000Z   347.345615      40.0       1.345506       0.381645   \n",
       "2019-10-24T15:30:00.000Z   344.744146      40.0       0.557173       0.230000   \n",
       "\n",
       "                          veloFric  distZaxsMeasDisp  distZaxsRgh  \\\n",
       "timeBgn                                                             \n",
       "2019-09-04T01:30:00.000Z  0.413994             42.41     3.882275   \n",
       "2019-08-30T03:30:00.000Z  0.200000             42.41     4.241000   \n",
       "2019-08-23T15:30:00.000Z  0.215207             42.41     4.241000   \n",
       "2019-09-24T06:00:00.000Z  0.200000             42.41     2.503394   \n",
       "2019-09-05T20:00:00.000Z  0.672338             42.41     4.241000   \n",
       "...                            ...               ...          ...   \n",
       "2019-08-24T11:30:00.000Z  0.200000             42.41     0.965674   \n",
       "2019-10-25T04:00:00.000Z  0.200000             42.41     0.088155   \n",
       "2019-07-05T08:00:00.000Z  0.200000             42.41     4.241000   \n",
       "2019-09-04T02:00:00.000Z  0.334950             42.41     4.241000   \n",
       "2019-10-24T15:30:00.000Z  0.200000             42.41     1.871786   \n",
       "\n",
       "                          distZaxsAbl  distXaxs90  distXaxsMax  distYaxs90  \\\n",
       "timeBgn                                                                      \n",
       "2019-09-04T01:30:00.000Z       1000.0       560.0        280.0       240.0   \n",
       "2019-08-30T03:30:00.000Z       1000.0       600.0        280.0       200.0   \n",
       "2019-08-23T15:30:00.000Z       1000.0       520.0        240.0       200.0   \n",
       "2019-09-24T06:00:00.000Z       1000.0       760.0        320.0       320.0   \n",
       "2019-09-05T20:00:00.000Z       1000.0       680.0        280.0       280.0   \n",
       "...                               ...         ...          ...         ...   \n",
       "2019-08-24T11:30:00.000Z       1000.0      1000.0        440.0       200.0   \n",
       "2019-10-25T04:00:00.000Z       1000.0      1640.0        640.0       200.0   \n",
       "2019-07-05T08:00:00.000Z       1000.0       600.0        280.0       320.0   \n",
       "2019-09-04T02:00:00.000Z       1000.0       600.0        280.0       520.0   \n",
       "2019-10-24T15:30:00.000Z       1000.0       840.0        360.0       320.0   \n",
       "\n",
       "                          sector   flux_CO2    flux_H2O   flux_Temp  label  \n",
       "timeBgn                                                                     \n",
       "2019-09-04T01:30:00.000Z    18.0   3.881349   57.971081  -38.694133      0  \n",
       "2019-08-30T03:30:00.000Z    18.0   2.401831  -20.109489  -26.950858      1  \n",
       "2019-08-23T15:30:00.000Z    36.0 -13.989782  185.337514  137.965709      0  \n",
       "2019-09-24T06:00:00.000Z    36.0  -5.419283   -0.765701    2.163597      1  \n",
       "2019-09-05T20:00:00.000Z    72.0  -8.367664  251.247036  179.738166      0  \n",
       "...                          ...        ...         ...         ...    ...  \n",
       "2019-08-24T11:30:00.000Z   324.0   0.074285    2.722504   -9.942192      1  \n",
       "2019-10-25T04:00:00.000Z   324.0   3.449447  -21.971639  -11.122484      1  \n",
       "2019-07-05T08:00:00.000Z   324.0   0.976421  -10.823095   -3.762834      1  \n",
       "2019-09-04T02:00:00.000Z   342.0  16.087157  106.646516  -88.122698      0  \n",
       "2019-10-24T15:30:00.000Z   342.0 -10.371942   44.400944   89.807029      1  \n",
       "\n",
       "[144 rows x 16 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQdUlEQVR4nO3df6zddX3H8edrBdyGJKi9Ij+9zBASNFKbm6phIziVtYWILmZrs0zmSKoGEk22ZHUmzv1Xt+gSxUDqaMDFFbcoStKqEOOCJiresgJliFRWQy2hF4kgwcxV3/vjfjuv13Pubc857TmXz/ORnJzv9/P9fL+f9/309pXv/d7z/d5UFZKkF77fGncBkqSTw8CXpEYY+JLUCANfkhph4EtSI04ZdwG9rF69uqanp8ddhiStGHv27HmqqqaW6jORgT89Pc3s7Oy4y5CkFSPJD5fr4yUdSWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqxETeaStpckxv3TWWcQ9su2os476QeYYvSY1Y9gw/yQ7gauBwVb2ma/sccHHX5UzgJ1W1pse+B4CfAr8AjlTVzEiqliQdt2O5pHMrcCPwmaMNVfWnR5eTfAx4Zon931RVTw1aoCRpNJYN/Kq6J8l0r21JAvwJ8IcjrkuSNGLDXsP/A+DJqnq0z/YC7kqyJ8mWpQ6UZEuS2SSzc3NzQ5YlSVps2MDfDOxcYvtlVbUW2ABcn+Tyfh2rantVzVTVzNTUks/wlyQNYODAT3IK8MfA5/r1qapD3fth4A5g3aDjSZKGM8wZ/luA71XVwV4bk5ye5Iyjy8CVwL4hxpMkDWHZwE+yE/gWcHGSg0mu6zZtYtHlnCTnJNndrZ4FfDPJ/cC9wK6q+sroSpckHY9j+ZTO5j7tf9Gj7RCwsVt+DLh0yPokSSPinbaS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJasSx/E1bqafprbsG3vfAtqtGWMnJ0+LXrBcOz/AlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWrEsoGfZEeSw0n2LWj7SJIfJdnbvTb22Xd9kkeS7E+ydZSFS5KOz7Gc4d8KrO/R/k9VtaZ77V68Mckq4FPABuASYHOSS4YpVpI0uGUDv6ruAZ4e4NjrgP1V9VhV/Ry4HbhmgONIkkZgmGv4NyR5oLvk85Ie288FHl+wfrBr6ynJliSzSWbn5uaGKEuS1MuggX8T8CpgDfAE8LEefdKjrfodsKq2V9VMVc1MTU0NWJYkqZ+BAr+qnqyqX1TVL4FPM3/5ZrGDwPkL1s8DDg0yniRpeAMFfpKzF6y+A9jXo9t3gYuSXJjkNGATcOcg40mShrfs45GT7ASuAFYnOQj8HXBFkjXMX6I5ALyn63sO8M9VtbGqjiS5AfgqsArYUVUPnYgvQpK0vGUDv6o292i+pU/fQ8DGBeu7gd/4yKYk6eTzTltJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUiGU/lqnJN71118D7Hth21QgrkTTJPMOXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqxLKBn2RHksNJ9i1o+8ck30vyQJI7kpzZZ98DSR5MsjfJ7AjrliQdp2M5w78VWL+o7W7gNVX1WuD7wAeX2P9NVbWmqmYGK1GSNArLBn5V3QM8vajtrqo60q1+GzjvBNQmSRqhUVzD/0vgy322FXBXkj1Jtix1kCRbkswmmZ2bmxtBWZKkhYYK/CQfAo4An+3T5bKqWgtsAK5Pcnm/Y1XV9qqaqaqZqampYcqSJPUwcOAnuRa4GvizqqpefarqUPd+GLgDWDfoeJKk4QwU+EnWA38DvK2qnu/T5/QkZxxdBq4E9vXqK0k68Y7lY5k7gW8BFyc5mOQ64EbgDODu7iOXN3d9z0myu9v1LOCbSe4H7gV2VdVXTshXIUla1inLdaiqzT2ab+nT9xCwsVt+DLh0qOokSSPjnbaS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWrEsoGfZEeSw0n2LWh7aZK7kzzavb+kz77rkzySZH+SraMsXJJ0fI7lDP9WYP2itq3A16rqIuBr3fqvSbIK+BSwAbgE2JzkkqGqlSQNbNnAr6p7gKcXNV8D3NYt3wa8vceu64D9VfVYVf0cuL3bT5I0BqcMuN9ZVfUEQFU9keTlPfqcCzy+YP0g8Pp+B0yyBdgCcMEFFwxYljS5prfuGnjfA9uuGtvYeuE4kb+0TY+26te5qrZX1UxVzUxNTZ3AsiSpTYMG/pNJzgbo3g/36HMQOH/B+nnAoQHHkyQNadDAvxO4tlu+FvhSjz7fBS5KcmGS04BN3X6SpDE4lo9l7gS+BVyc5GCS64BtwFuTPAq8tVsnyTlJdgNU1RHgBuCrwMPAv1XVQyfmy5AkLWfZX9pW1eY+m97co+8hYOOC9d3A7oGrkySNjHfaSlIjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiOW/YtXOvGmt+4adwmSGuAZviQ1YuDAT3Jxkr0LXs8m+cCiPlckeWZBnw8PXbEkaSADX9KpqkeANQBJVgE/Au7o0fUbVXX1oONIkkZjVJd03gz8oKp+OKLjSZJGbFSBvwnY2WfbG5Pcn+TLSV49ovEkScdp6MBPchrwNuDfe2y+D3hlVV0KfBL44hLH2ZJkNsns3NzcsGVJkhYZxRn+BuC+qnpy8YaqeraqnuuWdwOnJlnd6yBVtb2qZqpqZmpqagRlSZIWGkXgb6bP5Zwkr0iSbnldN96PRzCmJOk4DXXjVZLfBd4KvGdB23sBqupm4J3A+5IcAX4GbKqqGmZMSdJghgr8qnoeeNmitpsXLN8I3DjMGJKk0fBOW0lqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1YqiHp0k6Oaa37hp3CU0ZZr4PbLtqhJWMlmf4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1YqjAT3IgyYNJ9iaZ7bE9ST6RZH+SB5KsHWY8SdLgRnHj1Zuq6qk+2zYAF3Wv1wM3de+SpJPsRF/SuQb4TM37NnBmkrNP8JiSpB6GDfwC7kqyJ8mWHtvPBR5fsH6wa/sNSbYkmU0yOzc3N2RZkqTFhg38y6pqLfOXbq5Pcvmi7emxT/U6UFVtr6qZqpqZmpoasixJ0mJDBX5VHereDwN3AOsWdTkInL9g/Tzg0DBjSpIGM3DgJzk9yRlHl4ErgX2Lut0JvKv7tM4bgGeq6omBq5UkDWyYT+mcBdyR5Ohx/rWqvpLkvQBVdTOwG9gI7AeeB949XLmSpEENHPhV9RhwaY/2mxcsF3D9oGNIkkbHO20lqREGviQ1wsCXpEYY+JLUCANfkhph4EtSI0bxtEytYNNbd63IcQ9su2pElWhS+T0yep7hS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNWLgwE9yfpKvJ3k4yUNJ3t+jzxVJnkmyt3t9eLhyJUmDGubxyEeAv6qq+5KcAexJcndV/deift+oqquHGEeSNAIDn+FX1RNVdV+3/FPgYeDcURUmSRqtkVzDTzINvA74To/Nb0xyf5IvJ3n1EsfYkmQ2yezc3NwoypIkLTB04Cd5MfB54ANV9eyizfcBr6yqS4FPAl/sd5yq2l5VM1U1MzU1NWxZkqRFhgr8JKcyH/afraovLN5eVc9W1XPd8m7g1CSrhxlTkjSYYT6lE+AW4OGq+nifPq/o+pFkXTfejwcdU5I0uGE+pXMZ8OfAg0n2dm1/C1wAUFU3A+8E3pfkCPAzYFNV1RBjSpIGNHDgV9U3gSzT50bgxkHHkCSNzjBn+BNpeuuusYx7YNtVYxm3VcP8O/tv1YZxZcEkf2/6aAVJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEUMFfpL1SR5Jsj/J1h7bk+QT3fYHkqwdZjxJ0uAGDvwkq4BPARuAS4DNSS5Z1G0DcFH32gLcNOh4kqThDHOGvw7YX1WPVdXPgduBaxb1uQb4TM37NnBmkrOHGFOSNKBThtj3XODxBesHgdcfQ59zgScWHyzJFuZ/CgB4LskjQ9R20uWjrAaeGncdA2qq9nz0BFVy/Jqa9wkysbUfw/fmUrW/crmdhwn89GirAfrMN1ZtB7YPUc9YJZmtqplx1zEIax8Pax+Plmsf5pLOQeD8BevnAYcG6CNJOgmGCfzvAhcluTDJacAm4M5Ffe4E3tV9WucNwDNV9RuXcyRJJ97Al3Sq6kiSG4CvAquAHVX1UJL3dttvBnYDG4H9wPPAu4cveWKt2MtRWPu4WPt4NFt7qnpeUpckvcB4p60kNcLAl6RGGPgDSHIgyYNJ9iaZ7dpemuTuJI927y8Zd50ASXYkOZxk34K2vrUm+WD3KIxHkvzReKr+/1p61f6RJD/q5n5vko0Ltk1S7ecn+XqSh5M8lOT9XfvEz/0StU/83Cf57ST3Jrm/q/3vu/aVMO/9ah/dvFeVr+N8AQeA1Yva/gHY2i1vBT467jq7Wi4H1gL7lquV+Udk3A+8CLgQ+AGwasJq/wjw1z36TlrtZwNru+UzgO93NU783C9R+8TPPfP3/ry4Wz4V+A7whhUy7/1qH9m8e4Y/OtcAt3XLtwFvH18pv1JV9wBPL2ruV+s1wO1V9T9V9d/Mf7pq3cmos5c+tfczabU/UVX3dcs/BR5m/i7ziZ/7JWrvZ5Jqr6p6rls9tXsVK2Pe+9Xez3HXbuAPpoC7kuzpHgkBcFZ19xh07y8fW3XL61drv0dhTJobuqev7ljwo/nE1p5kGngd82dsK2ruF9UOK2Duk6xKshc4DNxdVStm3vvUDiOadwN/MJdV1VrmnwZ6fZLLx13QiBzzozDG6CbgVcAa5p/J9LGufSJrT/Ji4PPAB6rq2aW69mgba/09al8Rc19Vv6iqNczf2b8uyWuW6L4Sah/ZvBv4A6iqQ937YeAO5n+MejLdk0C798Pjq3BZ/Wqd+EdhVNWT3X+KXwKf5lc/wk5c7UlOZT4wP1tVX+iaV8Tc96p9Jc09QFX9BPgPYD0rZN6PWlj7KOfdwD9OSU5PcsbRZeBKYB/zj5G4tut2LfCl8VR4TPrVeiewKcmLklzI/N8xuHcM9fWVX3+89juYn3uYsNqTBLgFeLiqPr5g08TPfb/aV8LcJ5lKcma3/DvAW4DvsTLmvWftI533cfw2eiW/gN9j/jfj9wMPAR/q2l8GfA14tHt/6bhr7erayfyPgf/L/BnBdUvVCnyI+d/2PwJsmMDa/wV4EHig+4Y/e0Jr/33mf7x+ANjbvTauhLlfovaJn3vgtcB/djXuAz7cta+Eee9X+8jm3UcrSFIjvKQjSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1Ij/g+7mBXEzPwehgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(d.sector.values, bins=sc.shape[0]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save the  data, not big, just use csv\n",
    "d.to_csv(save_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "0b20d2e798ab320dcd9e3022b3f23321d95437ce0beeb6a34c6f05479d8e44a2"
  },
  "kernelspec": {
   "display_name": "Python 3.10.2 ('geo3')",
   "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.10.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}