Read_and_Plot_PebbleCounts_CSV.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "059277a2",
   "metadata": {},
   "source": [
    "This notebook is a short tutorial on how to read in a PebbleCounts .csv file. You could always open these in a spreadsheet software like Excel, but why not learn some Python. \n",
    "\n",
    "-Ben Purinton; August 18th, 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2ba9c01e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from ast import literal_eval\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb1d37b0",
   "metadata": {},
   "source": [
    "# The File\n",
    "\n",
    "The .csv files have some different sections for the different data output by PebbleCounts. These sections are:\n",
    "\n",
    "1. **PebbleCounts Parameters**: One row for the image containing the parameters used in the PebbleCounts run which generated the file.\n",
    "\n",
    "\n",
    "2. **Image Details**: One row for the image with the \"perc. not meas.\" (percent not measured) and \"perc. background color\" (percent background color; i.e., areal percentage of sand or smooth compacted surfaces).\n",
    "\n",
    "\n",
    "3. **Pebble Details**: N-rows, where N is the number of pebbles counted, containing the grain-size data.\n",
    "\n",
    "The first two pieces of information can be read in as a Python dictionary, and the last item (grain-size data) we will read in as a Pandas dataframe.\n",
    "\n",
    "Here is what the file looks like when you open it with a spreadsheet software:\n",
    "\n",
    "![](docs/figs/output_csv.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc1e0a00",
   "metadata": {},
   "source": [
    "# PebbleCounts (k-means manual selection; KMS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b7cd8ba",
   "metadata": {},
   "source": [
    "First we will open and plot output from the KMS (`PebbleCounts.py`) algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d7464b7",
   "metadata": {},
   "source": [
    "## Select CSV\n",
    "\n",
    "**Note:** If you don't have the file below in the `example_data` folder, then you need to run `PebbleCounts.py` to produce it. Or just follow along with the notebook without it and copy-paste the code you need into your own projects. This file was generated with the following command (and manual clicking):\n",
    "\n",
    "```\n",
    "python PebbleCounts.py -im example_data\\ortho_resolution_1.2mmPerPix.tif -ortho y -maxGS 0.3 \n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "11bcb11c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fn = 'example_data/ortho_resolution_1.2mmPerPix_PebbleCounts_CSV.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24f63b94",
   "metadata": {},
   "source": [
    "## Read PebbleCounts Parameters\n",
    "\n",
    "First we'll read the first line with the parameters from the run.\n",
    "\n",
    "**Note:** The `skiprows` and `nrows` parameter are very important and should not be changed. This is standardized across all PebbleCounts .csv outputs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3c561b6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = pd.read_csv(fn, skiprows=1, nrows=1).T.to_dict()[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad586869",
   "metadata": {},
   "source": [
    "This is now a Python dictionary, which you can read items from"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dd9e8c05",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'otsu_threshold': 50.0,\n",
       " 'maxGS': 0.3,\n",
       " 'cutoff': 20,\n",
       " 'min_sz_factors': '[50, 5, 1]',\n",
       " 'win_sz_factors': '[10, 3, 2]',\n",
       " 'improvement_ths': '[0.01, 0.1, 0.1]',\n",
       " 'coordinate_scales': '[0.5, 0.5, 0.5]',\n",
       " 'overlaps': '[0.5, 0.3, 0.1]',\n",
       " 'first_nl_denoise': 5,\n",
       " 'nl_means_chroma_filts': '[3, 2, 1]',\n",
       " 'bilat_filt_szs': '[9, 5, 3]',\n",
       " 'tophat_th': 90,\n",
       " 'sobel_th': 90,\n",
       " 'canny_sig': 2}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "43301257",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90\n"
     ]
    }
   ],
   "source": [
    "print(params['tophat_th'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22213535",
   "metadata": {},
   "source": [
    "One thing to note: We need the function imported above `literal_eval` to read the list-type items (e.g., `min_sz_factors`) in the above dictionary. When the algorithm saves the data to a .csv file, the list items are converted into strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ade50fa4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "str"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(params['min_sz_factors'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "650de1eb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[50, 5, 1]\n",
      "<class 'list'>\n"
     ]
    }
   ],
   "source": [
    "print(literal_eval(params['min_sz_factors']))\n",
    "print(type(literal_eval(params['min_sz_factors'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b0dd07e",
   "metadata": {},
   "source": [
    "## Read Image Details\n",
    "\n",
    "We can also read in the image details (percent unmeasured and percent background color). \n",
    "\n",
    "**Note:** The `skiprows` and `nrows` parameter are very important and should not be changed. This is standardized across all PebbleCounts .csv outputs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5ed78432",
   "metadata": {},
   "outputs": [],
   "source": [
    "tile_details = pd.read_csv(fn, skiprows=5, nrows=1).T.to_dict()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2a4fce95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'perc. not meas.': 80.56790000000001, 'perc. background color': 17.446}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tile_details"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4aa5a3d",
   "metadata": {},
   "source": [
    "## Read and Plot Pebbles\n",
    "\n",
    "And finally, the main event: grain-size data. Now we will use a pandas dataframe.\n",
    "\n",
    "**Note:** The `skiprows` and `nrows` parameter are very important and should not be changed. This is standardized across all PebbleCounts .csv outputs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6e7473d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "pebbles = pd.read_csv(fn, skiprows=9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a1c25436",
   "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>UTM X (m)</th>\n",
       "      <th>UTM Y (m)</th>\n",
       "      <th>a (px)</th>\n",
       "      <th>b (px)</th>\n",
       "      <th>a (m)</th>\n",
       "      <th>b (m)</th>\n",
       "      <th>area (px)</th>\n",
       "      <th>area (m2)</th>\n",
       "      <th>orientation</th>\n",
       "      <th>ellipse area (px)</th>\n",
       "      <th>perc. diff. area</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>836644.2878</td>\n",
       "      <td>7242472.866</td>\n",
       "      <td>87.024678</td>\n",
       "      <td>55.856531</td>\n",
       "      <td>0.104462</td>\n",
       "      <td>0.067049</td>\n",
       "      <td>3401</td>\n",
       "      <td>0.004900</td>\n",
       "      <td>-2.309977</td>\n",
       "      <td>3807.285687</td>\n",
       "      <td>10.671269</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>836644.5051</td>\n",
       "      <td>7242472.819</td>\n",
       "      <td>91.028929</td>\n",
       "      <td>55.045931</td>\n",
       "      <td>0.109269</td>\n",
       "      <td>0.066076</td>\n",
       "      <td>3433</td>\n",
       "      <td>0.004947</td>\n",
       "      <td>-2.527877</td>\n",
       "      <td>3924.675289</td>\n",
       "      <td>12.527795</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>836644.6023</td>\n",
       "      <td>7242472.809</td>\n",
       "      <td>46.619508</td>\n",
       "      <td>38.267768</td>\n",
       "      <td>0.055961</td>\n",
       "      <td>0.045936</td>\n",
       "      <td>1288</td>\n",
       "      <td>0.001856</td>\n",
       "      <td>-0.218997</td>\n",
       "      <td>1397.332929</td>\n",
       "      <td>7.824401</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>836644.5159</td>\n",
       "      <td>7242472.703</td>\n",
       "      <td>150.841046</td>\n",
       "      <td>69.113499</td>\n",
       "      <td>0.181065</td>\n",
       "      <td>0.082962</td>\n",
       "      <td>7258</td>\n",
       "      <td>0.010458</td>\n",
       "      <td>-2.214070</td>\n",
       "      <td>8165.475866</td>\n",
       "      <td>11.113570</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>836644.1702</td>\n",
       "      <td>7242472.628</td>\n",
       "      <td>114.262068</td>\n",
       "      <td>108.564897</td>\n",
       "      <td>0.137157</td>\n",
       "      <td>0.130318</td>\n",
       "      <td>8869</td>\n",
       "      <td>0.012779</td>\n",
       "      <td>-3.002590</td>\n",
       "      <td>9716.068830</td>\n",
       "      <td>8.718226</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",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64</th>\n",
       "      <td>836644.9420</td>\n",
       "      <td>7242472.026</td>\n",
       "      <td>68.150933</td>\n",
       "      <td>55.956297</td>\n",
       "      <td>0.081806</td>\n",
       "      <td>0.067168</td>\n",
       "      <td>2631</td>\n",
       "      <td>0.003791</td>\n",
       "      <td>-1.690621</td>\n",
       "      <td>2986.894275</td>\n",
       "      <td>11.915195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65</th>\n",
       "      <td>836645.0608</td>\n",
       "      <td>7242472.013</td>\n",
       "      <td>49.360133</td>\n",
       "      <td>36.526219</td>\n",
       "      <td>0.059251</td>\n",
       "      <td>0.043845</td>\n",
       "      <td>1106</td>\n",
       "      <td>0.001594</td>\n",
       "      <td>-2.726250</td>\n",
       "      <td>1412.147675</td>\n",
       "      <td>21.679579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>836644.9168</td>\n",
       "      <td>7242471.961</td>\n",
       "      <td>80.284876</td>\n",
       "      <td>42.840646</td>\n",
       "      <td>0.096372</td>\n",
       "      <td>0.051425</td>\n",
       "      <td>2194</td>\n",
       "      <td>0.003161</td>\n",
       "      <td>-2.755353</td>\n",
       "      <td>2693.945657</td>\n",
       "      <td>18.558120</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>836645.0932</td>\n",
       "      <td>7242471.847</td>\n",
       "      <td>37.748005</td>\n",
       "      <td>31.099214</td>\n",
       "      <td>0.045312</td>\n",
       "      <td>0.037331</td>\n",
       "      <td>819</td>\n",
       "      <td>0.001180</td>\n",
       "      <td>-2.054698</td>\n",
       "      <td>919.480448</td>\n",
       "      <td>10.927959</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>68</th>\n",
       "      <td>836645.2541</td>\n",
       "      <td>7242472.284</td>\n",
       "      <td>34.972655</td>\n",
       "      <td>22.441033</td>\n",
       "      <td>0.041980</td>\n",
       "      <td>0.026938</td>\n",
       "      <td>562</td>\n",
       "      <td>0.000810</td>\n",
       "      <td>-0.949405</td>\n",
       "      <td>614.710368</td>\n",
       "      <td>8.574830</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>69 rows × 11 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      UTM X (m)    UTM Y (m)      a (px)      b (px)     a (m)     b (m)  \\\n",
       "0   836644.2878  7242472.866   87.024678   55.856531  0.104462  0.067049   \n",
       "1   836644.5051  7242472.819   91.028929   55.045931  0.109269  0.066076   \n",
       "2   836644.6023  7242472.809   46.619508   38.267768  0.055961  0.045936   \n",
       "3   836644.5159  7242472.703  150.841046   69.113499  0.181065  0.082962   \n",
       "4   836644.1702  7242472.628  114.262068  108.564897  0.137157  0.130318   \n",
       "..          ...          ...         ...         ...       ...       ...   \n",
       "64  836644.9420  7242472.026   68.150933   55.956297  0.081806  0.067168   \n",
       "65  836645.0608  7242472.013   49.360133   36.526219  0.059251  0.043845   \n",
       "66  836644.9168  7242471.961   80.284876   42.840646  0.096372  0.051425   \n",
       "67  836645.0932  7242471.847   37.748005   31.099214  0.045312  0.037331   \n",
       "68  836645.2541  7242472.284   34.972655   22.441033  0.041980  0.026938   \n",
       "\n",
       "    area (px)  area (m2)  orientation  ellipse area (px)  perc. diff. area  \n",
       "0        3401   0.004900    -2.309977        3807.285687         10.671269  \n",
       "1        3433   0.004947    -2.527877        3924.675289         12.527795  \n",
       "2        1288   0.001856    -0.218997        1397.332929          7.824401  \n",
       "3        7258   0.010458    -2.214070        8165.475866         11.113570  \n",
       "4        8869   0.012779    -3.002590        9716.068830          8.718226  \n",
       "..        ...        ...          ...                ...               ...  \n",
       "64       2631   0.003791    -1.690621        2986.894275         11.915195  \n",
       "65       1106   0.001594    -2.726250        1412.147675         21.679579  \n",
       "66       2194   0.003161    -2.755353        2693.945657         18.558120  \n",
       "67        819   0.001180    -2.054698         919.480448         10.927959  \n",
       "68        562   0.000810    -0.949405         614.710368          8.574830  \n",
       "\n",
       "[69 rows x 11 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pebbles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d4626b5",
   "metadata": {},
   "source": [
    "We can pull out a column, in this case the b-axes of all pebbles measured."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "794fa8cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_axes_KMS = pebbles['b (m)'].to_numpy(copy=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6137c67",
   "metadata": {},
   "source": [
    "Let's convert that from meters to $\\psi$ units, where $\\psi=log_2(mm)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "60b601ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_axes_KMS = np.log2(b_axes_KMS * 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf7b864",
   "metadata": {},
   "source": [
    "And plot the CDF, or grain-size distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "dcc90e1c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'CDF')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAFBCAYAAAAR7ubGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAh4ElEQVR4nO3deZQU9b3+8fdHwICIooIbi6DBhURAHcElEiNqwA23n+KauISouORG0XiPXnMiGo3J1agoEjfUCCaKiivqBcElLANBFNA4EoQBIjvIooDz+f3x7dGm6Zmega6p6u7ndc6cXqpm5qGEx6ququ/X3B0REanZNnEHEBFJOhWliEgOKkoRkRxUlCIiOagoRURyUFGKiOTQOO4A9dWqVSvv0KFD3DFEpMhMmTJlibu3zras4IqyQ4cOlJeXxx1DRIqMmX1e0zIdeouI5KCiFBHJQUUpIpKDilJEJAcVpYhIDipKEZEcVJQiIjlEVpRm9qiZLTKzj2pYbmZ2r5lVmNl0Mzs4qiwiIlsjyj3Kx4HetSzvA3RKffUHHowwi4jIFovszhx3H29mHWpZpS/whIch1ieYWUsz28PdF0aVSUSKVFUVLFkC8+eHrwUL4OSTYY898vLj47yFsQ0wL+11Zeq9zYrSzPoT9jpp3759g4QTkYRYt+67AkwvwszXGzZs+n1t2xZFUVqW97JO4OPuQ4GhAGVlZZrkR6QYVFXB4sWbl2BmES5fvvn3Nm8ObdqEr6OO+u559deee+atJCHeoqwE2qW9bgssiCmLiOTT2rXZCzC9CBcu3HwvcJttYPfdQ9ntsw/07Ll5CbZpAzvs0KB/nDiLchRwpZmNAHoAK/X5pEiBWrQIbroJ3nsvFOHKlZuv06LFd0V39NHf7fmlF+Buu0Hj5A1qFlkiMxsOHA20MrNK4BagCYC7DwFeBU4AKoC1wEVRZRGRiFRVwaOPwvXXw+rVcMIJcMwxmx4CVz9v0SLutFssyrPe5+RY7sCAqH6/iERsxgz45S/DXmTPnjBkCBxwQNypIqE7c0SkftauhRtvhG7d4OOP4bHH4O23i7YkoQBHOBeRGFRVwbhx8NRT8OyzsGoV/PzncNdd0KpV3Okip6IUkZp9+GEox6efhspK2H57OOMM6N8fjjgi7nQNRkUpIpuqrIThw0NBTp8ezkL37g1//GO422W77eJO2OBUlCISvPZaOJR++21wh8MOg/vvh7POgtZZJycsGSpKEYG33gp7i+3bwy23wHnnwfe/H3eqxFBRipS6mTPhzDOhc2d4990Gv+ulEOjyIJFS9sUXcOKJ0KwZvPyySrIG2qMUKQXuoRQ//RT+9a/w+OmnMGkSLF0K48eHw27JSkUpUkyWLt28DKu/vvzyu/WaNIG99w4XjV9zDZSVxRa5EKgoRQrNypWbl2B1MaYPSbbNNtChA3TqBEceGR6rv/baK5GDTySVtpRIElVVhdsDZ87cvAwXLdp03fbtQ/mdffZ3RbjvvtCxI2y7bTz5i4yKUiQJ5s8PnxdOmgQTJ0J5+aaHynvsEQrw5JNDCVYX4j77hBMxEikVpUhDW7UKpkwJhVhdjvPnh2WNG4fPDS+4ALp3hy5dwvWMBTxEWTFQUYpEacOGcL90+t7irFnhLDSEEvzxj6FHj1CM3bpB06axRpbNqShFonDffTBiBEydCl99Fd5r1SqU4dlnh8dDD4Vddok3p9SJilIk3z77DK6+Gg48EC6/PJRijx7hDLRlm1NPkk5FKZJvI0eGxxdfDGeepeDpFkaRfBs5Eg4+WCVZRFSUIvk0fDhMmACnnx53EskjFaVIPqxeDRddBOeeC4cfDpddFnciySMVpcjWmjIlHGoPGxbmth4/Xmezi4yKUmRLucOf/hT2INetg7Fj4dZbdQ91EVJRimypu++G664L4zlOmxYuHJeipP/1iWyJsWPh+uvhtNPgued0fWSR0x6lSH3Nm/fdSD2PP66SLAEqSpH6cA+zEn71FTz/vKZOKBEqSpH6WLEiXCd5ww2w//5xp5EGoqIUqY9168KjLv8pKSpKkfqoLkoNlltSVJQi9VFRER633z7eHNKgVJQidbViBfTvHwa7OO64uNNIA9J1lCJ14R5KcsECePddne0uMSpKkbp48EH4+9/hjjvCILxSUnToLVKbhQvDdZMDBsDxx8PAgXEnkhioKEWyqaoKe5H77w+jRsGgQfDSS7CN/smUIh16i2T68MPweeSECdCrFwwZEmZLlJKl/z2KpPvb38LYkhUV8OST8OabKknRHqXIJkaOhF13henTdfeNfEt7lCLpli+Hdu1UkrKJSIvSzHqb2SdmVmFmv8myfEcze8nMPjCzGWZ2UZR5RHJatgx23jnuFJIwkRWlmTUCBgN9gM7AOWbWOWO1AcBMd+8KHA38ycy2jSqTSE4rVkDLlnGnkISJco+yO1Dh7rPdfT0wAuibsY4DLczMgO2BZcDGCDOJ1G71at3HLZuJsijbAPPSXlem3kt3P3AAsAD4ELjG3asizCRSuzVroHnzuFNIwkRZlNnGx/eM1z8FpgF7At2A+81ss5tozay/mZWbWfnixYvznVMkmDs3FOV228WdRBImyqKsBNqlvW5L2HNMdxEw0oMK4N/AZsNGu/tQdy9z97LWrVtHFlhK0Ndfh2sne/eGDh3CHTkHHhh3KkmYKItyMtDJzDqmTtD0A0ZlrDMX6AVgZrsB+wGzI8wkEkyfDtdcA3vuGSYKmzkTbr4ZZs+Gfv3iTicJE9kF5+6+0cyuBEYDjYBH3X2GmV2WWj4EuBV43Mw+JByq3+DuS6LKJCVuxQoYPhwefRTKy2HbbeHUU+Hii+HYY6FRo7gTSkJFemeOu78KvJrx3pC05wuA46PMIMI334Th0QYNCrMndukCf/4znHeeLiyXOtEtjFLcFi6E88+HMWPgzDPD7ImHHKK5uKVeVJRSvEaPhgsvhC+/hEcegYsuUkHKFtG93lJ8Nm6EG28MZ7Jbt4bJk8PnkCpJ2UIqSik+zzwTPpO89FKYNAl+8IO4E0mB06G3FJ+JE8PdNQ89pBHJJS/0t0iKz7Rp0LWrSlLyRn+TpLi4wwcfQLducSeRIqKilOIydSqsWqUpZSWvVJRSXF54IRxyn3hi3EmkiKgopbi88AL07Kk7biSvVJRSPObOhY8+gr6Z40OLbB0VpRSPuXPDY+fMGUdEto6KUorHihXhcaedYo0hxUdFKcVj+fLwqMnBJM9UlFI8qvcoVZSSZypKKR4rV4ZHFaXkmYpSisfKldCsGTRpEncSKTIqSikeK1bAjjvGnUKKkIpSisO6dTBqlO7xlkioKKU4PPwwLFoUBuwVyTMVpRS+9evhD3+Ao44Kty+K5JkG7pXC9/zzUFkZ9ipFIqA9Sil8U6aEObp79Yo7iRQpFaUUvhkzYP/9obEOkCQaKkopbF9+GUY01wRiEiEVpRSmWbPgqqugTRuYPx+OPz7uRFLEdKwihWPjRnjpJbj/fhgzJnwuefbZMGAAdO8edzopYipKSb5Fi+AvfwnTz86bB+3bw+23wyWXwK67xp1OSoCKUpLr66/DBeSDB4drJY89Fu67L8yHoxM30oD0t02SqaIiHFZPnQqXXgrXXhvObIvEQEUpyfPMM/CLX4S9xhdfhFNOiTuRlDid9ZbkWLcOfvlL6NcPfvhDmDZNJSmJoKKU5Lj2Whg6FG64AcaNCydtRBJAh96SDAsXwiOPhD3KO+6IO43IJrRHKclwzz3hOsmBA+NOIrIZFaXEb/lyePBBOOss2GefuNOIbEZFKfG79lpYs0aD7kpiqSglXi+/DI89Fk7gdOkSdxqRrFSUEp9ly8L1kgceCLfcEncakRrprLfEY9mycOfNkiXw2mvwve/FnUikRpHuUZpZbzP7xMwqzOw3NaxztJlNM7MZZjYuyjySEBMmwEEHhWslH3pIMydK4kVWlGbWCBgM9AE6A+eYWeeMdVoCDwCnuPsPgP8XVR5JAHe4++4wCVijRvD++3DxxXGnEskpyj3K7kCFu8929/XACKBvxjrnAiPdfS6Auy+KMI/EacUKOOMM+PWv4aSTwmAXZWVxpxKpkyiLsg0wL+11Zeq9dPsCO5nZ22Y2xcwujDCPxGX9+nDP9ksvhT3KkSOhZcu4U4nUWZQncyzLe57l9x8C9AKaAf8wswnu/q9NfpBZf6A/QHvd/1tY3MMI5O+8A08/DeecE3cikXqLco+yEmiX9rotsCDLOq+7+xp3XwKMB7pm/iB3H+ruZe5e1rp168gCSwTuuy/Mt/3f/62SlIIVZVFOBjqZWUcz2xboB4zKWOdF4Cgza2xm2wE9gFkRZpKG9Mor4TPJvn3h1lvjTiOyxSI79Hb3jWZ2JTAaaAQ86u4zzOyy1PIh7j7LzF4HpgNVwMPu/lFUmaQBDR8OF14Y7rZ58knYRvc2SOEy98yPDZOtrKzMy8vL444htbn3XrjmGujZE0aNgh13jDuRSE5mNsXds16Kof/NS/64w003hZI89VQYPVolKUVBRSn5c+utcNttYTKwv/8dmjaNO5FIXuheb8mPiRPhd7+Dc88N0zlYtqvDRAqT9ihl661ZAxdcAG3awAMPqCSl6GiPUrbewIFhHu4xY/SZpBQl7VHK1nn//TCNw3/9Fxx9dNxpRCKhopQtV1UVznC3aRM+nxQpUjr0li03bBiUl8NTT0Hz5nGnEYmM9ihly8yeHSYDO/zwcKZbpIipKKX+3n4buncPw6fpLLeUABWl1M+DD8Jxx0Hr1uHaSU3jICVARSl1s2EDXHFF+Dr++DDvTadOcacSaRC1FqWZPZ72/GeRp5HkuvnmsDd5/fUa6EJKTq49yvRBdK+JMogk2Ndfh8F3Tz8d7rwzTAwmUkJyFWVhjcEm0XjxRVi6FPr3jzuJSCxyXUfZ1szuJcx/U/38W+5+dWTJJDn+8hfYa69wEkekBOUqyoFpzzVabimaOBHeegsGDdIo5VKyai1Kdx/WUEEkgdzDyZtdd4WrdfAgpSvnLoKZ/czMpprZmtRXuebfLhEvvwzjx8NvfwstWsSdRiQ2te5RpgrxV8CvgamEzyoPBu4yM9z9icgTSjzcwy2K++4bRiwXKWG5PqO8AjjN3eekvTfGzM4ARgAqymL10UcwYwYMGQJNmsSdRiRWuQ69d8goSQBS7+0QRSBJiJdeCo+nnBJvDpEEyFWU67ZwmRS6UaPg0ENhjz3iTiISu1yH3geY2fQs7xuwdwR5JAneeAMmTQqzKopIzqLsCuwGzMt4fy9gQSSJJF4ffwxnnQUHHhhGLxeRnIfedwOr3P3z9C9gbWqZFJNly+Dkk+F73wuH3ttvH3cikUTItUfZwd03O/R293Iz6xBNJImFO5x9NsydC2PHhlsWRQTIXZRNa1nWLJ9BJGazZ4dbFe+8E444Iu40IomS69B7spn9IvNNM7sEmBJNJInFZ5+Fx8MOizeHSALl2qP8FfC8mZ3Hd8VYBmwLnBZhLmlo1UW5zz7x5hBJoFyDYnwBHGFmPwF+mHr7FXcfE3kyaVgzZ0LTprpuUiSLOs3r7e5jgbERZ5G4rFkDTz8NffpoKDWRLPSvQuCJJ8KlQddeG3cSkURSUZa6qiq4++4wT7fOdotkVadDbyliU6fCp5+GvUqzuNOIJJL2KEtd9dnugw6KN4dIgqkoS92cOeFRd+KI1EhFWcqWLg23K+68s6Z6EKmFirIUVVTAgAHQrh2MHg0XXxx3IpFE08mcUvKPf8Af/wjPPx+mdzj//HBJUOfOcScTSTQVZbH75pswrcNdd8H778NOO4VJw668UnfhiNRRpIfeZtbbzD4xswoz+00t6x1qZt+Y2ZlR5ilJAwbAaafBwoVw771hGLXbblNJitRDZHuUZtYIGAwcB1QSRiIa5e4zs6x3JzA6qiwla9UqGDYMLrgAHn0UGusAQmRLRLlH2R2ocPfZ7r6eML1t3yzrXQU8ByyKMEtpevZZ+OoruOIKlaTIVoiyKNuw6Vw7lan3vmVmbQjDtQ2JMEdpcg97k506QY8ecacRKWhRFmW2++E84/U9wA3u/k2tP8isv5mVm1n54sWL85WveK1dCxdeCOPHw6WX6tZEka0U5fFYJdAu7XVbNp+5sQwYYeEfcivgBDPb6O4vpK/k7kOBoQBlZWWZZSvp5syB00+HadPgd7+D666LO5FIwYuyKCcDncysIzAf6Aecm76Cu3esfm5mjwMvZ5ak1MNbb0G/frBxY7gk6MQT404kUhQiO/R2943AlYSz2bOAv7n7DDO7zMwui+r3lqxnn4Wf/hR23x0mT1ZJiuRRpKdC3f1V4NWM97KeuHH3n0eZpeg9/ngY2GLCBM3HLZJnute7GLjDpEnQs6dKUiQCKspiMHcuLF4cRikXkbxTURaDSZPC46GHxptDpEipKIvBu+9Cs2bQrVvcSUSKkoqyGLzzDhx2WBg6TUTyTkVZ6FauhA8+gKOOijuJSNFSURayqiq47LLweMIJcacRKVoqykJ2yy0wYgTccYcGvhCJkIqyUA0bBoMGwSWXwPXXx51GpKipKAvRzJnwi1/AMcfAgw9qdCCRiKkoC9Fbb8GGDeG2RZ3pFomcirIQTZsGu+4KbdvGnUSkJKgoC9EHH0DXrjrkFmkgKspCM358KMqysriTiJQMFWUh+fhjOPXUMA/OwIFxpxEpGSrKQvHFF+Gi8iZN4NVXYaed4k4kUjI0h2kh2LABTjkF/vMfGDcOOnbM/T0ikjcqykIwblwYSu3xxzWUmkgMdOhdCN54Ixxyn3FG3ElESpKKshCMHg0/+pGmeRCJiYoy6b74AqZPh+OPjzuJSMlSUSZd9TQPGm9SJDYqyqQrL4dtttE0DyIxUlEmXXk5dO4MzZvHnUSkZKkok27mTOjSJe4UIiVNRZl0K1ZAq1ZxpxApaSrKJHOHVatghx3iTiJS0lSUSfb552HiMBWlSKxUlEn1xhvQvXs4idOrV9xpREqaijJpNm6Em26C3r1ht93CWe+DD447lUhJ06AYSbJgAZx7bhgE45JL4N57Ybvt4k4lUvJUlEnSpw9UVMATT8AFF8SdRkRSdOidFOvXw4cfwnXXqSRFEkZFmRTz54fLgfbaK+4kIpJBRZkUs2aFx3bt4s0hIptRUcbJHcaMCROGnXwytGgBBx4YdyoRyaCijMOaNTB0aCjFXr3gvffgxhvDfd277x53OhHJoLPeDWnOHBg8GB55BJYvh4MOgsceg379oGnTuNOJSA1UlA3BHa6+Gh54AMzC3DdXXQVHHhlei0iiqSgbwj33wP33Q//+cPPN0LZt3IlEpB4i/YzSzHqb2SdmVmFmv8my/Dwzm576et/MukaZJxZvvw0DB8Jpp8GQISpJkQIUWVGaWSNgMNAH6AycY2adM1b7N/Bjd+8C3AoMjSpPLCor4eyzoVOnMCe3DrNFClKUh97dgQp3nw1gZiOAvsDM6hXc/f209ScAxbO7tXw5nHQSrF0b9io1VJpIwYry0LsNMC/tdWXqvZpcArwWYZ6Gs2pVGP1n1ix47jk44IC4E4nIVohyjzLbcaZnXdHsJ4Si/FENy/sD/QHat2+fr3zRWLMm7ElOnQojR2o+bpEiEOUeZSWQfj9eW2BB5kpm1gV4GOjr7kuz/SB3H+ruZe5e1rp160jC5sXnn0PfvuEC8r/+NdxtIyIFL8qinAx0MrOOZrYt0A8Ylb6CmbUHRgIXuPu/IswSrY8/hp//HL7/fRg/PlxEftZZcacSkTyJ7NDb3Tea2ZXAaKAR8Ki7zzCzy1LLhwD/A+wCPGDhjPBGdy+LKlPeTZkCv/99OMRu2hSuvBKuvVaXAIkUGXPP+rFhYpWVlXl5eXl8AdzhnXfg9tth9GjYccdwl83VV0OSPxYQkVqZ2ZSadtR0Z059XXQRDBsGu+4Kd9wBl1+uS39EipyKsj7Gjg0lec014ZC7WbO4E4lIA9Awa3VVVRWmaWjXTiUpUmK0R1lXw4eHayOffFIlKVJitEdZV//7v9ClS5hOVkRKioqyLlavhmnTwsXk22iTiZQa/auvi/Ly8BnlYYfFnUREYqCirItx48Jjjx7x5hCRWKgocxk1CgYNgmOOgV12iTuNiMRARVmbV16BM8+Ebt3CcGkiUpJUlDV5/XU4/fQwpewbb0DLlnEnEpGYqCizWbw4zHHTuTO8+SbstFPciUQkRirKbBYsgK++gptugp13jjuNiMRMRZnNxo3hsUmTeHOISCKoKLNRUYpIGhVlNmvWhMemTePNISKJoKLM5t//Do8dOsQaQ0SSQUWZzezZ0LhxGFJNREqeijKb2bNhr71CWYpIyVNRZjN7NnTsGHcKEUkIFWU2n30G++wTdwoRSQgVZaaVK2HpUhWliHxLRZmp+oy3Dr1FJEVFmemLL8LjHnvEm0NEEkNFmWnJkvDYunW8OUQkMVSUmRYvDo+tWsWbQ0QSQ0WZaeZMaNFC40+KyLdUlOmqquDll6F3b822KCLfUhukmzoVFi6EU06JO4mIJIiKstqyZXDbbWFPsk+fuNOISIKoKL/5Bh56CPbdN8y4ePPNmm1RRDZR2qM+vPsuXHUVTJsGP/4x3HsvdOkSdyoRSZjS3KNcsADOOw+OOipcN/nMMzB2rEpSRLIqvT3Kzz+HI44I93PffDPccAM0bx53KhFJsNIqyiVL4Kc/DVM9TJwIXbvGnUhECkDpFOWaNXDSSTBnTpirWyUpInVUGkXpDv36weTJ8Nxz4bNJEZE6Ko2i/PjjcMfN7bfDqafGnUZECkxpnPWeNCk89u0bbw4RKUilU5QtWsB++8WdREQKUGkU5cSJcMgh0KhR3ElEpABFWpRm1tvMPjGzCjP7TZblZmb3ppZPN7OD8x5iyZIw2MXRR+f9R4tIaYisKM2sETAY6AN0Bs4xs84Zq/UBOqW++gMP5j3Im2+Gs94a6EJEtlCUe5TdgQp3n+3u64ERQObZlL7AEx5MAFqaWX4nq3nttTDIxSGH5PXHikjpiLIo2wDz0l5Xpt6r7zqYWX8zKzez8sXVUzXUhTuMHh3uxtHnkyKyhaK8jtKyvOdbsA7uPhQYClBWVrbZ8poTWLjI/Ouv6/wtIiKZoizKSqBd2uu2wIItWGfrtG+f1x8nIqUnykPvyUAnM+toZtsC/YBRGeuMAi5Mnf0+DFjp7gsjzCQiUm+R7VG6+0YzuxIYDTQCHnX3GWZ2WWr5EOBV4ASgAlgLXBRVHhGRLRXpvd7u/iqhDNPfG5L23IEBUWYQEdlapXFnjojIVlBRiojkoKIUEclBRSkikoOKUkQkBxWliEgOKkoRkRwsXMpYOMxsMfB5lkWtgCUNHKe+CiEjKGc+FUJGKIycUWfcy91bZ1tQcEVZEzMrd/eyuHPUphAygnLmUyFkhMLIGWdGHXqLiOSgohQRyaGYinJo3AHqoBAygnLmUyFkhMLIGVvGovmMUkQkKsW0RykiEomCK0oza2Rm/zSzl7MsO9rMVprZtNTX/8SUcY6ZfZjKUJ5lefTT9OYnZ+zb08xamtmzZvaxmc0ys8MzlidlW+bKmYRtuV/a759mZqvM7FcZ68S6PeuYscG3ZaTjUUbkGmAWsEMNy99x95MaME9NfuLuNV3zlT5Nbw/CNL09GipYhtpyQvzb88/A6+5+Zmqk/O0ylidlW+bKCTFvS3f/BOgG304nPR94PmO1WLdnHTNCA2/LgtqjNLO2wInAw3Fn2UrRT9NbBMxsB6An8AiAu6939xUZq8W+LeuYM2l6AZ+5e+bNG7FvzzQ1ZWxwBVWUwD3A9UBVLescbmYfmNlrZvaDhom1GQfeMLMpZtY/y/I6TdPbAHLlhHi3597AYuCx1MctD5tZ84x1krAt65ITkvF3s1o/YHiW95OwPavVlBEaeFsWTFGa2UnAInefUstqUwm3IXUF7gNeaIhsWRzp7gcTDmMGmFnPjOV1mqa3AeTKGff2bAwcDDzo7gcBa4DfZKyThG1Zl5xxb8tvpT4aOAX4e7bFWd5r8L+bOTI2+LYsmKIEjgROMbM5wAjgGDN7Kn0Fd1/l7qtTz18FmphZq4YO6u4LUo+LCJ+vdM9YJfppeusgV84EbM9KoNLdJ6ZeP0sopMx14t6WOXMmYFum6wNMdfcvsixLwvaEWjLGsS0Lpijd/UZ3b+vuHQi75GPc/fz0dcxsdzOz1PPuhD/f0obMaWbNzaxF9XPgeOCjjNVin6a3Ljnj3p7u/h9gnpntl3qrFzAzY7XYt2Vdcsa9LTOcQ82HtLFvz5QaM8axLQvxrPcmbNPpb88ELjezjcA6oJ83/BX1uwHPp/47NgaedvfXLXnT9NYlZxK251XAX1OHYrOBixK4LeuSMwnbEjPbDjgO+GXae4nannXI2ODbUnfmiIjkUDCH3iIicVFRiojkoKIUEclBRSkikoOKUkQkBxWliEgOKkoRkRxUlLLFzKyDmWXedZSvn/1+PdZdneff3dLMrkh7Xec/p5k1M7NxqSHCalpnWzMbb2YFf8NHqVBRSiK5+xEx/vqWwBW5VqrBxcBId/+mphXcfT3wf8DZW/g7pIGpKGVrNTazYRZGw342dfvZJszshdRQbjPSh3Mzs0NT39c0de/5DDP7YWrZ6tR7r6SG0/rIzGotFjM738wmWRj1+qHqvbrUHuEsM/tL6ne8YWbNUstutjAq+ZtmNtzMrgPuAPZJ/Zy7Uj++Ubbvz+I84MW0TNPMbDczG2RmP7MwOvcIwog359V1I0u8VJSytfYDhrp7F2AV2ffELnb3Q4Ay4Goz2wXA3ScTBmEYBPwBeMrd0w9xewML3L2ru/8QeL2mEGZ2AGEP7Uh37wZ8w6ZF1AkY7O4/AFYAZ5hZGXAGcBBweiofhCHSPnP3bu4+sKbvz5JhW2Bvd5+Tet0Y2Dk1Ak5X4AOgS+rxI+DQmv48kiwqStla89z9vdTzp4AfZVnnajP7AJhAGMKrU9qy3xEGQCgjlGW6D4FjzexOMzvK3VfWkqMXcAgw2cympV7vnbb83+4+LfV8CtAhlfVFd1/n7l8CL9Xy87N9f6ZWhBKttj9h2hKAzoQRhboC01OH5uurR3CSZFNRytbKHFXFzWyAfTfx07nAscDhqYFW/wk0TVt/Z2B7oEXG+7j7vwjl9yHwe6t9EikDhqX2Aru5+37u/tu05V+nPf+GMGJStkFqa5Lt+zOtY9M/w37AJ2a2M7A69dlkGVA9kdv3gK/qkUFioqKUrdXevptx8BzgXXcfXF1YhNG+l7v7WjPbHzgs4/uHAjcDfwXuTF9gZnsCa939KeCPbD5ob7r/A840s11T37uzme2VI/u7wMmpz0i3J8zHBPAlobjrxd2XEz7LrC7L9YS9yjLgAzM7H5jj7l+kPn5Y7O4b6vt7pOHp8gTZWrOAn5nZQ8CnhFn70r0OXGZm04FPCIffAJjZhcBGd386deLlfTM7xt3HpFY5ELjLzKqADcDlNYVw95lmdhNhDqBtUusPAGqcmMrdJ5vZKMJnhp8T9vRWuvtSM3svdUnQa8DgOm8NeINwSP9W6s9+HOF/AuuAZcCFqfV+Qhj7UQqAxqOUkmZm27v76tTZ+vFAf3efuhU/7yDg1+5+Qdp7jxEGRn4z7b2RwI0epmeVhNOht5S6oamTP1OB57amJAHc/Z/A2IwLzrsA06tfpM6Ov6CSLBzaoxQRyUF7lCIiOagoRURyUFGKiOSgohQRyUFFKSKSg4pSRCQHFaWISA4qShGRHP4/li2Ak9/zqwkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(sorted(b_axes_KMS), np.linspace(0, 1, len(b_axes_KMS)), color='r')\n",
    "plt.xlabel('b-axis length ($\\psi$)')\n",
    "plt.ylabel('CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65c6631c",
   "metadata": {},
   "source": [
    "# PebbleCountsAuto (automatic with image filtering; AIF)\n",
    "\n",
    "Now we will repeat all of that for a `PebbleCountsAuto.py` (AIF algorithm) .csv file."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaa154d3",
   "metadata": {},
   "source": [
    "## Select CSV\n",
    "\n",
    "**Note:** If you don't have the file below in the `example_data` folder, then you need to run `PebbleCountsAuto.py` to produce it. Or just follow along with the notebook without it and copy-paste the code you need into your own projects. This file was generated with the following command (assuming you have already clicked and generated the sand mask .tif file in the previous step):\n",
    "\n",
    "```\n",
    "python PebbleCountsAuto.py -im example_data\\ortho_resolution_1.2mmPerPix.tif -ortho y -sand_mask example_data\\ortho_resolution_1.2mmPerPix_PebbleCounts_SandMask_TIFF.tif -otsu_th 50\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d6236f36",
   "metadata": {},
   "outputs": [],
   "source": [
    "fn = 'example_data/ortho_resolution_1.2mmPerPix_PebbleCountsAuto_CSV.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e32e88ca",
   "metadata": {},
   "source": [
    "## Read PebbleCounts Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b95f212f",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = pd.read_csv(fn, skiprows=1, nrows=1).T.to_dict()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "521f87d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'otsu_threshold': 45.0,\n",
       " 'cutoff': 20.0,\n",
       " 'percent_overlap': 15.0,\n",
       " 'misfit_threshold': 30.0,\n",
       " 'min_size_threshold': 10.0,\n",
       " 'first_nl_denoise': 5.0,\n",
       " 'tophat_th': 85.0,\n",
       " 'sobel_th': 85.0,\n",
       " 'canny_sig': 2.0}"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69c1e900",
   "metadata": {},
   "source": [
    "There are no list-type entries in the case of `PebbleCountsAuto.py` parameters, so we don't need to worry about the `literal_eval` function as we did above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0ff7925",
   "metadata": {},
   "source": [
    "## Read Image Details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "665fae56",
   "metadata": {},
   "outputs": [],
   "source": [
    "tile_details = pd.read_csv(fn, skiprows=5, nrows=1).T.to_dict()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8e361e30",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'perc. not meas.': 50.78, 'perc. background color': 22.8574}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tile_details"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2af1e48",
   "metadata": {},
   "source": [
    "## Read and Plot Pebbles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "bc94db40",
   "metadata": {},
   "outputs": [],
   "source": [
    "pebbles = pd.read_csv(fn, skiprows=9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "3bb3b010",
   "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>UTM X (m)</th>\n",
       "      <th>UTM Y (m)</th>\n",
       "      <th>a (px)</th>\n",
       "      <th>b (px)</th>\n",
       "      <th>a (m)</th>\n",
       "      <th>b (m)</th>\n",
       "      <th>area (px)</th>\n",
       "      <th>area (m2)</th>\n",
       "      <th>orientation</th>\n",
       "      <th>ellipse area (px)</th>\n",
       "      <th>perc. diff. area</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>836644.287800</td>\n",
       "      <td>7.242473e+06</td>\n",
       "      <td>91.099302</td>\n",
       "      <td>56.929474</td>\n",
       "      <td>0.109353</td>\n",
       "      <td>0.068337</td>\n",
       "      <td>3585</td>\n",
       "      <td>0.005166</td>\n",
       "      <td>-2.346842</td>\n",
       "      <td>4062.106433</td>\n",
       "      <td>11.745296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>836644.664717</td>\n",
       "      <td>7.242473e+06</td>\n",
       "      <td>87.562741</td>\n",
       "      <td>73.385072</td>\n",
       "      <td>0.105108</td>\n",
       "      <td>0.088089</td>\n",
       "      <td>4281</td>\n",
       "      <td>0.006168</td>\n",
       "      <td>-2.740125</td>\n",
       "      <td>5032.990990</td>\n",
       "      <td>14.941235</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>836644.431844</td>\n",
       "      <td>7.242473e+06</td>\n",
       "      <td>44.341632</td>\n",
       "      <td>23.999295</td>\n",
       "      <td>0.053226</td>\n",
       "      <td>0.028808</td>\n",
       "      <td>678</td>\n",
       "      <td>0.000977</td>\n",
       "      <td>-1.838869</td>\n",
       "      <td>833.506962</td>\n",
       "      <td>18.656948</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>836644.506268</td>\n",
       "      <td>7.242473e+06</td>\n",
       "      <td>92.439524</td>\n",
       "      <td>58.495033</td>\n",
       "      <td>0.110962</td>\n",
       "      <td>0.070216</td>\n",
       "      <td>3794</td>\n",
       "      <td>0.005467</td>\n",
       "      <td>-2.370185</td>\n",
       "      <td>4235.218002</td>\n",
       "      <td>10.417835</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>836644.113746</td>\n",
       "      <td>7.242473e+06</td>\n",
       "      <td>33.943588</td>\n",
       "      <td>24.522035</td>\n",
       "      <td>0.040745</td>\n",
       "      <td>0.029436</td>\n",
       "      <td>505</td>\n",
       "      <td>0.000728</td>\n",
       "      <td>-0.035510</td>\n",
       "      <td>651.948572</td>\n",
       "      <td>22.539902</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",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75</th>\n",
       "      <td>836645.099251</td>\n",
       "      <td>7.242472e+06</td>\n",
       "      <td>46.812721</td>\n",
       "      <td>32.803184</td>\n",
       "      <td>0.056193</td>\n",
       "      <td>0.039376</td>\n",
       "      <td>1090</td>\n",
       "      <td>0.001571</td>\n",
       "      <td>-2.607045</td>\n",
       "      <td>1202.759965</td>\n",
       "      <td>9.375101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>76</th>\n",
       "      <td>836645.214487</td>\n",
       "      <td>7.242472e+06</td>\n",
       "      <td>74.834690</td>\n",
       "      <td>36.298176</td>\n",
       "      <td>0.089829</td>\n",
       "      <td>0.043571</td>\n",
       "      <td>1910</td>\n",
       "      <td>0.002752</td>\n",
       "      <td>-0.057553</td>\n",
       "      <td>2127.584663</td>\n",
       "      <td>10.226839</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>77</th>\n",
       "      <td>836645.141264</td>\n",
       "      <td>7.242472e+06</td>\n",
       "      <td>64.393964</td>\n",
       "      <td>47.678347</td>\n",
       "      <td>0.077297</td>\n",
       "      <td>0.057232</td>\n",
       "      <td>2070</td>\n",
       "      <td>0.002983</td>\n",
       "      <td>-0.485613</td>\n",
       "      <td>2404.725066</td>\n",
       "      <td>13.919473</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>78</th>\n",
       "      <td>836644.123349</td>\n",
       "      <td>7.242472e+06</td>\n",
       "      <td>50.170835</td>\n",
       "      <td>25.534225</td>\n",
       "      <td>0.060224</td>\n",
       "      <td>0.030651</td>\n",
       "      <td>905</td>\n",
       "      <td>0.001304</td>\n",
       "      <td>-2.647203</td>\n",
       "      <td>1003.397672</td>\n",
       "      <td>9.806448</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>79</th>\n",
       "      <td>836645.182077</td>\n",
       "      <td>7.242472e+06</td>\n",
       "      <td>26.515197</td>\n",
       "      <td>24.846945</td>\n",
       "      <td>0.031828</td>\n",
       "      <td>0.029826</td>\n",
       "      <td>469</td>\n",
       "      <td>0.000676</td>\n",
       "      <td>-2.068539</td>\n",
       "      <td>516.020479</td>\n",
       "      <td>9.112134</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>80 rows × 11 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        UTM X (m)     UTM Y (m)     a (px)     b (px)     a (m)     b (m)  \\\n",
       "0   836644.287800  7.242473e+06  91.099302  56.929474  0.109353  0.068337   \n",
       "1   836644.664717  7.242473e+06  87.562741  73.385072  0.105108  0.088089   \n",
       "2   836644.431844  7.242473e+06  44.341632  23.999295  0.053226  0.028808   \n",
       "3   836644.506268  7.242473e+06  92.439524  58.495033  0.110962  0.070216   \n",
       "4   836644.113746  7.242473e+06  33.943588  24.522035  0.040745  0.029436   \n",
       "..            ...           ...        ...        ...       ...       ...   \n",
       "75  836645.099251  7.242472e+06  46.812721  32.803184  0.056193  0.039376   \n",
       "76  836645.214487  7.242472e+06  74.834690  36.298176  0.089829  0.043571   \n",
       "77  836645.141264  7.242472e+06  64.393964  47.678347  0.077297  0.057232   \n",
       "78  836644.123349  7.242472e+06  50.170835  25.534225  0.060224  0.030651   \n",
       "79  836645.182077  7.242472e+06  26.515197  24.846945  0.031828  0.029826   \n",
       "\n",
       "    area (px)  area (m2)  orientation  ellipse area (px)  perc. diff. area  \n",
       "0        3585   0.005166    -2.346842        4062.106433         11.745296  \n",
       "1        4281   0.006168    -2.740125        5032.990990         14.941235  \n",
       "2         678   0.000977    -1.838869         833.506962         18.656948  \n",
       "3        3794   0.005467    -2.370185        4235.218002         10.417835  \n",
       "4         505   0.000728    -0.035510         651.948572         22.539902  \n",
       "..        ...        ...          ...                ...               ...  \n",
       "75       1090   0.001571    -2.607045        1202.759965          9.375101  \n",
       "76       1910   0.002752    -0.057553        2127.584663         10.226839  \n",
       "77       2070   0.002983    -0.485613        2404.725066         13.919473  \n",
       "78        905   0.001304    -2.647203        1003.397672          9.806448  \n",
       "79        469   0.000676    -2.068539         516.020479          9.112134  \n",
       "\n",
       "[80 rows x 11 columns]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pebbles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "408de681",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_axes_AIF = pebbles['b (m)'].to_numpy(copy=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "e1534b5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_axes_AIF = np.log2(b_axes_AIF * 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "00b2673f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x249301739d0>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAFBCAYAAAAR7ubGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3FElEQVR4nO3dd3wU1frH8c9D6NICJIgECAiEDkroIiiogBQpUpRuFysKKJbr1XvFC4KKckVQQBREQUHwoiAKUhQhdOkI8UeooZcAScj5/XGWuAkbNkA2s5t93q/XvrI7M7v7ZAxfz8yZOUeMMSillMpYLqcLUEopf6dBqZRSXmhQKqWUFxqUSinlhQalUkp5oUGplFJe5Ha6gCtVsmRJExkZ6XQZSqkcZvXq1YeNMWGe1gVcUEZGRhITE+N0GUqpHEZE/sponR56K6WUFxqUSinlhQalUkp5EXDnKD1JSkoiLi6Oc+fOOV2K8kP58+cnIiKCPHnyOF2KClA5Iijj4uIoXLgwkZGRiIjT5Sg/YozhyJEjxMXFUaFCBafLUQEqRxx6nzt3jhIlSmhIqkuICCVKlNCjDXVNckRQAhqSKkP6t6Gulc+CUkQmisghEfkjg/UiImNEZKeIbBCRm31Vi1JKXQtftignA60vs74NUNn1eBj40Ie1BI3JkyfzxBNPXHZ5SkoKffv2ZcCAARhjiIyMpFmzZmm2r1u3LjVr1gQgISGB+++/n1q1alGzZk1uueUWTp8+fcl3GGO4/fbbOXny5DX/Hhs2bKBx48bUqFGDWrVqpR46f/nll9SuXZsaNWowZMiQ1O0/+OADJk2adM3fq5QnPuvMMcYsEZHIy2zSEZhi7BDrK0SkmIiUNsbs91VNyobZo48+SlJSEpMmTUo9LD116hR79uyhbNmybNmyJc173nvvPUqVKsXGjRsB2LZtm8ce5Hnz5lGnTh2KFClyTTUmJyfTq1cvPvvsM+rUqcORI0fIkycPR44cYfDgwaxevZqwsDD69u3LTz/9RMuWLRkwYABNmzalf//+1/Tdys+lpMDBg7B3L8TFwd69mAce5MS5fBw6BPHxpP5s1w5uuCFrvtbJc5RlgD1ur+Ncyy4hIg+LSIyIxMTHx2dLcVciNjaWqlWr8uCDD1KzZk3uv/9+Fi5cSNOmTalcuTIrV64EYOXKlTRp0oSbbrqJJk2asG3bNsC29jp37kzr1q2pXLlympZSoUKFUp/PnDmTfv36ATB37lwaNmzITTfdRKtWrTh48GCman366ac5cuQIU6ZMIVeuv//zd+vWjS+//BKAL774gp49e6au279/P2XK/P2fJioqinz58l3y2VOnTqVjx46p+6RatWo89NBD1KhRgzvvvJOzZ89mqsYFCxZQu3Zt6tSpA0CJEiUICQlh165dVKlShbAweztuq1at+PrrrwEoWLAgkZGRqftaBaDERIiNxSxdxunJM9k17GNW9HiXubf8h09ufJO3ir3FoNxj6HXDT9xV/wg3dSpPmSfuIV+RvISGQlQU3HILdO4MjzwCGzZkXWlOXh7k6Qy7xwl8jDHjgfEA0dHRl5/k55lnYN26aywtnbp14d13L7vJzp07mTFjBuPHj6d+/fpMmzaNZcuWMWfOHN58801mz55N1apVWbJkCblz52bhwoUMGzYs9R/6unXrWLt2Lfny5SMqKoonn3ySsmXLZvh9t9xyCytWrEBE+PjjjxkxYgSjRo26bI3Tpk2jWrVqLF68mNy50/6n79q1K/369eP5559n7ty5TJ06lc8++wyAAQMGcOeddzJz5kxatmxJ3759qVy58iWfv3z5cj766KPU1zt27OCLL75gwoQJdOvWja+//ppevXoxcuRIpk6desn7b731VsaMGcP27dsREe666y7i4+Pp0aMHQ4YMoVKlSmzdupXY2FgiIiKYPXs2iYmJqe+Pjo5m6dKlNGjQ4LL7QWW/s4dOEf/HQQ5tOUL8rlMcik3g0L5k2wI8lptDpwsSn1iUQ4RziHqco4DHz7kuz3nCi5wjrHgKEdcLN0fkJawshJeCsDAID7ePsDAoVSrr6ncyKOMA9ySIAPY5VMs1q1ChArVq1QKgRo0atGzZEhGhVq1axMbGAnDixAn69u3Ljh07EBGSkpJS39+yZUuKFi0KQPXq1fnrr78uG5RxcXF0796d/fv3k5iYmKlrBG+++Wa2bt3KypUradq0aZp1xYsXJzQ0lOnTp1OtWjUKFiyYuq5u3brs2rWLBQsWsHDhQurXr89vv/1GtWrV0nzG0aNHKVy4cJp9UrduXQDq1auXuh8GDx7M4MGDM6wzOTmZZcuWsWrVKgoWLEjLli2pV68eLVu25MMPP6R79+7kypWLJk2asGvXrtT3hYeHs3XrVq/7QWUhY1jxw3H+XHuSQ7tOcygukfiDFzh0JIRDJ/MTn3Adh5JCOU1hoDBQKc3b88l5wvOdILzIWcKKJVE9PJHw0gcIK1+A8BuLEFauAOGlJDX8ChbMB1x6NONrTgblHOAJEZkONAROZMn5SS8tP19xPxTNlStX6utcuXKRnJwMwCuvvMJtt93GrFmziI2NpUWLFh7fHxISkvoe90tb3K8FfPLJJxk0aBAdOnRg8eLFvPbaa15rrFq1Kq+//jrdunVj/vz51KhRI8367t27M3DgQCZPnnzJewsVKkTnzp3p3LkzuXLlYt68eZcEZe7cuUlJSUk9pE//O1089PbWooyIiKB58+aULFkSgLZt27JmzRpatmxJ+/btad++PQDjx48nJCQkzf4pUMBzS0RloZUr4dVXYft2Pt5zFw8lfwiEApCbJMKIJzzPMcILnOTG6+MJL7GLsFK5CC+bj/AK1xFWuRjh1UsSXi4/hQrlQyTc2d8nE3wWlCLyBdACKCkiccA/gDwAxphxwDygLbATSABy/Fn4EydOpJ7r8xRGnpQqVYotW7YQFRXFrFmzUlts7p/16aefZrqGJk2aMG7cOO6++26WLFlCuXLlUtd16tSJ/fv3c9ddd7Fv39+N++XLl1O9enVCQ0NJTExk8+bNaUL+oqioKHbt2kWlSpUuWefOW4vyrrvuYsSIESQkJJA3b15++eUXnn32WQAOHTpEeHg4x44d47///S9fffVV6vu2b99+SUtZZaHTp+Hll2HMGLj+evY06MKgvSNpUSmOcU9tJrxqcYpVK41cXwpyZ1Evip/wZa93Ty/rDTDQV9/vj4YMGULfvn0ZPXo0t99+e6be89Zbb9GuXTvKli1LzZo1Uy/Lee2117j33nspU6YMjRo1Yvfu3Zmuo127dsTHx9O6dWuWLl2aurxw4cIMHTr0ku3//PNPHnvsMYwxpKSkcPfdd9OlS5dLtrv77rtZvHix16D0JjQ0lEGDBlG/fn1EhLZt23L33XcDtjNq/fr1ALz66qtUqVIl9X3Lly/nH//4xzV9t8rAvHnw2GOwZw889hjmzeE80rMIF3LDJ/MjqFgxwukKfcsYE1CPevXqmfQ2b958yTKV/fbt22datWrlyHevWbPG9OrVK8P1+jdylQ4cMKZHD2PAmOrVjVm+3BhjzOzZdtF77zlcXxYCYkwGuZNjbmFUzitdujQPPfRQllxwfqUOHz7MG2+8ke3fmyMlJcF330G3blC+PHzzDfzzn7BmDTRpAsCIEVChAgwMkmPCHDF6kPIf3bp1c+R777jjDke+N8cwxl5W9+mnMG2avWK7ZEl7QeLjj9uLFF1+/dU+3n8f3PrScjQNSqWC2b59MHUqTJkCf/wBefNChw7Qpw+0bg3p7sBatMjmZvHiEEw3Qemht1LB6Jdf4K67oGxZGDIECheGDz+E/fthxgxo3z5NSO7YAR07wu23w9mzttF53XUO1p/NtEWpVLBZvRratIESJeCll6B3b/BwpxXAsWPw+uvwwQeQPz+8+SY8+6x9Hkw0KJUKJvv22UPr8HB74Xi454u9N22CiRNh0iQ4fhweeADeeAOuvz57y/UXeujtkNmzZ7N58+arem+dOnXSDFoB0K9fP2bOnAlAixYtiIqKom7dutStWzd1uTuTRUOiJSUl0bdvX2rVqkW1atUYPnx46rrVq1dTq1YtKlWqxFNPPYW9AkOHRHNMQoI9fj5xAubMuSQkT5yA8eOhYUOoWdNeV37bbbaze8KE4A1J0BalY2bPnk27du2oXr36Fb1vy5YtpKSksGTJEs6cOcN1GZwomjp1KtHR0Rl+TlYNiTZjxgzOnz/Pxo0bSUhIoHr16vTs2ZPIyEgee+wxxo8fT6NGjWjbti0//PADbdq00SHRssPJk/bE4vbt9rFtmz3k3rEDZs+G2rUB29m9ZIltPc6YYc8/1qgBo0dDr172/mqlLcosc88991CvXj1q1KjB+PHjU5d7Gibt119/Zc6cOQwePJi6devy559/sm7dOho1akTt2rXp1KkTx44d8/g906ZNo3fv3tx5553MmTPnquvNqiHRRIQzZ86QnJzM2bNnyZs3L0WKFGH//v2cPHmSxo0bIyL06dOH2bNnAzokWpZJTIStW23r8O234eGHoXlzKF0aihaF6Gi47z57DeRvv0FkpO3d7tABgC1boEoVaNHCZmefPvZofONGex5SQ/JvOa5F6dAoa0ycOJHixYtz9uxZ6tevT5cuXShRooTHbZs0aUKHDh1o164dXbt2BaB27dq8//77NG/enFdffZV//vOfvOvhS7/88kt+/PFHtm3bxgcffHDJIfhF999/f+oAET/99NMltWTVkGhdu3bl22+/pXTp0iQkJPDOO+9QvHhxYmJiiIj4+7a2iIgI9u7dm/pah0TLpJQUO0jtxZahewtx9267/qKwMJt8bdrY6x6rVLGPG2/02Pvy7ru2k/uzz+wYjm4DRql0clxQOmXMmDHMmjULgD179rBjx44MgzK9EydOcPz4cZo3bw5A3759uffeey/ZbtWqVYSFhVG+fHkiIiIYMGAAx44dIzQ09JJtvR16Z9WQaCtXriQkJIR9+/Zx7NgxmjVrRqtWrVLPR7pzHwlJh0RL5/x5WL/ethDdA3HHDntu8aKCBW341asHPXva51FRttfaw99BRpKT4euvbeOyVy8f/D45TI4LSidGWVu8eDELFy7kt99+o2DBgrRo0SJ1SLSMhkm7Gl988QVbt24lMjISgJMnT/L111/z4IMPXvFnZdWQaNOmTaN169bkyZOH8PBwmjZtSkxMDM2aNSMuLi51+7i4OG5wG5c/qIdEMwZiY2HFCvv4/XdYu9YeSoO93aViRRuCt9/+d8swKsrObXAVs0oaA0eO/H3a8rff7GuHbqQKODkuKJ1w4sQJQkNDKViwIFu3bmXFihWp6zIaJq1w4cKcOnUKgKJFixIaGsrSpUtp1qwZn332WWrr8qKUlBRmzJjBhg0bUodXW7RoEf/617+uKiizaki0cuXK8fPPP9OrVy8SEhJYsWIFzzzzDKVLl6Zw4cKsWLGChg0bMmXKFJ588snU9wXVkGinTsGqVWmD8dAhu65AAahf354zatDAdjdXrHjJHTGZdfq0DUP3fpyLD/fT3rlz297t1peb/k+l0qDMAq1bt2bcuHHUrl2bqKgoGjVqlLouo2HSevTowUMPPcSYMWOYOXMmn376KY8++igJCQlUrFjxkstnlixZQpkyZdLMXXPrrbeyefNm9u+/8vGOs2pItIEDB9K/f39q1qyJMYb+/ftT29Wj+uGHH9KvXz/Onj1LmzZtaNOmTer7cuyQaBcu2F6S33//Oxg3bbJNOoCqVe05xEaN7KNmTZtaVyAxEXbtSnt0fvG52zCigL3xpnJl6N7974ZplSq2X+cqszgoiadzSf4sOjraxMTEpFm2ZcuWS0bbVpe3f/9++vTpw48//pjt37127VpGjx6dOidPdvD538isWTB2rO02dh0pEBr6dyA2bGhbjJk8j3jypD06v/hwD8b0fTglS9owdA/CKlWgUiXtoLkSIrLaGOPxxL62KIOU+5Bo13ot5ZXKcUOiJSbaESJCQ+3tgBeDsXLlDM8npg/C9I/0V4ddd539OPc+nCpV7LLixX35yynQoAxqOiRaFlm82N7W8vnndjJpXEG4MfNBWKCAPRyOjLQ5e/H5xUdY2FX14agskmOC0hiTpodZqYt8eXppzx6Y9fpxYnOPIXZCG2Jf1SDMiXJEUObPn58jR45QokQJDUuVhjGGI0eOkN9Hw9081/8IM5Z3o2DIOSJ3hhAZCY0bXxqEJUtqEAayHBGUERERxMXFER8f73Qpyg/lz58/zV1CWcIYeOcdVv7Uma4FVvDV7+WRWjWz9juU38gRQZknTx4qVKjgdBkqWBw5Av36cfi73/iLQQwcGo7U0u7lnEwHxVDqSixfbm/+X7CARf3t5U31m2tI5nQalEpl1ubNdvqEfPngt9/46P/aUK4cNGvmdGHK1zQolcqMU6egSxd7QeMvv7C90M389JMd2SxYZiIMZjniHKVSPmUMPPigvS1m4ULMDWV4ubu98/CBB5wuTmUHDUqlvJkwAb76CoYPh9tu4+2RdjTw4cODe3qEYKKH3kp588MPdkSfIUOYPx9eeAHuvReGDnW6MJVdNCiV8ubcOShenBUrc9Gjh51TZtIkvYA8mGhQKuXN2bN8f7YFLVvaqbDnzrV9Oip4aFAqdTmJiUzdVJcOm4cTFWUvoyxf3umiVHbToFTqMj7vMote8e/QrMYxFi+GUqWcrkg5QYNSqQz831crGPhda269fhvzVoWRzcN2Kj+iQamUB+bgIR7sk8gFyc3kHyM8zfaqgogGpVLujMFMncZbFcfz4/lbGTE4ngo1tecm2GlQKnXRX39x+q4u3NdLGJbwMl3uOMGjwyOdrkr5AQ1KpS5cgDFj2FKtMw0WvslX0p1/v5HCVz8UJZf+C1HoLYwq2BkDbdvy9YJC9A1ZRsHieVnwZS5atnS6MOVP9P+XKqiZ4ycYtaAmXfmaWvXzs3Z9iIakuoQGpQpaFy7A088KzzOKrtG7+flnoUwZp6tS/sinQSkirUVkm4jsFJEXPKwvKiJzRWS9iGwSkf6+rEepi/bsga5d4f1Pi/Iso/nyxfUUKOB0Vcpf+SwoRSQEGAu0AaoDPUWkerrNBgKbjTF1gBbAKBHJ66ualPr9d+jRAypUsPdsvztwB6N5jlzFizldmvJjvmxRNgB2GmN2GWMSgelAx3TbGKCw2DlmCwFHgWQf1qSCUFISfPmlnUa2USP4/nt45hn48094+q6tdqNChRytUfk3X/Z6lwH2uL2OAxqm2+YDYA6wDygMdDfGpPiwJhVEkpPhnXdgzBiIi4NKleD996FvXyhc2LXRigT7U4cDUpfhyxalp9H6TLrXdwHrgBuAusAHInLJHbUi8rCIxIhIjM7drTLrk09gyBCoXBnmzIFt2+CJJ9xC0hjYscM+L6gzKaqM+bJFGQeUdXsdgW05uusPvGWMMcBOEdkNVAVWum9kjBkPjAeIjo5OH7ZKeTRzJkRFwU8/pRtk9+BBmDwZPv4Ydu60QwKFhTlVpgoAvmxRrgIqi0gFVwdND+xhtrv/A1oCiEgpIArY5cOaVJA4ehQWLYLOnV0hmZIC8+fbru6ICDufQ+nSMGUK7N6tLUp1WT5rURpjkkXkCWA+EAJMNMZsEpFHXevHAW8Ak0VkI/ZQfagx5rCvalLBY8kSe51ku4bx8MY4exz+119QsiQ8/bSdVbFqVafLVAHCp7cwGmPmAfPSLRvn9nwfcKcva1DB6eyJRCAvxTs1B7MFWrWCESOgY0fIl8/p8lSA0Xu9Vc6zfTvm1UnAcHLd3xNeuw9uvNHpqlQA01sYVc7y1VcQHc3ZI2cACHntFQ1Jdc00KFXOkJgITz4J3btDjRosb/smoaEQGel0YSon0KBUOcOkSfDBB/Dss5jFv7Dg10K0agUhIU4XpnICDUqVM6xcaa+FHDWK7bF52bsX7rjD6aJUTqFBqXKGDRugdm0Q4a+/7KJq1ZwtSeUcGpQq8F24AH/8AXXqAHDYdSWu3myjsooGpQp8K1fCuXPQoAEAF4cDKFnSwZpUjqJBqQLfnDmQOze0bg3Avn2QNy+Ehjpcl8oxNChV4JszB5o3h6JFATt6eUQEOoOiyjL6p6QC219/webN0L596qI9e6Bs2cu8R6krpEGpAltcnP3p1sW9cydUrOhQPSpH0qBUge3YMfvTdULy2DE4cEAvDVJZS4NSBbbjx+3PYsUA2LLFvqyefho7pa6BBqUKbBeD0tWRM3asHUWtXj3nSlI5jwalCmwnTtifRYuydClMmwaDB8P11ztblspZNChVYDt5EvLl40LufDz5pO3tfvFFp4tSOY0O3KsC2/HjULQomzbB+vUwfrxOf6OynrYoVeA6fx7mzoU6dVLv765c2dmSVM6kLUoVuKZMgf37YcoUjh61i0qUcLYklTNpUKrAdOGCnSwsOpqzTVryYXs7SO8NNzhdmMqJNChVYJo/H3bu5PwX39Clq7BoEXz6qbYolW9oUKrAtGYNSeSm++ft+f57mDABevd2uiiVU2lnjgpIm347QdO8q/j2f7n54AN48EGnK1I5mbYoVUBJToaRbybx2vf/pkies8yYAV27Ol2Vyuk0KFXA2PzTfvr1TmbV/rJ0ZQZjXzpCeNdHnS5LBQE99Fb+zRjM0mWMveljbm4Vyu79+fmywShmLA4n/JVHnK5OBQkNSuWfzp2DyZM5Xqc5XW89yBPrHuT2irFsijlHt9+fsyOaizhdpQoSeuit/M/nn8Ozz7Li8I30yPM1e3OVZuQbiQx6oapO76AcoX92yn+cOwePPAK9e7O8dFeahfwKN9zAsl9z8fywvBqSyjHaolT+Ydcu2329di0MHco7O96k6L5crF2rsykq52lQKud9+y307WvPOc6Zw4H67fm2LDz9tIak8g8alMpZa9ZAp05w880wYwbnSlfgxcfs9ZIPPeR0cUpZGpTKWSNGQOHC8NNPrNhSlP5tYetWGDQIoqKcLk4pS0+PK+fs3g0zZnB2wECef6MoTZvCmTN2vItRo5wuTqm/aYtSOWfUKLblqkaHOf9k+y7b4T1iBBQp4nRhSqWlQamcsW4dceO+484Cqzl7Kg8LF0LLlk4XpZRnGpQq+yUlcaT3M9wl8zlGcRb/YPtylPJXeo5SZbtDL4+h3R/D2Uklvp0jGpLK72mLUmUbY+Cr535n4Dt9OZWrKNO/CuG225yuSinvfNqiFJHWIrJNRHaKyAsZbNNCRNaJyCYR+cWX9SjnHIxLomvUBnq805CK1x1k7dIzdOrkdFVKZY7PglJEQoCxQBugOtBTRKqn26YY8F+ggzGmBnCvr+pRzoiPh49GnKB6hQT+t6MK/2kym18PVaZ6k2JOl6ZUpvny0LsBsNMYswtARKYDHYHNbtvcB3xjjPk/AGPMIR/Wo7JBSgqsWwf/+599rFxpMKYoDXOtZNJ/DlJtyD1Ol6jUFfNlUJYB9ri9jgMaptumCpBHRBYDhYH3jDFTfFiT8oFTp+DHH20wzpsHBw7Y5fVrnOEfIaO5+4Z13Dz3n+Sq3d7ZQpW6Sr4MSk+jqhoP318PaAkUAH4TkRXGmO1pPkjkYeBhgHLlyvmgVHWljh6FyZNtOC5dCklJ9kLxu+6Cu++G1tGHKdW2HpQ2sCoGwsOdLlmpq+bLoIwDyrq9jgD2edjmsDHmDHBGRJYAdYA0QWmMGQ+MB4iOjk4ftiqbGQP33GMDsnp1eOYZG45NmkCePNjUvPNeOHQIli3TkFQBz5dBuQqoLCIVgL1AD+w5SXffAh+ISG4gL/bQ/B0f1qSywJQpNiQ/+ggefjjdSmPs+GiLF9sN69VzokSlspTPgtIYkywiTwDzgRBgojFmk4g86lo/zhizRUR+ADYAKcDHxpg/fFWTunZHj8LgwdC4sYe5tC9cgEcfhY8/thv17u1IjUplNZ9ecG6MmQfMS7dsXLrXI4GRvqxDZZ3Bg+HIEViwgLRTMyQm2mD86it46SV44w3HalQqq+mdOSrTvvsOJk6EYcOgbl23FQkJ0KUL/PADjBwJzz/vVIlK+YQGpcqUI0fsiOO1a8Orr7qtMAa6dbODSE6Y4OF4XKnAp0GpvDIGHnvMhuUPP0C+fG4rJ0yw1wi9956GpMqxNCiVV59/DjNmwL//DXXquK3YvRueew5uvx2eeMKx+pTyNR1mTV1WbKzNwFtugaFD3VakpED//nbmxIkT0Um3VU6mLUp1WU8/bQ+9P/sMQkLcVnz2Gfzyiz30Ll/esfqUyg7aDFCXtWkTtGsHkZFuC8+csV3fDRrAgAFOlaZUttEWpbqshAS47rp0C0eOhH377DWTesitgoD+lasMpaTYkYHSBOWGDXaqxG7doGlTx2pTKjtpUKoMTZ0Kp0/bwS4Ae8V506YQGmrDUqkgoUGpPDp3zt6JWK8edO1iYPRo6NABoqJg5UrtwFFBRc9RKo/eeQf27IFPP0km12OP297tLl3siEAFCzpdnlLZ6rItShGZ7Pa8r8+rUX5h6VJ7m2LnznDbz6/YkBw2zHbeaEiqIOTt0Nv9PoynfVmI8g9xcdC1K1SoAJ98mAiffAKdOtnbcrSHWwUpb4feOpp4EDl82LYiExJg0SIotvx/dhpFvYdbBTlvQRkhImOw899cfJ7KGPOUzypT2cYYe1T95JNw7BjMnGmneGDoRLjhBrjzTqdLVMpR3oJysNvzGF8WopwRFwePPw5z50L9+rBwoR1KjY0b7ahAw4ZBbu3zU8Htsv8CjDGfZlchKvtNmWJbkUlJMGqUva879X7uYcPstIqDBjlao1L+wOvZeRHpKyJrROSM6xEjIn2yozjlOydO2MF/ata0jcdBg9xCctkye3H50KFQvLijdSrlDy7bonQF4jPAIGAN9lzlzcBIEcEYM8XnFSqfWL/e3qL48stw443pVr7yCpQubZuYSimvLcrHgU7GmEXGmBPGmOPGmJ+BLq51KkCtXWt/ppn7BuDAATvV7OOP6zWTSrl4C8oixpjY9Atdy4r4oiCVPdavh/Bw23BM4/vv7c/27bO9JqX8lbegPHuV65Sf27gRatXysOK77yAiwtX1rZQC75cHVRORDR6WC1DRB/WobJCSYgfkfeSRdCvWr4d586BfPzvFg1IK8B6UdYBSwJ50y8sD+3xSkfK5PXvg7FmoVs1t4eHD0LEjlCiRbj5apZS3Q+93gJPGmL/cH0CCa50KQAkJ9meRi2eZk5LsQLwHDsCsWR5OXCoV3Ly1KCONMZccehtjYkQk0jclKV9LSrI/8+Z1LRg61N7cPWWKvT1HKZWGtxZl/susK5CVhajsc+6c/Zk3L/ZG7wkT4L77oHdvR+tSyl95C8pVIvJQ+oUi8gCw2jclKV/74w/7s0oV7OhAp09Dw4aO1qSUP/N26P0MMEtE7ufvYIwG8gKdfFiX8qGVK6FoUahUCVi5yy6sqBcxKJURb4NiHASaiMhtQE3X4v+57s5RASomBqKjXePwbtliF2pQKpWhTI2fZYxZBCzycS0qmxw8CHUujl0/aRJERtpJw5RSHunY/kHo9GkoVAh7w/fSpXastdShg5RS6WlQBhljbFBedx3w3nv2yYABTpellF/ToAwy+/dDcjKUvj7Fzv9w//1QrJjTZSnl1zQog8ymTfZnjeuP2vsY69S5/BuUUhqUwSY1KK+LtU/Kl3esFqUChQZlkPn9dyhd2hC+fJZdEBnpaD1KBQKdXi+ImFOn+fk74Y6U+cjwN6FJE6hc2emylPJ72qIMBocOwcsvszniTg6dvo7bInbaqWiXLXMbGUMplRFtUeZku3bB22/bi8rPn2dDvbchBhp+PeTv+6yUUl75tEUpIq1FZJuI7BSRFy6zXX0RuSAiXX1ZT1DZsMGOzPvJJ/YSoC1bONrPztEdFuZwbUoFGJ+1KEUkBBgL3AHEYUcimmOM2exhu/8A831VS1D65BP7c/v21J7tYzPsotBQh2pSKkD5skXZANhpjNlljEkEpgMdPWz3JPA1cMiHtQSX5GSYPh3atUtz+c/mzXamBz0tqdSV8WVQliHtXDtxrmWpRKQMdri2cT6sI/gsXGg7cHr1Sl10+jR8+y106eJgXUoFKF8Gpadp/Ey61+8CQ40xFy77QSIPi0iMiMTEx8dnVX050+LFNiBLlYK2bVMXf/utnSvHLTuVUpnky6CMA8q6vY7g0pkbo4HpIhILdAX+KyL3pP8gY8x4Y0y0MSY6THsiPDMGxo6FVq1sb82SJZAvHwAnT8K//gUVKkDTpg7XqVQA8uXlQauAyiJSAdgL9ADuc9/AGFPh4nMRmQx8Z4yZ7cOacqbz5+GJJ+Djj+15yalTU6dYTEmxnd47d9oj8lx65axSV8xnQWmMSRaRJ7C92SHARGPMJhF51LVez0tmBWNsOC5cCMOGweuvpxlb8tVX4bvv4IMPoHlzB+tUKoCJMelPG/q36OhoExMT43QZ/uP//s/2bP/znzYVXbZts7n5zTfw4IMwfjyIp7PGSikARGS1MSba0zo9EAt0F/+n0bo1APv2wSOPQI0asGCBzc+xYzUklboWegtjoIuJgdy5OVG+NiNegnfesZdRDhwIL70E4eFOF6hU4NOgDHSrVpFYvS51G+UnNhbuuw/eeEMnVVQqK+mhdyC7cAF+/50l5XsTGwuff247vDUklcpaGpSBbONGOHWK/yXdSb58cM89ThekVM6kQRnIli3jOEX5ZmMlWrRwzayolMpyeo4yUO3Zw+nXR9O2wCL2Hwph0hCnC1Iq59KgDERnznC23b10ODqJldTlq6+E2293uiilci4NykCTksKZng/SdcNrLJZbmTJF6NzZ6aKUytk0KAPMoX9PoP3cZ4iR+owfLzoakFLZQIMygOzcCa2Ht2OvlOCbb4SO9zhdkVLBQXu9A8SWLXZ22ePn8vNz21F0vEfvSVQqu2hQBohx4+D0acOvpjGNW+l1QEplJw3KALFoEdxSehdV2AH16jldjlJBRc9RBoD4eHsTTk8+hj594JZbnC5JqaCiLcoAMOvtPwG4o+5hHVhSKQdoUPq7vXuZMPokNfNuo96C4anz4Cilso8GpZ9b95/5xCTfxEODiiBhJZ0uR6mgpEHpx/buhfs+aUmhXGfoNbi00+UoFbQ0KP1UbCzc2iyFuIRQ/td1MsWLO12RUsFLg9IPbd8OzZrBscMXWEgrbu1/o9MlKRXUNCj90FNPwdmzsLjvZBqwCho3drokpYKaBqUf2r0bWrWC2nu/h8qVoWhRp0tSKqhpUPqh/fvh+uuBNWvg5pudLkepoKdB6WdiY+HUKbihZCL89RfUquV0SUoFPQ1KPzNkCBQoAD1bH7MLSpRwtiCllAalP/nlF5gxA154AcoWPm4XFiniaE1KKQ1Kv7FlCzzwAJQrB88/j+3RAQ1KpfyABqUfmD4d6te35yanTYOC33wOXbrYHp369Z0uT6mgp0HpoMREe81kz55Qty6s+e08TT9/DHr3huho2+tdqpTTZSoV9HQ8SoecPg1t2sCyZfDss/Cfx2LJ0/1eiImxPTr//jfk1v88SvkD/ZfogPPnoXNn+O03e6jdsydQqZUdoXfWLLjnHqdLVEq50aDMZhcu2CPrH3+EyZNdIXnyJPz5J7z1loakUn5Iz1Fms2eftZcAjRoFffu6FsbF2Z/lyjlWl1IqYxqU2WjePHj/fXjmGRg0yG3F1q32Z9myTpSllPJCgzKbHD0KDz4INWvaI2yMsVeYd+oEXbtCsWJQtarTZSqlPNCgzCZPPWX7aqZ8kkS+L6fYKWdbtIClS2HYMNi0CUrqVA9K+SPtzMkG69bB1KnwcvOl3NSxGxw4ANWr2xkVe/WyN3crpfyWBmU2GD5gB4UpxXO/tIfWje1Jyjvv1GlnlQoQGpQ+tumTFcxY24ChFWdSbO6vtiWplAooPj1HKSKtRWSbiOwUkRc8rL9fRDa4Hr+KSB1f1pPdfv4qnhYPVaJorlM8u/BuDUmlApTPglJEQoCxQBugOtBTRNInxW6guTGmNvAGMN5X9WQnY2D0yAvc0b04YRxmxdzDhFe4zumylFJXyZeH3g2AncaYXQAiMh3oCGy+uIEx5le37VcAET6sJ1ucOwf9+xmmfxlCZ75m8ieGwm27Ol2WUuoa+PLQuwywx+11nGtZRh4AvvdhPdnio3GG6V8Kb/IiM1/dSOH+GpJKBTpftig9dekajxuK3IYNylsyWP8w8DBAOT+/zW/G6D3U4SgvPnMOXnvT6XKUUlnAly3KOMD9nrwIYF/6jUSkNvAx0NEYc8TTBxljxhtjoo0x0WFhYT4p9pqdO8e+50axfE85ut60y97MrZf/KJUj+DIoVwGVRaSCiOQFegBz3DcQkXLAN0BvY8x2H9biO2fOwOjRUKEC746+AEDXKR0gl970pFRO4bNDb2NMsog8AcwHQoCJxphNIvKoa/044FWgBPBfsa2vZGNMtK9qylKnTsHYsTYk4+NZetNTvH1wMA8/aKhaUy9PVSonEWM8njb0W9HR0SYmJsa5Ao4fhzFj4N134dgxaN2ak8/+gzqPNCIkxN6uWKiQc+Uppa6OiKzOqKGmTZ8rsWEDNG9uw7JDB3j5Zahfn5efgv/7Pzutg4akUjmPBuWVeP55e+5x7Vo7Gxiwdy989BEMGACNGztbnlLKNzQoM2vhQjt/w+jRqSEJMHKknd5h2DDnSlNK+ZZ2zWZGSgq88AKULw+PP566+MQJ25rs0wcqVHCwPqWUT2mLMjOWL4fVq2HiRMiXL3Xx+vX2lsV773WwNqWUz2mLMjN+dd2S3r59msWbXXet16iRzfUopbKVBmVm/P473HjjJVM1rFlje7l1TjClcjYNSm+SkmDFCmjYMM3iLVvsvNydO+udikrldBqUl5OcDPffD/v3Q5cuqYtTUuCRR2xrcuRIB+tTSmUL7czJSEoKPPAAzJgBb79tm44u06bZyRM//hjCwx2sUSmVLbRF6Ykx9jKgKVPg9dfhuefSrF6+HIoXh/79HapPKZWtNCg9mTzZXiD54ov2NsV0kpLsDLM6QJBSwUH/qXvy55+2h+bf//bYU5OcDLn1pIVSQUOD0pPkZMiTx2NIHj8O339vrxZSSgUHDUpPLgalBy+9BIcP2/4dpVRw0KD05MwZyJ//ksWLFsGHH8ITT8BNNzlQl1LKERqUnsTG2gEwXLZssfdz3347lCsHb7zhXGlKqeynQenJ7t1QoQK7d0O/flCzJvzwA7z6qh0Io0gRpwtUSmUn7btNLyWFY7uP83KBQUyIgpAQGDQIhg695FZvpVSQ0KBM78ABhia+zsQNjXj4Edt5U6aM00UppZykQZlOys5dzKEDXZsd5L//Le10OUopP6DnKNNZ9dNJDnI97ds5XYlSyl9oULoxBkZ9XZ7cJNHm/uJOl6OU8hMalG7GjoUZm2rwRoHhFC+dz/sblFJBQYPSZcUK27vd/oYYhkRMc7ocpZQf0c4c4OxZ6NkTIiLg04g3yJWi1wEppf6mLUrgnXfszTgTPzGE/hljE1MppVyCPigPHIDhw6FTJ2hRbB3s2wetWztdllLKjwR9UL78Mpw/DyNGAN99Z4dWa9PG6bKUUn4kqINy7VqYOBGefBIqJW2BTz+FBg2gVCmnS1NK+ZGgDUpj4NlnoXio4ZXkf0Dt2nDkCPzjH06XppTyM0EblHPnGH75BV6/8BLFxrwOffrA9u162K2UukRwXh60bRvv9U+gPKE8HLkAPvwVGjd2uiqllJ8KvqBcv54dTfvx85m1/KvdCnLP/t2OpaaUUhkIrkPvPXugbVvGyaPkzm0YML6RhqRSyqvgaVEePw5t2rDpRATvn3+I++4TSusoakqpTAiOoLxwATp3Jnnbn/S/8SDFjuZi1Cini1JKBYrgCMq1a2HRIt67ZymrZhdh+nSd1kEplXnBcY4yJgaAiZsa0Lw5dOvmcD1KqYASHEG5ahXnit/Atl15uOUWe5eiUkplVnAE5cqVbI7qxIULQp06ThejlAo0Pg1KEWktIttEZKeIvOBhvYjIGNf6DSJyc5YXcfAg/PEHW8u2AqBGjSz/BqVUDuezoBSREGAs0AaoDvQUkerpNmsDVHY9HgY+zPJCFiwAYF/paECnnlVKXTlftigbADuNMbuMMYnAdKBjum06AlOMtQIoJiJZe3Xj/PkQFsb+XGUoWBCKFMnST1dKBQFfBmUZYI/b6zjXsivdBhF5WERiRCQmPj4+8xUYAz/+CHfcwf4D9gJz7chRSl0pX15H6SmSzFVsgzFmPDAeIDo6+pL1GVcgsHIlJCYyooAdRU0ppa6UL4MyDijr9joC2HcV21yb8uVTP1inwlFKXQ1fHnqvAiqLSAURyQv0AOak22YO0MfV+90IOGGM2e/DmpRS6or5rEVpjEkWkSeA+UAIMNEYs0lEHnWtHwfMA9oCO4EEoL+v6lFKqavl03u9jTHzsGHovmyc23MDDPRlDUopda2C484cpZS6BhqUSinlhQalUkp5oUGplFJeaFAqpZQXGpRKKeWFBqVSSnkh9lLGwCEi8cBfHlaVBA5nczlXwt/rA60xK/h7faA1ZqS8MSbM04qAC8qMiEiMMSba6Toy4u/1gdaYFfy9PtAar4YeeiullBcalEop5UVOCsrxThfghb/XB1pjVvD3+kBrvGI55hylUkr5Sk5qUSqllE8EXFCKSIiIrBWR7zysayEiJ0RknevxqgP1xYrIRtf3x3hY7/speq+9Rkf3o4gUE5GZIrJVRLaISON06/1hH3qr0el9GOX23etE5KSIPJNuG8f2Yybrc/zf80U+HY/SR54GtgAZzae41BjTLhvr8eQ2Y0xG14C5T9HbEDtFb8PsKszN5WoEZ/fje8APxpiurtHxC6Zb7w/70FuN4OA+NMZsA+pC6tTRe4FZ6TZzbD9msj7wj3/PgdWiFJEI4G7gY6druQa+n6I3gIlIEeBW4BMAY0yiMeZ4us0c3YeZrNGftAT+NMakv1HDX/4WM6rPbwRUUALvAkOAlMts01hE1ovI9yJSI3vKSsMAC0RktYg87GF9pqbo9TFvNYJz+7EiEA9Mcp1i+VhErku3jdP7MDM1gvN/ixf1AL7wsNzp/XhRRvWBn+zDgAlKEWkHHDLGrL7MZmuwtyHVAd4HZmdHbek0NcbcjD2sGSgit6Zbn6kpen3MW41O7sfcwM3Ah8aYm4AzwAvptnF6H2amRn/4W8R1WqADMMPTag/LsvVv0Ut9frEPIYCCEmgKdBCRWGA6cLuIfO6+gTHmpDHmtOv5PCCPiJTMziKNMftcPw9hz7k0SLeJ76fo9cJbjQ7vxzggzhjzu+v1TGwopd/GyX3otUZ/+Ft0aQOsMcYc9LDO6f0Il6nPj/Zh4ASlMeZFY0yEMSYS21T/2RjTy30bEbleRMT1vAH29zuSXTWKyHUiUvjic+BO4I90mzk6RW9manRyPxpjDgB7RCTKtaglsDndZo7uw8zU6PTfopueZHxY6w/TRWdYnx/tw4Ds9U5D0k5/2xV4TESSgbNAD5O9V9SXAma5/tvmBqYZY34Q/5qiNzM1Or0fnwSmug7LdgH9/WwfZqZGp/chIlIQuAN4xG2Z3+zHTNTn+D5MrUvvzFFKqcsLmENvpZRyigalUkp5oUGplFJeaFAqpZQXGpRKKeWFBqVSSnmhQamUUl5oUKqrJiKRIpL+zqOs+uxfr2Db01n83cVE5HG315n+PUWkgIj84ho6LKNt8orIEhEJ+Bs+goUGpfJLxpgmDn59MeBxbxtlYADwjTHmQkYbGGMSgZ+A7lf5HSqbaVCqa5VbRD4VO0L2TNdtaWmIyGzXkG6b3Id1E5H6rvfld92DvklEarrWnXYt+59rmK0/ROSywSIivURkpdjRsD+62KpztQi3iMgE13csEJECrnWviB2l/EcR+UJEngfeAm50fc5I18eHeHq/B/cD37rVtE5ESonIv0Skr9hRu6djR8K5P7M7WTlLg1JdqyhgvDGmNnASzy2xAcaYekA08JSIlAAwxqzCDszwL2AE8Lkxxv0QtzWwzxhTxxhTE/ghoyJEpBq2hdbUGFMXuEDaIKoMjDXG1ACOA11EJBroAtwEdHbVB3bItD+NMXWNMYMzer+HGvICFY0xsa7XuYHirpFx6gDrgdqun38A9TP6fZR/0aBU12qPMWa56/nnwC0etnlKRNYDK7DDelV2W/c6dmCEaGxYutsItBKR/4hIM2PMicvU0RKoB6wSkXWu1xXd1u82xqxzPV8NRLpq/dYYc9YYcwqYe5nP9/T+9EpiQ/SiqthpSwCqY0cYqgNscB2aJ14cyUn5Nw1Kda3Sj6piRGSg/D0h1H1AK6CxawDWtUB+t+2LA4WAwumWY4zZjg2/jcBwufzkUgJ86moF1jXGRBljXnNbf97t+QXsyEmeBq7NiKf3p3eWtL9DFLBNRIoDp13nJqOBixO65QPOXUENyiEalOpalZO/ZyDsCSwzxoy9GFjY0b+PGWMSRKQq0Cjd+8cDrwBTgf+4rxCRG4AEY8znwNtcOoCvu5+AriIS7npvcREp76X2ZUB71znSQtj5mABOYYP7ihhjjmHPZV4My0RsqzIaWC8ivYBYY8xB1+mHeGNM0pV+j8p+enmCulZbgL4i8hGwAzuTn7sfgEdFZAOwDXv4DYCI9AGSjTHTXB0vv4rI7caYn12b1AJGikgKkAQ8llERxpjNIvIydi6gXK7tBwIZTlhljFklInOw5wz/wrb0ThhjjojIctclQd8DYzO9N2AB9pB+oet3vwP7P4GzwFGgj2u727DjQaoAoONRqqAmIoWMMaddvfVLgIeNMWuu4fNuAgYZY3q7LZuEHSD5R7dl3wAvGjttq/Jzeuitgt14V+fPGuDrawlJAGPMWmBRugvOawMbLr5w9Y7P1pAMHNqiVEopL7RFqZRSXmhQKqWUFxqUSinlhQalUkp5oUGplFJeaFAqpZQXGpRKKeWFBqVSSnnx/5OqUulz4+IzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "plt.plot(sorted(b_axes_KMS), np.linspace(0, 1, len(b_axes_KMS)), color='r', label='manual KMS (n={})'.format(len(b_axes_KMS)))\n",
    "plt.plot(sorted(b_axes_AIF), np.linspace(0, 1, len(b_axes_AIF)), color='b', label='auto AIF (n={})'.format(len(b_axes_AIF)))\n",
    "plt.xlabel('b-axis length ($\\psi$)')\n",
    "plt.ylabel('CDF')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77ed7802",
   "metadata": {},
   "source": [
    "Pretty similar results, looks like we might be able to trust out automatic measurements for this area. Read more on that topic here: https://doi.org/10.1029/2021JF006260."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}