07_wss_imd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Within Sum of Squares and IMD analysis\n",
    "\n",
    "This notebook analyses resulting clusters using WSS and a relationship between clusters and the IMD.\n",
    "\n",
    "The `darribas/gds_py:6.1` Docker container was used to create a computational environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "massive-ladder",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dask.dataframe as dd\n",
    "import geopandas as gpd\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "arabic-sacrifice",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('sentinel_gb_encoder_map_3_20210317_smoothed.csv', index_col=0)\n",
    "labels = pd.read_csv(\"sentinel_gb_encoder_map_3_20210317_smoothed_10k_labels.csv\", index_col=0,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unlikely-tradition",
   "metadata": {},
   "outputs": [],
   "source": [
    "centroids = data.groupby(labels.cluster.values).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "environmental-boards",
   "metadata": {},
   "outputs": [],
   "source": [
    "squared_error = pd.concat([(data.loc[labels.cluster == i] - centroids.loc[i]) ** 2 for i in range(7)]).sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "numerous-shadow",
   "metadata": {},
   "outputs": [],
   "source": [
    "squared_error_sum = squared_error.sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "numerous-guatemala",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    1.239086\n",
       "1    0.256577\n",
       "2    0.359320\n",
       "3    1.427487\n",
       "4    0.852287\n",
       "5    0.220248\n",
       "6    0.322644\n",
       "dtype: float64"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean_squared_error = squared_error_sum.groupby(labels.cluster.values).mean()\n",
    "mean_squared_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "knowing-fisher",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.set(context=\"paper\", style=\"ticks\", rc={'patch.force_edgecolor': False})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "presidential-cable",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels[\"error\"] = squared_error_sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32e15742-d6e6-4a39-88da-bedb508dd328",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_counts = labels.cluster.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "297be53a-410a-4594-bd4b-c9264cf2e1d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_names = \"ABCDEFG\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "perfect-validation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAICCAYAAADbH1mbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1jElEQVR4nO3de3xU9Z3/8XcmASKEXAEV5BchLAEbQKxCtVwUETBgAS/LNRSrIK4EBaoJRVBArgvog6DcBGpBgXUFuhbBqmgRSxWqQOUSlkgQiQi5J4CZZHJ+f7AZGfMNQwIzJwmv5+ORh8ycOef7mY9nJu+ca4BlWZYAAADgwWF3AQAAANURIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAxqXEiKjY21uwQAAHANqHEhCQAAwB8ISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADQhIAAIBBkN0FVBfLly9XWlpalefPyMhQ06ZNqzRvTEyMRo8eXeWxAQDA1RdgWZZldxGVERsbq9TUVLvLKCcpKUlz5861uwwAAHCVsLsNAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADn11M8ocfftDYsWNVr149lZSU6MUXX1SbNm3c05OTk3XkyBE1aNBALVq00PTp031VCgAAQKX5LCQ1atRIGzZskMPh0K5du/T6669r/vz5Hq+ZNm2a2rVr56sSAAAAqsxnu9sCAwPlcFxY/NmzZ3XLLbeUe8306dOVkJCgnTt3+qoMAACAKvHpvduOHj2q559/Xt9//71SUlI8piUlJSkiIkJZWVkaOXKkOnTooIYNG/qyHAAAgMvm0wO3W7VqpfXr12vp0qWaMWOGx7SIiAhJUlRUlOLi4nTs2LFy86ekpCg2NtbjBwAAwB98FpKcTqf736GhoQoODvaYXlBQIEkqKirSwYMH1axZs3LLSExMVGpqqscPAACAP/hsd9v+/fv1yiuvKCAgQNKFs9l27NihwsJCxcfHa8KECSooKFBJSYkSEhIUFRXlq1IAAAAqzWch6fbbb9fatWsrnL5ixQpfDQ0AAHDFuJgkAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADDw2Q1uAeBKLF++XGlpaVWePyMjQ02bNq3SvDExMRo9enSVxwZQOxCSAFRLVxpSkpKSNHfu3KtUDYBrEbvbAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAwISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADQhIAAICBz0LSDz/8oEceeUTDhw/X4MGDdfjwYY/px48fd09bv369r8oAAACokiBfLbhRo0basGGDHA6Hdu3apddff13z5893T1+wYIGSk5PVpk0bDRkyRH369FF4eLivygEAAKgUn21JCgwMlMNxYfFnz57VLbfc4jE9PT1dcXFxCgoKUqdOnbR//35flQIAAFBpPtuSJElHjx7V888/r++//14pKSke0yzLcv87LCxMeXl5viwFAACgUnx64HarVq20fv16LV26VDNmzPAc2PHT0Pn5+QoLCys3f0pKimJjYz1+AAAA/MFnIcnpdLr/HRoaquDgYI/p0dHROnDggFwul3bv3q327duXW0ZiYqJSU1M9fgAAAPzBZ7vb9u/fr1deeUUBAQGSpOTkZO3YsUOFhYWKj4/XxIkTNXnyZJWUlGjgwIEctA0AAKoVn4Wk22+/XWvXrq1wenR09CWnAwAA2ImLSQIAABgQkgAAAAwISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAwISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAyCfLXgtLQ0Pf/883I4HHI4HJo1a5aaN2/unp6cnKwjR46oQYMGatGihaZPn+6rUgAAACrNZyEpIiJCy5YtU2hoqHbs2KGlS5dq5syZHq+ZNm2a2rVr56sSAAAAqsxnu9siIyMVGhoqSQoKClJgYGC510yfPl0JCQnauXOnr8oAAACoEp9tSSpz/vx5LVq0qNxWpKSkJEVERCgrK0sjR45Uhw4d1LBhQ4/XpKSkaPHixb4uEQAAoByfHrhdUlKi8ePH6/HHH1dMTIzHtIiICElSVFSU4uLidOzYsXLzJyYmKjU11eMHAADAH3wWkizL0uTJk9WtWzf17Nmz3PSCggJJUlFRkQ4ePKhmzZr5qhQAAIBK89nutk8//VTbtm1TRkaGtm7dqjZt2qhr164qLCxUfHy8JkyYoIKCApWUlCghIUFRUVG+KgUAAKDSfBaSunXrpn379lU4fcWKFb4aGgAA4IpxMUkAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAwISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAQZDdBQAAro7ly5crLS2tSvNmZGSoadOmVZo3JiZGo0ePrtK8QHVGSAKAWuJKgkpSUpLmzp17FasBaj52twEAABgQkgAAAAwISQAAAAaVCkmlpaUqKCjwVS0AAADVhteQ9Mwzz6iwsFCFhYWKj4/XwIEDtXTpUn/UBgAAYBuvISk9PV0hISH64IMP1KtXL/31r3/Vli1b/FEbAACAbbyGpOLiYp08eVLvvvuuevToIYeDw5gAAEDt5zXxPPXUUxozZoyaN2+uW2+9VSdOnFB0dLQ/agMAALCN14tJxsfHKz4+XllZWZKk5s2ba/HixT4vDAAAwE5etyTt3LlTffv21SOPPCJJOnz4sJ577jmfFwYAAGAnryFp4cKFevPNNxUaGipJatOmjQ4dOuTzwgAAAOzkNSQFBgYqPDzc/diyLF/WAwAAUC14PSapffv2Wrt2rYqLi7V371699dZb6tKliz9qAwAAsI3XLUmTJk2Sw+FQTEyMVq5cqVtvvVXPPvusP2oDAACwzSW3JLlcLj388MPavHmzhg4d6q+aAAAAbHfJLUmBgYFq27at0tLS/FUPAABAteD1mKTU1FT1799fLVu2VP369WVZlgICArR+/Xp/1AcAAGALryEpJSXFH3UAAABUK15DUrNmzeRyuXTq1Ck5nU5/1AQAAGA7ryFp69atevXVV5WRkaGbb75ZR44cUVxcHLvbAABAreb1EgCvvfaa1q9fr+bNm2vjxo16++23ddNNN/mjNgAAANt4DUl169ZVSEiIJKm4uJiz3QAAwDXB6+62xo0bKz8/Xz169NCYMWMUHh6uqKgof9QGAABgG68haenSpZKkp59+Wp9//rkKCwvVtWtXnxcGAABgJ68h6dixY+5/N2nSRE2aNNHJkyfVokULnxYGAABgJ68haerUqe5/FxcXKzU1Va1bt9aGDRt8WhgAAICdvIakNWvWeDzOyMjQwoULfVYQAABAdeD17Lafu/HGG3X48GFf1AIAAFBteN2SNGHCBAUEBEiSSktLdfToUf3qV7/yeWEAAAB28hqSBg8e7P63w+FQ06ZN1bRpU58WBQAAYDevIalTp05VWnBaWpqef/55ORwOORwOzZo1S82bN3dPP378uCZPnqySkhINGDDAI4wBAADYzWtI6tKli7KyshQcHCxJ+vHHHz0uJrlz507jfBEREVq2bJlCQ0O1Y8cOLV26VDNnznRPX7BggZKTk9WmTRsNGTJEffr0UXh4+BW+HQAAgKvDa0jq1auXunXrprvvvluS9PHHH2vPnj169tlnLzlfZGTkT4MEBSkwMNBjenp6uuLi4iRd2Fq1f/9+devWrbL1AwAA+ITXs9t2797tDkiSdM8992jXrl2XPcD58+e1aNEi/fa3v/V43rIs97/DwsKUl5dXbt6UlBTFxsZ6/AAAAPiD1y1J0dHRmjdvnuLj4yVJW7Zs0Y033nhZCy8pKdH48eP1+OOPKyYmxmOaw/FTPsvPz1ebNm3KzZ+YmKjExESP5whKAADAH7xuSfrP//xPRUZGasmSJXrttdfUqFEjzZ8/3+uCLcvS5MmT1a1bN/Xs2bPc9OjoaB04cEAul0u7d+9W+/btq/YOAAAAfMDrlqTrrrtOjz/+uCQpNzdXp06d0nXXXed1wZ9++qm2bdumjIwMbd26VW3atFHXrl1VWFio+Ph4TZw40X1228CBAzloGwAAVCteQ9LQoUO1fPlyOZ1OPfzww2rWrJlatWqlF1544ZLzdevWTfv27atwenR0tNauXVv5igEAAPzA6+628+fPKyQkRNu2bdMjjzyiNWvW6Msvv/RHbQAAALbxGpJKSkr0xRdfaNOmTbrnnnskSS6Xy+eFAQAA2MlrSJoyZYpWrVql++67T23atNGJEyfUuXNnf9QGAABgm8u6LcnFtyZp3ry5pkyZ4tOiAAAA7OZ1SxIAAMC1iJAEAABgUGFIGj58uCR53JQWAADgWlHhMUm5ubnavn27/va3v6l79+7lpnfp0sWnhQEAANipwpA0ZcoUvffee8rLy9OWLVvKTSckAQCA2qzCkNS5c2d17txZHTt21IABA/xYEgAAgP28XgKga9eumjJlivsq27/85S/1zDPPKDIy0ufFAQAA2MXr2W1JSUlq27atNmzYoA0bNqht27Z69tln/VEbAACAbbyGpNOnT2vo0KEKCQlRSEiIhgwZojNnzvijNgAAANt4DUk33HCD1qxZo+zsbOXk5Gjt2rW6/vrr/VEbAACAbbwekzRnzhylpKTod7/7nSTptttu09y5c31eGAAAgJ28hqTIyEi98MIL/qgFAACg2vAakgCgqqYmJSnnTKYtY2fm5yvxd4/5fdyIxo00na3tQK1ASALgMzlnMnV/Sak9g9cPkWwYe6tNoRDA1ccNbgEAAAy8bkk6cOCAVq1apVOnTsnlcsmyLAUEBGj9+vX+qA8AAMAWXkPShAkT9Nxzz6l169ZyONjwBAAArg1eQ1JUVJTuvfdef9QCAABQbXgNSb///e81btw4derUSXXq1HE/P2jQIJ8WBgAAYCevIem1115TcHCwsrOz2d0GAACuGV5D0g8//KB3333XH7UAAABUG143Dd1333167733dPbsWTmdTvcPAABAbeZ1S9LmzZslSfPnz3c/FxAQoI8++shnRQEAANjNa0javn27P+oAAACoVryGpA0bNhif5+w2AABQm3kNSWfOnHH/2+l06rPPPlOLFi0ISQAAoFbzGpLGjh1b7vF//Md/+KwgAACA6sBrSPq5nJwcffvtt76o5Yo9l/y8MrNybBk7PzdLvxv1lN/HbRQVoXlzXvL7uAAA1HZeQ1KXLl08Hjdo0ECJiYk+K+hKZGbl6MeIu2wZu26E9KMN42Zm/d2GUQEAqP28hqSdO3f6ow4AAIBqxevFJHfv3q1z585JkjZt2qTZs2frxIkTPi8MAADATl5D0owZM1S/fn3961//0tq1a9W6dWslJSX5ozYAAADbeA1JAQEBkqQPP/xQCQkJeuihh9xblgAAAGorr8ckNW7cWJMmTdLnn3+ujRs3yul0yuVy+aM2AAAA23jdkpSSkqJ77rlHq1atUnh4uHJzc9ndBgAAaj2vW5Kuu+469erVy/24SZMmatKkiU+LAgAAsJvXLUkAAADXIkISAACAwWXdluTbb7/VyZMnPQ7Y/vmVuAEAAGoTryFp2rRp+uqrrxQbGyuH46cNT4QkAABQm3kNSf/4xz+0ZcsWj4AEAABQ23lNPnFxcTp58qQ/agEAAKg2vG5JSktLU79+/RQTE6O6devKsiwFBARo/fr1/qgPAADAFl5DUkpKij/qAAAAqFYqDEkZGRlq2rSpnE6nP+sBAACoFioMSUuWLNGMGTM0derUctMCAgL0pz/9yaeFAQAA2KnCkDRjxgxJ0po1a6q0YKfTqREjRujo0aN66aWX1KdPH4/pycnJOnLkiBo0aKAWLVpo+vTpVRoHAADAF7wek3T27Fm99957ysjI8LiY5IQJEy694KAgLVq0SBs2bKjwNdOmTVO7du0qUS4AAIB/eA1JY8aMUWxsrNq2bVupayU5HA6vN8KdPn26goOD9cQTT3BxSgAAUK14DUn5+fl6/vnnr/rASUlJioiIUFZWlkaOHKkOHTqoYcOGHq9JSUnR4sWLr/rYAAAA3lS4acjpdMrpdKpnz5569913VVhY6H7uapzxFhERIUmKiopSXFycjh07Vu41iYmJSk1N9fgBAADwhwq3JPXp00cBAQGyLEuS9PLLL7sfBwQE6KOPPrqigQsKCtSwYUMVFRXp4MGDatas2RUtDwAA4GqqMCRt375d0oUDtxs0aOAx7ezZs5e18MTERB08eFD169fX3r17ddddd6mwsFDx8fGaMGGCCgoKVFJSooSEBEVFRV3B2wAAALi6vB6TNHz4cG3atMnrcyaXulr3ihUrLqM8AAAAe1QYkk6cOKH09HSdPXtWO3fudD9fWFiooqIivxQHAABglwpDUmpqqj766CPl5eVpy5Yt7ufr16/PhR8BAECtV2FI6tmzp3r27Kmvv/5acXFx/qwJAADAdhWGpIULF2rChAlatWqVAgICyk1fsGCBTwsDAACwU4UhqVevXpKkwYMH+60YAACA6qLCkFS2iy0/P18dO3bkFH0AAHBN8XoJgK+++kqrV69Wfn6+2rZtq44dO6pjx45q06aNP+oDAACwhdeQ9Oyzz0q6cOr/hg0btHTpUmVmZurAgQM+Lw4AAMAuXkPSvHnzdOjQIQUGBqpdu3aaMWOGOnbs6I/aAAAAbFPhDW7LHDt2TE6nU82aNdPNN9+sFi1aqGHDhv6oDQAAwDZetyQtWbJEkpSenq6PPvpIw4cPl8vl8rgKNwAAQG3jNSRt3rxZX375pfbt26e6deuqd+/euu222/xRGwAAgG28hqRTp06pX79++sMf/qDg4GB/1AQAAGA7ryFpzJgx/qgDAACgWvF64DYAAMC1iJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAwISQAAAAaEJAAAAIMguwsAUHtl5udro9Npdxl+dbZuXbtLAHCVEJIA+Eyj0FDdX1Jqdxl+tTWIDfRAbcGnGQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgwG1JAKAa+cOkKcrKyvH7uDm5mXpi9Fi/jytJUVERmjV7hi1jA5dCSAKAaiQrK0fRje/2+7jRjf0+pNvxM5/YNzhwCexuAwAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAgc9CktPp1ODBg3X77bdr27Zt5aYfP35cw4cP1+DBg7V+/XpflQEAAFAlPgtJQUFBWrRokX77298apy9YsEDJyclau3at3nnnHeXm5vqqFAAAgErzWUhyOBxq0qRJhdPT09MVFxenoKAgderUSfv37/dVKQAAAJVm2zFJlmW5/x0WFqa8vDy7SgEAACjHttuSOBw/5bP8/Hy1adOm3GtSUlK0ePFif5YFAAAgycYtSdHR0Tpw4IBcLpd2796t9u3bl3tNYmKiUlNTPX4AAAD8wadbkhITE3Xw4EHVr19fe/fu1V133aXCwkLFx8dr4sSJmjx5skpKSjRw4ECFh4f7shQAAIBK8WlISklJqXBadHS01q5d68vhAQAAqoyLSQIAABgQkgAAAAwISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGQXYXcDXl52apKOcDu8vwK2fAj3aXgMu0fPlypaWlVWnejIwMNW3atMpjx8TEaPTo0VWeHwCuRbUqJIWGR+nHiLvsLsOvgnP+bncJuExXElKSkpI0d+7cq1iNf0Q0bqStZzJtGTszP1+NQkP9Pm5E40Z+HxOAb9SqkASgepluY7CrqcESQPXBMUkAAAAGhCQAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAwISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAMCEkAAAAGhCQAAACDILsLAAD8JCc3U3m579ldhl+V6pzdJQBGhCQAqEYiwhspuvHddpfhV8fPfGJ3CYCRT0PS+vXrtXnzZgUFBWnmzJmKjo52T0tOTtaRI0fUoEEDtWjRQtOnT/dlKQAAAJXis5CUm5urd955R+vXr9fBgwe1YMECLVq0yOM106ZNU7t27XxVAgAAQJX57MDtffv2qXPnzgoMDFS7du2Unp5e7jXTp09XQkKCdu7c6asyAAAAqsRnW5Ly8/MVGhrqfmxZlsf0pKQkRUREKCsrSyNHjlSHDh3UsGFDj9ekpKRo8eLFvioRV2j58uVKS0ur0rwZGRlq2rRpleaNiYnR6NGjqzQvAKBqruQ7X6qZ3/s+C0mhoaE6cuSI+7HD4bnRKiIiQpIUFRWluLg4HTt2TO3bt/d4TWJiohITEz2ei42N9VHFqKwrWWGTkpI0d+7cq1gNAMCXrjSk1MTvfZ/tbuvQoYO++OILuVwuHThwwOOgbUkqKCiQJBUVFengwYNq1qyZr0oBAACoNJ9tSQoPD9eAAQM0bNgw99ltO3bsUGFhoeLj4zVhwgQVFBSopKRECQkJioqK8lUpAAAAlebTSwAMGTJEQ4YMcT++eGvSihUrfDk0AADAFeG2JAAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMPDpJQD8rVFUhDKz/m7L2Pm5WQoN9/+1nhpFRfh9TACoDey6tZLE7ZVqiloVkubNecm2sWvi5dYB4FrGrZXgTa0KSYCvJT2frMycLL+Pm5+Vq8eeGuX3cSWpUUSU5r40x5axAcBOhCSgEjJzsuS6vb7fx22g+nL5fdQLMvf4PxQCQHXAgdsAAAAGhCQAAAADQhIAAIABIQkAAMCAkAQAAGBASAIAADAgJAEAABgQkgAAAAy4mCQAVCNRURE6fuYTv4+bk5upiPBGfh9XuvCegeqIkAQA1cis2TNsGZd7kQHlsbsNAADAgJAEAABgwO62a9gLf0hSXnamLWOfzsnXM2Mes2XssMhGmjaL3QoAgEsjJF3D8rIzldAywKbRw2waV1rzjT3BEADs9uLkqcrLzrNl7NPZZzT+yaf9Pm5YZJhenDm9SvMSkgAAuEbkZedpZLuBdpfhV3/816Yqz8sxSQAAAAaEJAAAAAN2twEAaqw/TEpSdlaW38fNyc3TmNGP+31cSYqMitKs2Zx84g+EJABAjZWdlaXWUYH+Hzgq0v9j/p8jNoTCaxUhCaiE/KxcOT/JtrsMv6pbxF55ANcmQhJQCaFR4XLdXt/uMvwqcM85u0sAAFvwJyIAAIABIQkAAMCA3W3XsNM5+Vr5VbHdZfhdbmkdu0sAANQAhKRrWJOIUBtvS2KfNd9YdpcAAKgBCEkAAFwjTmef0fK/r7O7DL/KLS6o8ryEJABAjZWTm6c9OSV2l+FXroCq/+puEtmYe7dVAiEJAFBjRYSH2XMxSRsdyXLZXcI1g7PbAAAADNiSBFRCo4goZe7x/y0B8rNyFRoV7vdxpQvvGQCuRYSka1hYZCOt+SbTlrFP5+SrSUSoLWOHRTaq8rxzX5pzFSu5fElJSZo7lxtaAoA/EZKuYdNm2fdLl1/6AOB/YZFhV3Qg85U4nX1GTSIb+33csMiwKs9LSAIA1FiRUVE6kuX/XeA5uXmKCK/6L98rERlV9V3gL86cfhUrqZya+McxIQlAtbR8+XKlpaVVef6MjAwlJSVVad6YmBiNHj26ymPDf2bNtueXbk38hY/KIyQBqJYIKQDsxiUAAAAADAhJAAAABoQkAAAAA45J+j8cJAoAAC5GSPo/hBQANd2V/LHHH3pAeT4NSevXr9fmzZsVFBSkmTNnKjo62j3t+PHjmjx5skpKSjRgwAANHjzYl6UAQK1HUAGuLp8dk5Sbm6t33nlHb775ppKSkrRgwQKP6QsWLFBycrLWrl2rd955R7m5ub4qBQAAoNJ8tiVp37596ty5swIDA9WuXTulp6d7TE9PT1dcXJwkqVOnTtq/f7+6devmq3LgA2zarxy7+iXV3J4BvsRnsnKuxWN3fRaS8vPzFRr60w1MLcvymH7x47CwMOXl5ZVbRkpKihYvXuyrEnGFatoH3G70C6he+ExWzrXYL5/tbgsNDVVBQcFPAzk8h7r4cX5+vsLCyt8DJzExUampqR4/AAAA/uCzkNShQwd98cUXcrlcOnDggMdB25IUHR2tAwcOyOVyaffu3Wrfvr2vSgEAAKg0n+1uCw8P14ABAzRs2DD32W07duxQYWGh4uPjNXHiRPfZbQMHDlR4eLivSgEAAKi0AOvnBwtVc7Gxsex2AwAAPsdtSQAAAAwISQAAAAaEJAAAAANCEgAAgAEhCQAAwICQBAAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwCLK7gKqIjY21uwQAAFBLpKamGp8PsCzL8nMttVJsbGyFTUZ59Kty6Ffl0bPKoV+VQ78qryb2jN1tAAAABoQkAAAAA0ISAACAASHpKhk7dqzdJdQo9Kty6Ffl0bPKoV+VQ78qryb2jAO3AQAADNiSBAAAYEBIAgAAMKiRF5Osbp555hmVlpZq0aJFdpdS7X333Xd66KGH1Lp1a507d0733HNPjdxP7U8HDx7U/Pnz5XQ65XQ6NXToUA0YMMDusqqlsvUrNjZWAQEBuuGGGzRx4kQ1adLE7tKqrYs/k5IUERHBd5kXP+9Zv379NGjQIJurqr6+/vprLVy4UE6nUy6XS506ddL48ePtLuuyEJKu0NmzZ5WdnS3LslRYWKiQkBC7S6r2br31Vi1btkyWZWnAgAEaMWKEQkND7S6rWiooKFBSUpKWLFmim266SSUlJdq9e7fdZVVrZeuXJH3yySeaOHGi1qxZY3NV1dvFPcPloWeXJz8/X5MmTdLSpUvVrFkzSdKuXbtsrurysbvtCn3wwQfq3bu3+vTpo7/+9a92l1OjFBcXS5KCgsjqFfn444/Vo0cP3XTTTZIu9OrOO++0uaqa4+6775bL5dKpU6fsLgW4Jn3yySe699573QFJUo36DiMkXaH3339f999/v/r27UtIukx79+5VQkKC+vTpo1/+8peqX7++3SVVW6dOndINN9xgdxk12vXXX6/Tp0/bXUa1VvaZTEhI0IwZM+wup0a4uGeHDh2yu5xq6+LvsOPHj7u/+/Py8myu7PLwJ/wVyM7O1tdff61nn31WknTkyBFlZmaqUaNGNldWvZVtpi4tLdW4ceP01VdfqWPHjnaXVS3dcMMNSktLs7uMGu3777/X9ddfb3cZ1Rq7jiqPnl2ei7/DoqOjtWbNGiUkJMjlctlc2eVhS9IV2Lp1q8aNG6eVK1dq5cqVmjBhgrZu3Wp3WTWGw+FQw4YNlZ2dbXcp1dY999yj7du367vvvpMklZSU1Kj9+Xb76KOPFBgYSEgCbHL33Xd7fIdJF77Hagq2JF2BLVu2aP78+e7HXbp0UWJiohISEmysqvor20xdUlKipk2bqlu3bnaXVG01bNhQc+fO1dSpU91ntw0fPtzusqq1vXv3asSIEZIu/BW7cOFCmyuq/so+k5JUp04drVq1yuaKUFuEhoZq9uzZmjJlikpKSuRwOHTbbbepYcOGdpd2WbjiNgAAgAG72wAAAAwISQAAAAYck/R/du/erb1792rUqFFyuVzq27evHnnkET322GOSLpzJ0K5dO0nSiBEjdN999+mTTz7RvHnzdPr0ae3Zs8e9rE2bNmnNmjUKDAzUfffdp9GjR+vs2bNKSkpSXl6emjRpolmzZqlevXpKSUlR9+7d1b59e1ved2VUpUcbN27Ua6+9phtvvFGStGLFCgUHB2vfvn2aPXu2LMvS6NGjde+9916yR9u2bVNkZKRCQkK0ZMkSlZSU6Omnn1ZKSoocjuqb9ct6dvToUR0/flwul0tDhw7VwIEDJUnr16/X5s2bFRQUpJkzZyo6OrrC9SohIUHFxcWqU6eO7rjjDo0bN849bfPmzZo2bZq++uorj/FHjhypli1baurUqTWmZ2XKehcbG6tXXnlFdevW1S233KKpU6eqtLRUkydP1okTJxQSEqLZs2crIiJCGzZs0KZNmyRJzZs31+zZsxUUFKQxY8aooKBA58+f12OPPaa+fftKunANl9WrV8vlcunhhx/WgAEDatRnUqraOlaZPn3//fd66qmndOzYMf3pT39yf8bL1LR1rKxfJ06c0NGjR/Xjjz/qN7/5jUaOHKlDhw5p1qxZki5cKNiyLG3atEnnz59XcnKyzpw5o1atWunFF1+Uw+FQVlaWpk2bptzcXDVq1EgLFy7U7t27NX/+fNWpU0fBwcGaP3++wsPD3ePXpn6V2bNnj4YNG6Zdu3YpMjKywu99U78qWr9SUlL097//XcXFxbrjjjuUlJRkT78sWJZlWU888YSVl5dnWZZlbdiwwRo1apT1+uuvu6f37du33Dw5OTnW+fPny03r3bu3VVhYaLlcLuuBBx6w8vPzrZUrV1pvvPGGZVmWtXr1auvNN9+0LMuyTp8+bSUmJvrqbV1VVenRO++84/GaMoMGDbJOnTplFRYWWgMGDLCKi4sr7NGiRYusrVu3llvGihUrrA8//PCqvDdfKevZsWPHLMuyrKKiIqt3796W0+m0cnJyrIcfftgqKSmx9u/f714PKlqvhg8fbmVlZZUbo6ioyBozZox1//33ezz/2WefWaNGjbKmTZvmfq4m9KxMWe8GDRpknTx50rIsyxo9erR16NAha9u2bdasWbMsy7KsDz74wJo3b55lWRd6UebZZ5+1PvvsM4/nCwoKrN69e1uWZVmZmZnW2LFjrZKSEo9xa9Jn0rKqto5Vpk8//vijlZ2dbSUlJVn79+/3GLsmrmNl/Sp7r8XFxVavXr08emJZF76DlixZYlmWZa1Zs8ZauXKlZVmWNW3aNOvjjz+2LMuyfv/731vHjx/3mO/777+3zp8/b1mWZb311lvWq6++6p5WW/s1duxY68EHH3R/P1X0vW/qV0Xr18XLHzZsmJWenm5Zlv/7VT2jq5+V/eUUGhqqoqIibd++Xb179/Z4zffff6/hw4dr4sSJ7lPWw8PDFRwcXG55rVq10vnz5+V0OlWnTh3VrVtXx48f1y9+8QtJUlxcnHsLQePGjXX69Gk5nU4fv8srU9UeSdLbb7+toUOHus+YKSoqUmlpqa6//no1aNBA0dHROn78eIU9kqQlS5Zo6NCh+vOf/+x+rmvXrvrggw98+bavyMU9u/nmmyX9dHVxh8Ohffv2qXPnzgoMDFS7du2Unp4uqeL1SpLGjh2r3/3ud/r666/dz23YsEEPPvhgub+s/vSnP2nYsGEez1X3npW5uHetW7dWfn6+XC6XioqKFB4eXuG6UrduXUmSZVkKCAjQ//t//8/j+XPnzik2NlaS9Le//U0NGjTQqFGj9OSTT7qvyl1TPpNS1dexyvSpXr16ioiIMI5f09axi/tV9l6LiooUHR2tOnXqeLz2vffec29x3LNnj7p37y5J6t69u/bs2aOSkhJ98803SklJUUJCgvvyLzfccIP78xsUFKTAwED3Mmtjv7Zv3268KPDPv/cr6ldF61fZeMXFxQoNDVVUVJQk//eLkCTp2LFj7ts+rFmzRoMGDVJAQIDHaz788EOtXbtW9957r+bMmXPJ5fXs2VMDBgxQ7969FR8fr3r16ql169bauXOnJOnTTz/1uNrojTfeqG+//fYqv6urq6o96tmzp7Zs2aI33nhDe/bs0a5du5Sbm+tx+mdoaKjy8vIq7NHw4cP15z//WStWrNC6devcX/Q333yzjh496uu3XmUX96zMihUr1LdvXwUGBio/P9/jnnWWlxNNFy1apLfeeksvvPCCJk2aJMuydO7cOe3cuVP33Xefx2u3bt2qLl26lAtb1b1nZS7uXe/evfX444+rT58+atmypW644Qa1bt1an332maTyn6fVq1erT58+ysnJcX+xStKjjz6q/v37q0uXLpKkM2fO6OTJk1q+fLkSEhI0d+5c92trwmdSurJ17HL7VJGauI79vF/PPfecevXqpXbt2nl8n3377bdyOBxq3ry5JHn0MSwsTHl5ecrOztbhw4c1ZswYLVu2TMuWLVNubq57GVlZWXrrrbf0yCOPSKqd/SotLdW6des0ZMgQj/lM3/ve+mUyb9489erVSxEREe4Q5u9+EZIuUlBQoC+++EL33HNPuWllSff+++/X4cOHK1xGYWGhli9frm3btumDDz7Qxx9/rIyMDD388MPKysrSiBEj5HQ6Pe5KXvbXXE1Q2R6FhoYqMDBQderUUa9evXTgwAGFh4ersLDQY5lhYWEV9qhsuQ0aNFCXLl3ctwCoSX2TpP/5n//RoUOH9NRTT0m60JuCggL3dG/72Mv6EB0drYiICOXk5OiPf/xjub9MXS6X3n77bf37v/97uWXUtJ5J0owZM7Rx40a9//77ysnJ0T//+U91795dTZo0UUJCgr777juPi0U++uijev/999WiRQtt3LjR/fzq1av1/vvv6/XXX1dBQYFCQ0PVuXNnBQUF6a677vK4snlN7JNUuXXscvtkUlvWsXnz5mn79u367LPPPH7x/uUvf1G/fv3cjy/uY35+vsLCwhQWFqbrr79eMTExql+/vn7xi1+4g3VhYaHGjRunF198UZGRkbW2X++++6569OihevXqebzO9L1/qX5V5LnnntOHH36ogoICffrpp5L83y8O3JbUokULfffdd/rmm2+Uk5Ojxx57TKdPn1ZxcbHatGmjjh07ql69egoMDNSePXs8btT3cw6HQ3Xq1FH9+vXlcDgUHByswsJCNW3aVNOnT5d0YdfRHXfc4Z7n1KlT5f4arG6q2qOCggL3VqOycFWvXj05HA6dPn1aISEhSk9PV3R0tIKCgow9KluGy+XSl19+qR49eki68Ndey5YtbejG5SnrmSTt2LFDmzZt0rJly9y/qDp06KDXXntNLpdLhw8fVnR09CWXV1hYqJCQEOXk5OiHH35QeHi4vvnmG/3zn//UG2+8oZMnTyopKUkTJkxQTk6OnnzySeXl5SkrK0sdOnRQ//79q33Pylzcu6CgIDVs2FAOh0OhoaHKz89XQECAJk6cKOlCMCi7FZDT6XRvpg8NDVVwcLBKS0tVWlqqoKAgBQcHKzg4WHXr1lWnTp3cF4M9dOiQx+e6JnwmpaqvY5Xpk0lmZmaNXMcu7ldZD+rWravrrrvO4xf9tm3btHr1avfjO+64Q3/729/UsmVL7dixQ126dFG9evXUtGlTnT59WlFRUfrf//1fNWvWTE6nU08//bQee+wxdejQQVLt7deRI0d04MABffjhh0pNTdX48eP1xhtvVPi9b+pXRcrGCwwMVEhIiPv/j7/7RUjShasaX3fddWrZsqXefvttSdLGjRuVk5OjX//61/r66681ZcoU1a9fX4GBge5f5Pv27dPLL7+skydPauTIkXr88cfVpUsXDRw40L076tZbb1Xr1q2Vmpqql156SQ6HQ7fffrv7F31mZqaaNGlSLolXN1Xt0apVq7Rz504FBgbqlltuUc+ePSVJycnJGjdunEpLSzV27FgFBQVV2KM5c+YoLS1NLpdLPXr00C233CLpwm6Wn+9mqk7KelZQUKBJkyapSZMm7jMBFy5cqMaNG2vAgAEaNmyY+8wjybxe3XXXXRoxYoTq1aunkpISTZo0SQ6Hw+OK7/369XPvMio7c+nzzz/X+++/r/79+0uq/j0rc3HvnnzySSUkJKhOnTpq2rSpunbtquzsbD399NNyOBxq1aqVkpOTJUmvvvqqvvzyS5WWluqmm27SqFGjdO7cOT355JOSLhzf8Nhjj6levXqKiYlR+/btNXz4cFmW5V5na8pnUqr6OlaZPhUVFemJJ57Q0aNHlZaW5t79WRPXsYv7NXHiRPexo71793bvWjt8+LAaN27ssQvywQcfVHJysoYOHapWrVq57xKQnJys8ePHq7i4WP3791dUVJTWrVunf/3rX1q9erVWr16trl27avTo0bWyX2X3LZUunH378ssvS7r09/7P+1XR+vXiiy/qxIkTKi4uVseOHfWrX/1Kkg398tsh4tXc559/bi1btszv46akpFh79+71+7hVYVePTIqLi60nn3yy3JlJ1Q09qzo+k5eHdaxy6FflXOv94rYkAAAABhy4DQAAYEBIAgAAMCAkAQAAGBCSAAAADAhJAGzRo0cPFRUVVWqe/Px89yUoqiolJUXr1q275GsqW9uvf/3rK6oJQPVESAJQY1QlJJWUlPioGgC1HSEJgM9t2LBBDzzwgH7zm99o8eLFHtM+//xzjR8/3v04ISFBaWlpOn36tAYPHqz+/fvrgQceUFpamhYtWqQjR46of//+ev3111VcXKzp06froYce0oABA7Rjxw5JF7YWTZ06VSNGjHBf4M5k0aJFeuihhzwuxFkmJSVF/fr106hRo9y30fnqq680ZMgQDRw4UBMmTKj0ljAANQtX3AbgU4cPH9Zbb72ldevWKSQkxOtNLcv85S9/0V133aVx48apuLhYLpdL48aNU3p6uv7rv/5LkvTmm2/q5ptv1tSpU5WTk6Nhw4apa9eukqS0tDT98Y9/LHd394uNGDHC48rv+/btc99KolmzZvrLX/6iRYsW6Y033tCoUaM0f/58LVu2TKGhoVq8eLHefvttDR8+/MoaBKDaIiQB8KkvvvhCffv2VUhIiCQpPDz8suaLi4vTpEmTFBgYqD59+igmJqbca/7xj3/o6NGjeueddyRJ586dU2ZmpiTp3nvvvWRAkqRdu3Zp5cqVKioqUnZ2to4ePeoOSb1793b/d968eTp27JiOHDmihIQESRdu3dG9e/fLei8AaiZCEgCf8nZR/8DAQJWWlrofO51OSVKnTp305ptvavv27Xrqqaf0wgsvuO+vdfGyZ86cqdtuu63ccq+77rpLjltUVKQ5c+bov//7v9W4cWPNmTPHPfbPBQQEyLIstWvXTqtWrbrkcgHUHhyTBMCnOnfurPfee899XM/Pd7fdeOONOnr0qFwulzIyMpSamipJOnnypBo3bqyhQ4fq/vvvV2pqqurXr6+zZ8+6573zzju1bt06d8g6fPjwZddVVFQkh8Oh8PBw5efn66OPPvKY/v7777v/27FjR7Vs2VLfffedu77CwkKdOHGics0AUKOwJQmAT7Vp00aDBg3SoEGDFBgYqN69e+upp55yT2/WrJnuvPNOPfDAA2rbtq3+7d/+TdKF3XSvv/666tSpo/DwcC1cuFCRkZFq1aqVHnjgAfXv31+PPvqoTpw4of79+6u0tFRt27bV/PnzL6uu0NBQ/eY3v1G/fv3UrFkz9262MllZWXrwwQcVGRmpV155RXXr1tX8+fP1wgsv6Ny5cwoICNAf/vCHclu3ANQe3OAWAADAgN1tAAAABoQkAAAAA0ISAACAASEJAADAgJAEAABgQEgCAAAwICQBAAAYEJIAAAAM/j9KK3N1MYIflgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 7), constrained_layout=True)\n",
    "\n",
    "sns.boxplot(x=\"cluster\", y=\"error\", hue='cluster', dodge=False, data=labels, showfliers=False, linewidth=.75, ax=ax)\n",
    "plt.legend([],[], frameon=False)\n",
    "sns.despine()\n",
    "ax.set_ylabel(\"within sum of squares\")\n",
    "ax.set_xlabel(\"cluster label\")\n",
    "\n",
    "x_labels = [cluster_names[i] + f\"\\n({obs_counts[i]})\" for i, item in enumerate(ax.get_xticklabels())]\n",
    "ax.set_xticklabels(x_labels)\n",
    "\n",
    "plt.savefig(f\"figs/squared_error_sum.pdf\", dpi=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "healthy-breakdown",
   "metadata": {},
   "outputs": [],
   "source": [
    "sector = labels.PCDS.str.split(\" \").apply(lambda x: x[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "median-civilization",
   "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>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PCDS</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>AB10</th>\n",
       "      <td>675.0</td>\n",
       "      <td>0.201004</td>\n",
       "      <td>0.388954</td>\n",
       "      <td>0.023556</td>\n",
       "      <td>0.081900</td>\n",
       "      <td>0.123259</td>\n",
       "      <td>0.183171</td>\n",
       "      <td>3.247271</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AB11</th>\n",
       "      <td>642.0</td>\n",
       "      <td>0.260221</td>\n",
       "      <td>0.372207</td>\n",
       "      <td>0.045780</td>\n",
       "      <td>0.132138</td>\n",
       "      <td>0.186119</td>\n",
       "      <td>0.251538</td>\n",
       "      <td>3.485213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AB12</th>\n",
       "      <td>797.0</td>\n",
       "      <td>0.597891</td>\n",
       "      <td>0.831370</td>\n",
       "      <td>0.041540</td>\n",
       "      <td>0.150103</td>\n",
       "      <td>0.260930</td>\n",
       "      <td>0.607934</td>\n",
       "      <td>4.599404</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AB13</th>\n",
       "      <td>86.0</td>\n",
       "      <td>0.627698</td>\n",
       "      <td>0.579101</td>\n",
       "      <td>0.059035</td>\n",
       "      <td>0.268618</td>\n",
       "      <td>0.452741</td>\n",
       "      <td>0.683968</td>\n",
       "      <td>3.254864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AB14</th>\n",
       "      <td>142.0</td>\n",
       "      <td>0.443040</td>\n",
       "      <td>0.710041</td>\n",
       "      <td>0.051355</td>\n",
       "      <td>0.105250</td>\n",
       "      <td>0.153745</td>\n",
       "      <td>0.313562</td>\n",
       "      <td>3.430795</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",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>YO8</th>\n",
       "      <td>1294.0</td>\n",
       "      <td>0.504283</td>\n",
       "      <td>0.715602</td>\n",
       "      <td>0.037709</td>\n",
       "      <td>0.117974</td>\n",
       "      <td>0.208944</td>\n",
       "      <td>0.527891</td>\n",
       "      <td>6.540000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>YO90</th>\n",
       "      <td>10.0</td>\n",
       "      <td>0.574485</td>\n",
       "      <td>0.355448</td>\n",
       "      <td>0.184361</td>\n",
       "      <td>0.184361</td>\n",
       "      <td>0.706791</td>\n",
       "      <td>0.898455</td>\n",
       "      <td>0.898455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZE1</th>\n",
       "      <td>327.0</td>\n",
       "      <td>0.576843</td>\n",
       "      <td>0.690216</td>\n",
       "      <td>0.061873</td>\n",
       "      <td>0.255441</td>\n",
       "      <td>0.325530</td>\n",
       "      <td>0.557516</td>\n",
       "      <td>4.791093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZE2</th>\n",
       "      <td>311.0</td>\n",
       "      <td>1.552007</td>\n",
       "      <td>0.971910</td>\n",
       "      <td>0.301914</td>\n",
       "      <td>0.827429</td>\n",
       "      <td>1.272075</td>\n",
       "      <td>2.016887</td>\n",
       "      <td>6.439161</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZE3</th>\n",
       "      <td>12.0</td>\n",
       "      <td>1.961530</td>\n",
       "      <td>1.760545</td>\n",
       "      <td>0.550234</td>\n",
       "      <td>0.932386</td>\n",
       "      <td>1.116277</td>\n",
       "      <td>2.107792</td>\n",
       "      <td>5.773155</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2871 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       count      mean       std       min       25%       50%       75%  \\\n",
       "PCDS                                                                       \n",
       "AB10   675.0  0.201004  0.388954  0.023556  0.081900  0.123259  0.183171   \n",
       "AB11   642.0  0.260221  0.372207  0.045780  0.132138  0.186119  0.251538   \n",
       "AB12   797.0  0.597891  0.831370  0.041540  0.150103  0.260930  0.607934   \n",
       "AB13    86.0  0.627698  0.579101  0.059035  0.268618  0.452741  0.683968   \n",
       "AB14   142.0  0.443040  0.710041  0.051355  0.105250  0.153745  0.313562   \n",
       "...      ...       ...       ...       ...       ...       ...       ...   \n",
       "YO8   1294.0  0.504283  0.715602  0.037709  0.117974  0.208944  0.527891   \n",
       "YO90    10.0  0.574485  0.355448  0.184361  0.184361  0.706791  0.898455   \n",
       "ZE1    327.0  0.576843  0.690216  0.061873  0.255441  0.325530  0.557516   \n",
       "ZE2    311.0  1.552007  0.971910  0.301914  0.827429  1.272075  2.016887   \n",
       "ZE3     12.0  1.961530  1.760545  0.550234  0.932386  1.116277  2.107792   \n",
       "\n",
       "           max  \n",
       "PCDS            \n",
       "AB10  3.247271  \n",
       "AB11  3.485213  \n",
       "AB12  4.599404  \n",
       "AB13  3.254864  \n",
       "AB14  3.430795  \n",
       "...        ...  \n",
       "YO8   6.540000  \n",
       "YO90  0.898455  \n",
       "ZE1   4.791093  \n",
       "ZE2   6.439161  \n",
       "ZE3   5.773155  \n",
       "\n",
       "[2871 rows x 8 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sector_data = squared_error_sum.groupby(sector).describe()\n",
    "sector_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "talented-liability",
   "metadata": {},
   "outputs": [],
   "source": [
    "sector_data.to_csv(\"../urbangrammar_samba/postcode_classification/data/fit_per_sector.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IMD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "regulation-copying",
   "metadata": {},
   "outputs": [],
   "source": [
    "onspd = pd.read_csv(\"../urbangrammar_samba/postcode_classification/data/ONSPD_202102_utm.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "indie-appraisal",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['pcds', 'usertype', 'oseast1m', 'osnrth1m', 'osgrdind', 'lat', 'long',\n",
       "       'ctry', 'current', 'dointr', 'doterm', 'doupdt', 'oslaua', 'osward',\n",
       "       'oscty', 'pcon', 'oa11', 'lsoa11', 'msoa11', 'parish', 'imd',\n",
       "       'UTM_ZONE', 'UTM_X', 'UTM_Y'],\n",
       "      dtype='object')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "onspd.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89ea0085-b821-49f1-9bb5-ef80c26b94e6",
   "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>PCDS</th>\n",
       "      <th>geometry</th>\n",
       "      <th>cluster</th>\n",
       "      <th>error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AB10 1AB</td>\n",
       "      <td>POINT (394235.441 806530.784)</td>\n",
       "      <td>0</td>\n",
       "      <td>0.562764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AB10 1AF</td>\n",
       "      <td>POINT (394181.440 806430.784)</td>\n",
       "      <td>0</td>\n",
       "      <td>0.607934</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AB10 1AG</td>\n",
       "      <td>POINT (394230.026 806469.006)</td>\n",
       "      <td>2</td>\n",
       "      <td>0.535169</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AB10 1AH</td>\n",
       "      <td>POINT (394371.467 806360.721)</td>\n",
       "      <td>2</td>\n",
       "      <td>0.189136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AB10 1AL</td>\n",
       "      <td>POINT (394296.014 806581.022)</td>\n",
       "      <td>0</td>\n",
       "      <td>0.615652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710710</th>\n",
       "      <td>ZE3 9JU</td>\n",
       "      <td>POINT (438872.230 1111806.351)</td>\n",
       "      <td>4</td>\n",
       "      <td>0.958750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710711</th>\n",
       "      <td>ZE3 9JW</td>\n",
       "      <td>POINT (438975.215 1110037.332)</td>\n",
       "      <td>3</td>\n",
       "      <td>1.116277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710712</th>\n",
       "      <td>ZE3 9JX</td>\n",
       "      <td>POINT (438872.255 1110218.355)</td>\n",
       "      <td>3</td>\n",
       "      <td>1.515320</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710713</th>\n",
       "      <td>ZE3 9JY</td>\n",
       "      <td>POINT (438498.246 1112028.331)</td>\n",
       "      <td>4</td>\n",
       "      <td>1.611978</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710714</th>\n",
       "      <td>ZE3 9JZ</td>\n",
       "      <td>POINT (438662.260 1112121.367)</td>\n",
       "      <td>4</td>\n",
       "      <td>0.853294</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1710715 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             PCDS                        geometry  cluster     error\n",
       "0        AB10 1AB   POINT (394235.441 806530.784)        0  0.562764\n",
       "1        AB10 1AF   POINT (394181.440 806430.784)        0  0.607934\n",
       "2        AB10 1AG   POINT (394230.026 806469.006)        2  0.535169\n",
       "3        AB10 1AH   POINT (394371.467 806360.721)        2  0.189136\n",
       "4        AB10 1AL   POINT (394296.014 806581.022)        0  0.615652\n",
       "...           ...                             ...      ...       ...\n",
       "1710710   ZE3 9JU  POINT (438872.230 1111806.351)        4  0.958750\n",
       "1710711   ZE3 9JW  POINT (438975.215 1110037.332)        3  1.116277\n",
       "1710712   ZE3 9JX  POINT (438872.255 1110218.355)        3  1.515320\n",
       "1710713   ZE3 9JY  POINT (438498.246 1112028.331)        4  1.611978\n",
       "1710714   ZE3 9JZ  POINT (438662.260 1112121.367)        4  0.853294\n",
       "\n",
       "[1710715 rows x 4 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "geom = gpd.read_parquet(\"postcodes_clustered.pq\")\n",
    "geom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "corporate-plant",
   "metadata": {},
   "outputs": [],
   "source": [
    "imd = geom.drop(columns=\"error\").merge(onspd[[\"pcds\", \"imd\"]], left_on=\"PCDS\", right_on=\"pcds\", how=\"inner\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "advisory-russell",
   "metadata": {},
   "outputs": [],
   "source": [
    "imd['cluster'] = imd.cluster.replace({0: \"A\", 1: \"B\", 2: \"C\", 3: \"D\", 4: \"E\", 5: \"F\", 6: \"G\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "played-happening",
   "metadata": {},
   "outputs": [],
   "source": [
    "imd.imd = imd.imd.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "romantic-booth",
   "metadata": {},
   "outputs": [],
   "source": [
    "imd['cluster'] = pd.Categorical(imd.cluster, ordered=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0df1708b-81fd-4aa6-93ca-7eda322b4a56",
   "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>PCDS</th>\n",
       "      <th>geometry</th>\n",
       "      <th>cluster</th>\n",
       "      <th>pcds</th>\n",
       "      <th>imd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AB10 1AB</td>\n",
       "      <td>POINT (394235.441 806530.784)</td>\n",
       "      <td>A</td>\n",
       "      <td>AB10 1AB</td>\n",
       "      <td>4907</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AB10 1AF</td>\n",
       "      <td>POINT (394181.440 806430.784)</td>\n",
       "      <td>A</td>\n",
       "      <td>AB10 1AF</td>\n",
       "      <td>4907</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AB10 1AG</td>\n",
       "      <td>POINT (394230.026 806469.006)</td>\n",
       "      <td>C</td>\n",
       "      <td>AB10 1AG</td>\n",
       "      <td>4907</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AB10 1AH</td>\n",
       "      <td>POINT (394371.467 806360.721)</td>\n",
       "      <td>C</td>\n",
       "      <td>AB10 1AH</td>\n",
       "      <td>3290</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AB10 1AL</td>\n",
       "      <td>POINT (394296.014 806581.022)</td>\n",
       "      <td>A</td>\n",
       "      <td>AB10 1AL</td>\n",
       "      <td>4907</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710710</th>\n",
       "      <td>ZE3 9JU</td>\n",
       "      <td>POINT (438872.230 1111806.351)</td>\n",
       "      <td>E</td>\n",
       "      <td>ZE3 9JU</td>\n",
       "      <td>4187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710711</th>\n",
       "      <td>ZE3 9JW</td>\n",
       "      <td>POINT (438975.215 1110037.332)</td>\n",
       "      <td>D</td>\n",
       "      <td>ZE3 9JW</td>\n",
       "      <td>4187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710712</th>\n",
       "      <td>ZE3 9JX</td>\n",
       "      <td>POINT (438872.255 1110218.355)</td>\n",
       "      <td>D</td>\n",
       "      <td>ZE3 9JX</td>\n",
       "      <td>4187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710713</th>\n",
       "      <td>ZE3 9JY</td>\n",
       "      <td>POINT (438498.246 1112028.331)</td>\n",
       "      <td>E</td>\n",
       "      <td>ZE3 9JY</td>\n",
       "      <td>4187</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710714</th>\n",
       "      <td>ZE3 9JZ</td>\n",
       "      <td>POINT (438662.260 1112121.367)</td>\n",
       "      <td>E</td>\n",
       "      <td>ZE3 9JZ</td>\n",
       "      <td>4187</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1710715 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             PCDS                        geometry cluster      pcds   imd\n",
       "0        AB10 1AB   POINT (394235.441 806530.784)       A  AB10 1AB  4907\n",
       "1        AB10 1AF   POINT (394181.440 806430.784)       A  AB10 1AF  4907\n",
       "2        AB10 1AG   POINT (394230.026 806469.006)       C  AB10 1AG  4907\n",
       "3        AB10 1AH   POINT (394371.467 806360.721)       C  AB10 1AH  3290\n",
       "4        AB10 1AL   POINT (394296.014 806581.022)       A  AB10 1AL  4907\n",
       "...           ...                             ...     ...       ...   ...\n",
       "1710710   ZE3 9JU  POINT (438872.230 1111806.351)       E   ZE3 9JU  4187\n",
       "1710711   ZE3 9JW  POINT (438975.215 1110037.332)       D   ZE3 9JW  4187\n",
       "1710712   ZE3 9JX  POINT (438872.255 1110218.355)       D   ZE3 9JX  4187\n",
       "1710713   ZE3 9JY  POINT (438498.246 1112028.331)       E   ZE3 9JY  4187\n",
       "1710714   ZE3 9JZ  POINT (438662.260 1112121.367)       E   ZE3 9JZ  4187\n",
       "\n",
       "[1710715 rows x 5 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "imd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "general-looking",
   "metadata": {},
   "outputs": [],
   "source": [
    "england = gpd.read_file(\"../urbangrammar_samba/postcode_classification/data/Regions_(December_2017)_Boundaries.geojson\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "general-above",
   "metadata": {},
   "outputs": [],
   "source": [
    "england = england.to_crs(27700).unary_union"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "engaged-handle",
   "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>PCDS</th>\n",
       "      <th>geometry</th>\n",
       "      <th>cluster</th>\n",
       "      <th>pcds</th>\n",
       "      <th>imd</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1554107</th>\n",
       "      <td>TR21 0HE</td>\n",
       "      <td>POINT (90785.972 10170.990)</td>\n",
       "      <td>D</td>\n",
       "      <td>TR21 0HE</td>\n",
       "      <td>22165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1554147</th>\n",
       "      <td>TR21 0NH</td>\n",
       "      <td>POINT (91410.724 10209.391)</td>\n",
       "      <td>E</td>\n",
       "      <td>TR21 0NH</td>\n",
       "      <td>22165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1554148</th>\n",
       "      <td>TR21 0NJ</td>\n",
       "      <td>POINT (91450.760 10239.418)</td>\n",
       "      <td>E</td>\n",
       "      <td>TR21 0NJ</td>\n",
       "      <td>22165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1554161</th>\n",
       "      <td>TR21 0PA</td>\n",
       "      <td>POINT (91339.761 10262.423)</td>\n",
       "      <td>E</td>\n",
       "      <td>TR21 0PA</td>\n",
       "      <td>22165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1554149</th>\n",
       "      <td>TR21 0NL</td>\n",
       "      <td>POINT (91529.731 10309.366)</td>\n",
       "      <td>E</td>\n",
       "      <td>TR21 0NL</td>\n",
       "      <td>22165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1559379</th>\n",
       "      <td>TS10 3PT</td>\n",
       "      <td>POINT (461759.026 523988.577)</td>\n",
       "      <td>G</td>\n",
       "      <td>TS10 3PT</td>\n",
       "      <td>18345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1559378</th>\n",
       "      <td>TS10 3PS</td>\n",
       "      <td>POINT (461864.991 524001.589)</td>\n",
       "      <td>G</td>\n",
       "      <td>TS10 3PS</td>\n",
       "      <td>18345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1559375</th>\n",
       "      <td>TS10 3PP</td>\n",
       "      <td>POINT (461898.033 524032.649)</td>\n",
       "      <td>G</td>\n",
       "      <td>TS10 3PP</td>\n",
       "      <td>18345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1559363</th>\n",
       "      <td>TS10 3NW</td>\n",
       "      <td>POINT (461816.001 524074.607)</td>\n",
       "      <td>G</td>\n",
       "      <td>TS10 3NW</td>\n",
       "      <td>18345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1559358</th>\n",
       "      <td>TS10 3NP</td>\n",
       "      <td>POINT (461781.021 524081.585)</td>\n",
       "      <td>G</td>\n",
       "      <td>TS10 3NP</td>\n",
       "      <td>18345</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1458467 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             PCDS                       geometry cluster      pcds    imd\n",
       "1554107  TR21 0HE    POINT (90785.972 10170.990)       D  TR21 0HE  22165\n",
       "1554147  TR21 0NH    POINT (91410.724 10209.391)       E  TR21 0NH  22165\n",
       "1554148  TR21 0NJ    POINT (91450.760 10239.418)       E  TR21 0NJ  22165\n",
       "1554161  TR21 0PA    POINT (91339.761 10262.423)       E  TR21 0PA  22165\n",
       "1554149  TR21 0NL    POINT (91529.731 10309.366)       E  TR21 0NL  22165\n",
       "...           ...                            ...     ...       ...    ...\n",
       "1559379  TS10 3PT  POINT (461759.026 523988.577)       G  TS10 3PT  18345\n",
       "1559378  TS10 3PS  POINT (461864.991 524001.589)       G  TS10 3PS  18345\n",
       "1559375  TS10 3PP  POINT (461898.033 524032.649)       G  TS10 3PP  18345\n",
       "1559363  TS10 3NW  POINT (461816.001 524074.607)       G  TS10 3NW  18345\n",
       "1559358  TS10 3NP  POINT (461781.021 524081.585)       G  TS10 3NP  18345\n",
       "\n",
       "[1458467 rows x 5 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "imd_england = imd.iloc[imd.sindex.query(england, \"contains\")]\n",
    "imd_england"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "continuing-subcommittee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAIBCAYAAABdiys1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzV0lEQVR4nO3de1hVdaL/8c8GFATdXBQpL4dSEGtqSCDTSYFGU7s4ZuV4SW0mL9PvJFB6FOziOFjhLQvwjKmpTaZoZpeni9WZLlhPOoma3TRGUx/LGkTkbiDs/fuj0z5DfHGjtfcCfL+exyf32mut72evkj5+19pr2ZxOp1MAAABowMfqAAAAAC0RJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAwumJIUExNjdQQAANCKXDAlCQAA4FxQkgAAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYEBJAgAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADP6sDtHarVq3SoUOHzmvb48ePq1u3bue1be/evTV9+vTz2tYKHKfmO99jxXFqvgvpWHGcms+KY8Vxaj5LjpXzAtGnTx+rIzQyZ84cqyO0Chyn5uE4NR/Hqnk4Ts3HsWqe1nacON0GAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYEBJAgAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMDAz1M7/te//qUZM2bI399fdXV1mj9/viIjI5WRkaETJ04oKipK8+fPl4+Pj/bt26esrCw5nU5Nnz5dQ4YMkSRlZ2dr586dCgoK0uLFixUWFqaSkhLNmTNHVVVVGjBggNLS0jz1EQAAwAXMYzNJXbp00ebNm/Xss88qLS1NTz31lLZu3arY2Fht3LhRfn5+2r59uyQpKytL2dnZWrt2rZYvX666ujoVFhZq//79ysvL0+23367Vq1dLklavXq0xY8YoLy9P+/fvV2Fhoac+AgAAuIB5rCT5+vrKx+eH3VdVVenyyy9XQUGBkpKSJElJSUkqKChQTU2NHA6HIiIiFBQUpMjISB09erTBusnJydqzZ48kaffu3a7liYmJ2r17t6c+AgAAuIB59JqkgwcPaty4cVqwYIESEhJUXl4uu90uSQoODlZZWZlKS0vVqVMn1zZ2u11lZWUN1g0ICFB1dbUk6fTp0woICGiwj5/Kzc1VTExMg18AAADnwmPXJElSVFSUNm3apP3792vevHnq3r27KioqFB4ervLycgUHByskJESVlZWubSoqKhQcHCy73a6KigpJUk1NjQIDAyVJHTp0UE1Njfz9/V37+KmUlBSlpKQ0WEZRAgAA58JjM0m1tbWu39vtdgUEBOjqq69Wfn6+JGn79u1KSEiQv7+/fHx8VFRUpOrqah05ckSRkZEN1s3Pz1dcXJwkKT4+vsHyhIQET30EAABwAfPYTNInn3yiJ554QjabTZKUkZGhXr16KSMjQxMmTFBUVJQSExNd76WmpsrhcGjGjBny8/NTdHS0+vbtq/HjxyswMFBLliyRJE2bNk2zZ8/WunXrNGDAAEVHR3vqIwAAgAuYx0pSQkKCnn322UbLs7OzGy2LjY3Vpk2bGi1PS0tr9BX/sLAwrVmz5pcLCgAAYMDNJAEAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYEBJAgAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAz8rA6A1uPP96errKTY6+MWnSrXvXdP8fq4wWFd9JdHF3l9XABAy0BJQrOVlRRrUi+bBSMHWzCmtP4r7xdCAEDLwek2AAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADDwszoAAODs7p/7kE6ePOX1cU+VFutP02d4fdzOnUP1aNYCr48L/BQlCQBauJMnTykyPNnr40aGe31ISdLRE+9ZMzDwE5xuAwAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMDAz1M7PnTokB588EH5+PjIx8dHjz76qHbt2qW//vWvuvjiiyVJq1evVkBAgPbt26esrCw5nU5Nnz5dQ4YMkSRlZ2dr586dCgoK0uLFixUWFqaSkhLNmTNHVVVVGjBggNLS0jz1EQAAaJPmPzBPZSVlXh+3qOSE7vt/3v//dnBYsOY/knnO23msJIWGhmrlypWy2+3avn27nnzyScXHx2v8+PGaMmVKg3WzsrKUnZ2tjh07auLEiUpKStJXX32l/fv3Ky8vT2+88YZWr16t9PR0rV69WmPGjNHw4cN19913q7CwUH369PHUxwAAoM0pKynTH64cbXUMr3n60xfPazuPnW4LCwuT3W6XJPn5+cnX11eStGXLFk2YMEFr166VJNXU1MjhcCgiIkJBQUGKjIzU0aNHVVBQoKSkJElScnKy9uzZI0navXu3a3liYqJ2797daOzc3FzFxMQ0+AUAAHAuPDaT9KPTp08rJydHjzzyiMLDwzVq1Cg5HA6lpaXpsssuU69evdSpUyfX+na7XWVlZSovL1fPnj0lSQEBAaqurnbtLyAgQJIUHBysY8eONRozJSVFKSkpDZZRlAAAwLnw6IXbdXV1uu+++zR16lT17t1bdrtdvr6+ateunYYNG6bPP/9cISEhqqysdG1TUVGh4OBg2e12VVRUSPphtikwMFCS1KFDB9XU1EiSysvLFRwc7MmPAAAALlAeK0lOp1MPPPCAEhMTNXToUElylR5J+uijjxQZGSl/f3/5+PioqKhI1dXVOnLkiCIjI3X11VcrPz9fkpSfn6+4uDhJUnx8fIPlCQkJnvoIAADgAuax023vv/++3njjDR0/flzbtm1T37591bFjR33wwQfy9fXV5Zdf7ipPGRkZSk1NlcPh0IwZM+Tn56fo6Gj17dtX48ePV2BgoJYsWSJJmjZtmmbPnq1169ZpwIABio6O9tRHAAAAFzCPlaTExETt27ev0XLTV/ZjY2O1adMm47o/XT8sLExr1qz55YICAAAYcDNJAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAwoSQAAAAaUJAAAAAM/qwMAuHDNS0/XqRPFXh+3uLxcKXdN8fq4oeFdlLlokdfHBXB+KEkALHPqRLFuqHN4f+DAjpIF426zoBACOH+cbgMAADCgJAEAABhQkgAAAAyadU3SV199peLiYvn7++vSSy+V3W73dC4AAABLNVmSjhw5orVr1+of//iHevToodDQUNXW1urYsWMKCAjQbbfdptGjR8vX19ebeQEAALyiyZK0cuVKjRs3TpmZmY3eKysr07Zt2/T6669r5MiRHg0IAABghSZLUlZWliSpvr6+0WyRr6+vxo0b59lkAAAAFnJ74fZ9992n+vp61+vi4mLdddddHg0FAABgNbclafDgwUpLS1NdXZ2OHTumKVOmKC0tzRvZAAAALOP2221jxoyRzWbTn/70J504cUKPPvqorrjiCm9kAwAAsEyTJemxxx6TzWaTJDmdThUVFemSSy7RW2+9pbfeekszZ870WkgAAABva7Ik9erV66yvAQAA2rImS9Lo0aO9mQMAAKBFcXtN0ocffqjc3Fx9++23qqurcy3/4IMPPBoMAADASm5L0vz58/X444/r8ssvd12jBAAA0Na5LUldunTRr371K29kAQAAaDHclqTevXtrypQpSk5OVvv27V3Lx44d69FgAAAAVnJbkiIiIhQREaGysjJv5AEAAGgR3JakGTNmeCMHAABAi+K2JP3rX//SihUrdOjQIZ05c8a1fNOmTR4NBgAAYCW3z26bPXu2rr32Wp06dUqZmZnq16+frrnmGm9kAwAAsIzbklRZWanrr79eNptNffr0UXp6unbt2uWNbAAAAJZxe7qtffv2cjqd6tmzp1544QVFRESopKTEG9kAAAAs47YkzZ07V9XV1XrooYeUnZ2tqqoqLVy40BvZAAAALHPWklRfX6+///3vio2NVVBQEOUIAABcMM56TZKvr68++eQTORwOb+UBAABoEdyebrvkkkt05513asiQIerQoYNrOXfcBgAAbZnbkhQeHq7w8HBVVlaqsrLSG5kAADgv989NV8nJk14f91Rpme6ePtWrY4Z17qxHsxZ5dcwLDXfcBoAW7lRpscpKX7c6htc4VH3e25acPKk+nX1/wTTN1DnM60MWWlAGLzRNlqQ5c+bo97//vRISEhq9d+rUKb366quy2+0aNWqURwMCwIUuNKSLIsOTrY7hNUdPvGd1BEDSWUpSamqqnn76ac2bN09dunRRaGioamtr9c033yg0NFRjxozRjTfe6M2ssFjRqXKt2XvG/YptRKmjndURAAAWarIk9ejRQw8++KAk6dixYyouLpa/v78uueQSBQYGei0gWo6uoXZN6mWzOobXrP/KaXUEAICF3F6TJEk9e/ZUz549PZ0FAACgxXD77DYAAIALESUJAADAoFmn26QfHlFSX1/vet2+fXuPBAIAAGgJ3Jakl156STk5OZIkm80mp9Mpm82mt99+2+PhAAAArOK2JP31r3/Vc889py5dupzTjg8dOqQHH3xQPj4+8vHx0aOPPqouXbooIyNDJ06cUFRUlObPny8fHx/t27dPWVlZcjqdmj59uoYMGSJJys7O1s6dOxUUFKTFixcrLCxMJSUlmjNnjqqqqjRgwAClpaWd3ycHAAA4C7fXJHXr1k12u/2cdxwaGqqVK1dqw4YNmjZtmp588klt3bpVsbGx2rhxo/z8/LR9+3ZJUlZWlrKzs7V27VotX75cdXV1Kiws1P79+5WXl6fbb79dq1evliStXr1aY8aMUV5envbv36/CwsJzzgYAAOCO25mksLAw3X777frNb37T4DqkmTNnut3ONYifn3x9fVVQUKCUlBRJUlJSknbt2qWBAwfK4XAoIiJCkhQZGamjR4+qoKBASUlJkqTk5GStW7dOkrR7927X7FFiYqJ2796tPn36nMtnBgAAcMttSRo8eLAGDx583gOcPn1aOTk5euSRR/TII4+4ZqWCg4NVVlam0tJSderUybW+3W5XWVmZysvLXfdmCggIUHV1tWt/AQEBrn0cO3as0Zi5ublavnz5eWcGAABwW5JGjx6t06dP65///KckKTo6Wh06dGjWzuvq6nTfffdp6tSp6t27t+x2uyoqKhQeHq7y8nIFBwcrJCRElZWVrm0qKioUHBzsWleSampqXHf57tChg2pqauTv7+/ax0+lpKS4Zqx+FBMT06zMAAAAUjOuSXr33Xd100036amnntKqVas0cuRI5efnu92x0+nUAw88oMTERA0dOlSSdPXVV7u23b59uxISEuTv7y8fHx8VFRWpurpaR44cUWRkZIN18/PzFRcXJ0mKj49vsNz0AF4AAICfy+1M0hNPPKHNmzcrPDxcknTixAlNnTrVdb1QU95//3298cYbOn78uLZt26a+fftq5syZysjI0IQJExQVFaXExERJUkZGhlJTU+VwODRjxgz5+fkpOjpaffv21fjx4xUYGKglS5ZIkqZNm6bZs2dr3bp1GjBggKKjo3/uMQAAAGjEbUlyOBwNLsIODQ2Vw+Fwu+PExETt27ev0fLs7OxGy2JjY7Vp06ZGy9PS0hp9xT8sLExr1qxxOz4AAMDP4bYk3XDDDbrjjjt0/fXXy2az6a233tJNN93kjWwAAACWcVuS/vM//1ODBg3S3r175XQ6df/99+vXv/61N7IBAABYpsmSdPz4cXXr1k2HDx9Wp06dXNcPSdLhw4d16aWXeiUgAACAFZosSStWrNCCBQs0b9481zPbfmSz2fTMM894JSAAAIAVmixJCxYskCRlZmY2mjU6fPiwZ1MBAABYzO19kkyPH3H3SBIAAIDWrsmZpM8//1yffvqpSktLtXnzZtfyyspK1dXVeSUcAACAVZosSVVVVSouLtaZM2d04sQJ1/LAwEDjvY4AAADakiZLUv/+/dW/f3/9/ve/V9euXb2ZCQAAwHJu75PkdDo1f/58HTp0SGfOnHEtN90hGwAAoK1we+H27Nmzde211+rUqVPKzMxUv379dM0113gjGwAAgGXcziRVVlbq+uuvV05Ojvr06aP09HRNmDDBG9kAtHHF5eV6obbW6hheU9W+vdURAJwDtyWpffv2cjqd6tmzp1544QVFRESopKTEG9kAtHFd7HbdUOf+gdltxTY/t5P3AFoQtyVp7ty5qq6u1kMPPaTs7GxVVVVp4cKF3sgGAABgGbclSZKCgoIUFBREOQIAABcMtyXpqaee0j//+U8lJiZq+PDhio+P90YuAAAAS7k9QZ6bm6sXX3xRcXFx2rBhg0aMGKGHH37YG9kAAAAs06zTbR06dFBSUpIcDoeqq6v15ptv6sEHH/R0Nq+Zk/Ggik+e8vq45aUndde0e7w+bpfOoVq8kKILAMDZuC1JL7/8st566y0dPHhQgwcP1pQpU5SQkOCNbF5TfPKUvg/9jdfHbR8qfe/1UaXikx9aMCoAAK2L25L0xRdfaMqUKYqLi/NGHgAAgBahyZJ0/PhxdevWTePGjZMkHT58uMH7l156qWeTAQAAWKjJkrRixQotWLBA8+bNk81mk9PpdL1ns9n0zDPPeCUgAACAFZosSQsWLFB9fb0mTpyo4cOHezMTAACA5c56CwBfX19mjAAAwAXJ7YXbCQkJevzxxzV8+HB16NDBtZxrkgAAQFvmtiTt2bOnwT8lrkkCAABtn9uStH79em/kAAAAXlJUckKrPsyzOobXlJ6pOK/t3JakoqIiZWdna+/evbLZbOrXr59SU1PVtWvX8xoQAABYq2tYuP5w5WirY3jN05++eF7buX1227333qv4+Hht2bJFzz33nPr166e0tLTzGgwAAKC1cFuSqqqqdOuttyooKEhBQUG67bbbVF9f741sAAAAlnF7um3AgAFatmyZhg0bJpvNpjfffFMDBw503YGbb7kBAIC2qFnPbpOkvXv3Nli+Z88evuUGAADaLL7dBgAAYOC2JH399ddavHixiouLtXHjRh08eFA7duzQpEmTvJEPaJXSH8xQ8amTXh2z/GSpptwzzatjSlKX0M5a9PBCr48LAJ7mtiTNnTtXaWlpevjhhyVJvXv31r333ktJAs6i+NRJ1ScEenXMIAXKiq9UFBd4twwCgLe4/Xbb6dOnlZCQ4Hpts9nk6+vr0VAAAABWc1uSwsPDdeDAAdlsNknSli1b1KNHD48HAwAAsJLbkpSZmalVq1apqKhIgwcP1ocffqjMzExvZAMAALCM22uSwsPDtWzZMm9kAQAAaDGaLElz584964ZZWVm/eBgAAICWosnTbTfeeKNuvPFGBQUFqbq6WklJSUpMTFR1dbXCwsK8mREAAMDrmpxJGjx4sCQpOztbzz//vGv5DTfcoAkTJng+GQAAgIXcXrj9/fffq7Cw0PW6sLBQ5eXlHg0FAABgNbcXbi9YsECzZs2SJDmdTvn4+LhuLAkAANBWuS1J/fr10yuvvKKKigo5nU7Z7XZv5AIAALCU25L0o06dOnkyBwAAQIvi9pokAACACxElCQAAwMBtSfr666+Vmprq+tr/wYMHtX79eo8HAwAAsJLba5Lmzp2rtLQ01zfaevfurXvvvVeTJk3yeDgAgNS5c6iOnnjP6+OeKi1WaEgXr4/buXOo18cETNyWpNOnTyshIcH12mazydfX16OhAAD/59GsBZaMm56erkWLFlkyNtASuD3dFh4ergMHDshms0mStmzZoh49eng8GAAAgJXclqTMzEytWrVKRUVFGjx4sD788ENlZmZ6IxsAAIBl3J5uCw8P17Jly7yRBQAAoMVosiTNnDnTdYrN5LHHHvNIIAAAgJagyZI0btw4b+YAAABoUZosSf3793f9/ptvvtG+fftks9n061//Wt27d/dKOAAAAKu4vXD76aef1t13361PP/1Un3zyie6++2797W9/c7vj2tpajRs3TgkJCXrjjTckSS+88IKGDh2qSZMmadKkSfr+++8lSfv27dO4ceM0duxYvf322659ZGdna/z48Zo6dapKSkokSSUlJZo6darGjx+v7Ozs8/rQAAAA7ri9cHvz5s166aWX5O/vL0lKS0vTrbfeqjvvvPPsO/bzU05OjjZv3txg+fjx4zVlypQGy7KyspSdna2OHTtq4sSJSkpK0ldffaX9+/crLy9Pb7zxhlavXq309HStXr1aY8aM0fDhw3X33XersLBQffr0OdfPDQAAcFZuZ5K6deumqqoq1+vq6mr16tXL/Y59fNS1a9dGy7ds2aIJEyZo7dq1kqSamho5HA5FREQoKChIkZGROnr0qAoKCpSUlCRJSk5O1p49eyRJu3fvdi1PTEzU7t27m/ExAQAAzo3bmSSn06nhw4crPj5ekrRnzx5dddVVmjVrlqRz+5bb0KFDNWrUKDkcDqWlpemyyy5Tr1691KlTJ9c6drtdZWVlKi8vV8+ePSVJAQEBqq6ulvTDHcADAgIkScHBwTp27FijcXJzc7V8+fJm5wIAAPgptyXp7rvvbvD6rrvuOu/B7Ha7JMnX11fDhg3T559/rri4OFVWVrrWqaioUHBwsOx2uyoqKiT9MNsUGBgoSerQoYNqamrk7++v8vJyBQcHNxonJSVFKSkpDZbFxMScd24AAHDhabIk1dXVyc/PT1dddZXx/fbt25/zYBUVFa5Zo48++kjXXXed/P395ePjo6KiInXs2FFHjhxRZGSkHA6Hli1bpnHjxik/P19xcXGSpPj4eOXn52vYsGHKz893zWgBAAD8kposSffcc49WrlypESNGNLippNPplM1ma/AttKakpKToiy++UGBgoD7++GN16NBBH3zwgXx9fXX55Zdr6NChkqSMjAylpqbK4XBoxowZ8vPzU3R0tPr27avx48crMDBQS5YskSRNmzZNs2fP1rp16zRgwABFR0f/3GMAwCKh4V207USx18ctLi9Xl/+d2fam0PAuXh8TwPlrsiStXLlSkvTKK68oKCiowXv/fiH32eTm5jZalpaW1mhZbGysNm3aZFz3p+uHhYVpzZo1zRofQMuWadET5nm6PYDmcPvttokTJzZrGQAAQFvS5EzSsWPHdOTIEVVVVemDDz5wLa+srFRNTY1XwgEAAFilyZL05Zdf6u2331ZZWZlee+011/LAwEBlZmZ6JRwAAIBVmixJQ4cO1dChQ/XZZ5/piiuu8GYmAAAAyzVZkmbOnNngW20/dS43kQQAAGhtmixJ48aN82YOAACAFqXJktS/f39JUm1trdfCAAAAtBRuH0vy7zeTrKmpUUlJiS6++OJm3UwSAACgtXJbkt55550Grz/99FO9+OKLHgsEAADQEri9meRPXXnlldqzZ48nsgAAALQYbmeSHnvsMdfpNofDoQMHDug//uM/PB4MAADASm5LUq9evVy/9/HxUXJysuLj4z0aCgAAwGpuS9Lo0aO9kQMAAKBFcVuS8vPztWLFCn333Xeqq6tzLf/357kBAAC0NW5L0sMPP6zc3FzFxMSc9Q7cAAAAbYnbb7d1795dvXv3piABAIALituZpP/6r//SpEmT1K9fP7Vr1861fObMmR4NBgAAYCW3Jekvf/mL4uLiFBUVJR+fc76tEgAAQKvktiTV19dr7ty53sgCAADQYrgtSQMHDtSqVauUnJzc4HTbpZde6tFgAAAAVnJbkj755BNJ0vvvv+9aZrPZ9Mwzz3guFQAA5+FUaZkKTtW5X7ENqLe5/V84fia3R3j9+vXeyAEAwM8WGhKsPp19rY7hFYUn662O0OZxJTYAAIABJQkAAMCgyZK0du1aSdKOHTu8FgYAAKClaLIkvfjii5KkxYsXey0MAABAS9HkhdsDBgzQoEGDVFpaqkGDBjV6nwfcAgCAtqzJkvTAAw/ogQceUGpqqnJycryZCQAAwHJubwGQk5OjvXv3avfu3ZKk+Ph49evXz+PBAAAArOT2223Lly/X4sWL5e/vr/bt22vp0qX67//+b29kAwAAsIzbmaQ333xTL730knx9f7g514QJEzR69Gjdc889Hg8HAABgFbczST4+Pjp58qTr9alTp+Tjw+2VAABA2+Z2JmnmzJkaP368oqKiJEmHDh3SvHnzPB4MAADASm5LUlJSkl5//XUdPnxYTqdTvXr1kr+/vzeyAQAAWKZZjxD29/dX3759PZ0FLVxwWBet/6rY6+MWnSpX11C718cNDuvi9TEBAC1Hs0oSIEl/eXSRJeOmp6dr0SJrxgYAXLi4AhsAAMDAbUnKysrS999/73pdUlKimTNnejQUAACA1dyWpK5du2rMmDHasWOHXn75Zd1xxx0aMmSIN7IBAABYxu01SVOmTNEVV1yhu+66SyEhIXrhhRcUERHhjWwAAACWcTuTtHHjRi1YsEA5OTmaNGmSpk2bpoKCAm9kAwAAsIzbmaQDBw5o06ZN6tixo4YMGaJhw4Zp3rx5evbZZ72RDwAAwBJuZ5IyMzNVV1enzz77TJJ06aWX6umnn/Z0LgAAAEu5LUkvvfSS/vjHP+ree++V9MNjSf70pz95OhcAAICl3JakdevWafPmzerYsaMkKSoqSkVFRR4PBgAAYCW3Jaldu3Zq3769bDabJKm2ttbjoQAAAKzm9sLt6667TgsXLlRVVZVee+01Pf/88xo1apQ3sgEAAFjGbUm65557lJ+fL5vNpk8//VR/+MMflJSU5I1sAAAAlmnWA26TkpIoRgAA4ILSZEkaNGiQ6/cnT55UQECAJOn7779X586d9cEHH3g+HQAAgEWaLEk/lqDMzEwlJiYqOTlZkvTuu+9yx20AANDmuf12265du1wFSfrhQu4dO3Z4MhMAAIDl3F6TFBkZqcWLF+vGG2+UJL322mu6+OKLPR4MAADASm5L0pIlS7RhwwatWLFCTqdTcXFxSk1N9UY2AAAAy7gtSR06dNDUqVO9kQUAAKDFcFuSPvzwQ+Xm5urbb79VXV2daznfbgMAAG2Z25I0f/58Pf7447r88stdjyZpa8pLT6rm1P9YHcNram3fWx0BAIAWz21J6tKli371q1+d845ra2s1efJkHTx4UA8//LBGjBih06dPKyMjQydOnFBUVJTmz58vHx8f7du3T1lZWXI6nZo+fbqGDBkiScrOztbOnTsVFBSkxYsXKywsTCUlJZozZ46qqqo0YMAApaWlnfun/gl7SGd9H/qbn72f1iLg1IdWRwAAoMVzewuA3r17a8qUKVq/fr02b97s+uWOn5+fcnJydOedd7qWbd26VbGxsdq4caP8/Py0fft2SVJWVpays7O1du1aLV++XHV1dSosLNT+/fuVl5en22+/XatXr5YkrV69WmPGjFFeXp7279+vwsLC8/3sAAAATXI7kxQREaGIiAiVlZWd0459fHzUtWvXBssKCgqUkpIi6YdHnezatUsDBw6Uw+FQRESEpB9uOXD06FEVFBS4HoWSnJysdevWSZJ2797tmj1KTEzU7t271adPn3PKBnha+clS1b5XYnUMr2hf4/bvWgDQKrktSTNmzPjFBisvL5fdbpckBQcHq6ysTKWlperUqZNrHbvdrrKyMpWXl6tnz56SpICAAFVXV0uSTp8+7XpESnBwsI4dO9ZonNzcXC1fvvwXyw2cK3vnENUnBFodwyt8C6qtjgAAHtFkScrIyNDChQs1duxY4wXbmzZtOufB7Ha7KioqFB4ervLycgUHByskJESVlZWudSoqKhQcHOxaV5JqamoUGPjD/3A6dOigmpoa+fv7u/bxUykpKa4Zqx/FxMScc14AAHDharIkzZw5U5K0bNmyX2ywq6++Wvn5+erVq5e2b9+uQYMGyd/fXz4+PioqKlLHjh115MgRRUZGyuFwaNmyZRo3bpzy8/MVFxcnSYqPj1d+fr6GDRum/Px8zZo16xfLBwAA8KMmS9KP1xN17979vHeekpKiL774QoGBgfr444+VlpamjIwMTZgwQVFRUUpMTJT0w6xVamqqHA6HZsyYIT8/P0VHR6tv374aP368AgMDtWTJEknStGnTNHv2bK1bt04DBgxQdHT0eecDAABoittrkn6O3NzcRsuys7MbLYuNjTWevktLS2v0Ff+wsDCtWbPmlwsJAABgwNdSAAAADChJAAAABpQkAAAAA0oSAACAgUcv3AYAAC1PcFiwnv70Ra+PW1RyQl3Dwr0+bnBY43sqNgclCQCAC8z8RzItGTc9PV2LFi2yZOzzwek2AAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYEBJAgAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAICBn9UBAAD4pYR17qzCkye9Pu6p0jKFhgR7dcywzp29Ot6FiJIEAGgzHs1aZMm46enpWrTImrHhOZxuAwAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMDAz4pBr7rqKl155ZWSpMmTJ2vQoEHKyMjQiRMnFBUVpfnz58vHx0f79u1TVlaWnE6npk+friFDhkiSsrOztXPnTgUFBWnx4sUKCwuz4mMAAIA2zJKZpB49emj9+vVav369rr/+em3dulWxsbHauHGj/Pz8tH37dklSVlaWsrOztXbtWi1fvlx1dXUqLCzU/v37lZeXp9tvv12rV6+24iMAAIA2zpKS9O2332rixImaNWuWSkpKVFBQoKSkJElSUlKSCgoKVFNTI4fDoYiICAUFBSkyMlJHjx5tsG5ycrL27NnTaP+5ubmKiYlp8AsAAOBcWHK67e9//7tCQ0P1+uuva+HChSovL5fdbpckBQcHq6ysTKWlperUqZNrG7vdrrKyMpWXl6tnz56SpICAAFVXVzfaf0pKilJSUhosoygBAIBzYclMUmhoqCTphhtu0IEDB2S321VRUSFJKi8vV3BwsEJCQlRZWenapqKiQsHBwQ3WrampUWBgoPc/AAAAaPO8XpKqq6tVX18vSSooKFD37t119dVXKz8/X5K0fft2JSQkyN/fXz4+PioqKlJ1dbWOHDmiyMjIBuvm5+crLi7O2x8BAABcALx+uu2rr77SQw89pMDAQPn6+iozM1MRERHKyMjQhAkTFBUVpcTERElSRkaGUlNT5XA4NGPGDPn5+Sk6Olp9+/bV+PHjFRgYqCVLlnj7IwBudQntrOKCk14ds/xkqeydQ7w6pvTDZwWAtsjrJemKK67Qiy++2Gh5dnZ2o2WxsbHatGlTo+VpaWlKS0vzSD7gl7Do4YVeHzM9PV2LFi3y+rgA0FZZcuF2S9Olc6iKT37o9XHLS0/KHuL9v4V36Rzq9TEBAGhtKEmSFi982JJx+Zs/AAAtF48lAQAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYEBJAgAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSAACAASUJAADAwM/qAABwPlatWqVDhw6d17bHjx9Xenr6eW3bu3dvTZ8+/by2BdC6UJIAtEoUFQCexuk2AAAAA0oSAACAASUJAADAgJIEAABgQEkCAAAwoCQBAAAYUJIAAAAMKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYEBJAgAAMGiVJWnTpk0aN26cJk6cqKNHj1odBwAAtEGtriSVlpZq69at2rBhg9LT0/XYY49ZHQkAALRBra4k7du3T9dcc418fX115ZVX6siRI1ZHAgAAbZDN6XQ6rQ5xLl555RV9++23mj59uiRp5MiReuWVVxqsk5ubq+XLlzfa9ssvv/zF86xatUqHDh06r22PHz+ubt26nde2vXv3dh2D1oDj1Hzne6wutOOE5uHPXvNZcaw4Ts1nxbFqdSUpPz9fBQUFmjVrliRp1KhRevnll91uFxMT45GSBAAA2qZWd7otNjZWH330kerr6/X5558rMjLS6kgAAKAN8rM6wLkKCQnRLbfcojvuuEN+fn565JFHrI4EAADaoFZ3uu18cboNAACci1Z3ug0AAMAbKEkAAAAGlCQAAAADShIAAIABJQkAAMCAkgQAAGBASQIAADCgJAEAABhQkgAAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYEBJAgAAMKAkAQAAGFCSAAAADPysDuBNMTExVkcAAAAtzJdffmlcbnM6nU4vZ8H/iomJafJfDP4Px6l5OE7Nx7FqHo5T83Gsmqe1HSdOtwEAABhQkgAAAAwoSQAAAAaUJAvNmDHD6gitAsepeThOzcexah6OU/NxrJqntR0nLtwGAAAwYCYJAADA4IK6T1JLcu+998rhcCgnJ8fqKC3W119/rdtuu019+vRRdXW1rrvuulY3VestX3zxhZYuXara2lrV1tZqwoQJuuWWW6yO1aL8+N9TTEyMbDabLrroIs2aNUtdu3a1OlqL8+9/9iQpNDSUn1VN+OmxuvnmmzV27FiLU7U8n332mZYtW6ba2lrV19erf//+uu+++6yO5RYlyQJVVVUqKSmR0+lUZWWlOnbsaHWkFuuqq67SypUr5XQ6dcstt2jy5Mmy2+1Wx2pRKioqlJ6erhUrVqhHjx6qq6vTrl27rI7VIv3435Mkvffee5o1a5bWr19vcaqW6d+PFc6OY3V25eXlmjt3rp588kl1795dkrRjxw6LUzUPp9ss8D//8z8aPny4RowYobfeesvqOK3CmTNnJEl+fvT6n3r33Xf129/+Vj169JD0wzEaOHCgxalavuTkZNXX1+u7776zOgrQpr333nsaMmSIqyBJajU/oyhJFnjzzTd1ww036KabbqIkufHxxx9r0qRJGjFihOLj4xUYGGh1pBbnu+++00UXXWR1jFYpIiJCRUVFVsdokX78szdp0iQtWLDA6jgt2r8fq/3791sdp8X5959RR48edf1MLysrsziZe/y13MtKSkr02Wefafbs2ZKkwsJCFRcXq0uXLhYna5l+nMZ2OBxKTU3V3r171a9fP6tjtSgXXXSRDh06ZHWMVunbb79VRESE1TFaJE4hNR/H6uz+/WdUZGSk1q9fr0mTJqm+vt7iZO4xk+Rl27ZtU2pqqtasWaM1a9Zo5syZ2rZtm9WxWjwfHx916tRJJSUlVkdpca677jq98847+vrrryVJdXV1reZ8v5Xefvtt+fr6UpIAD0tOTm7wM0r64edUa8BMkpe99tprWrp0qev1oEGDlJKSokmTJlmYquX6cRq7rq5O3bp1U2JiotWRWpxOnTpp0aJFmjdvnuvbbRMnTrQ6Vov08ccfa/LkyZJ++NvtsmXLLE7Ucv34Z0+S2rVrp7Vr11qcCK2V3W5XVlaWHnroIdXV1cnHx0dxcXHq1KmT1dHc4maSAAAABpxuAwAAMKAkAQAAGFCSAAAADChJAAAABpQkAAAAA0oSgBbnt7/9rWpqas5pm/Lycm3ZsuVnjZubm6u8vLyzrnOu2a699tqflQmAdShJANqE8ylJreWGdgCsQUkCYKnNmzdr5MiR+t3vfqfly5c3eO8f//iH7rvvPtfrSZMm6dChQyoqKtK4ceM0atQojRw5UocOHVJOTo4KCws1atQoPfXUUzpz5owyMzN122236ZZbbtH27dsl/TBbNG/ePE2ePFmPP/54k7lycnJ022236eabb9aiRYsavJebm6ubb75Z06ZNU2VlpSRp7969Gj9+vEaPHq2ZM2ee80wYgJaHO24DsMyBAwe0ceNG5eXlqWPHjiotLW3Wdq+++qp+85vfKDU1VWfOnFF9fb1SU1N15MgRPffcc5KkDRs26JJLLtG8efN06tQp3XHHHRo8eLAk6dChQ3r66afVrl27JseYPHmyUlNT5XA4NGPGDO3bt0+xsbGSpO7du+vVV19VTk6O/va3v2natGlaunSpVq5cKbvdruXLl2vLli3c+Rxo5ShJACzz0Ucf6aabblLHjh0lSSEhIc3a7oorrtDcuXPl6+urESNGqHfv3o3W2blzpw4ePKitW7dKkqqrq1VcXCxJGjJkyFkLkiTt2LFDa9asUU1NjUpKSnTw4EFXSRo+fLjrn4sXL9bhw4dVWFjoeozHmTNnlJSU1KzPAqDloiQBsIy7pyL5+vrK4XC4XtfW1kqS+vfvrw0bNuidd97RPffcoz//+c/q2bNno30/8sgjiouLa7TfDh06nHXcmpoaLVy4UM8//7zCw8O1cOFC19g/ZbPZ5HQ6deWVV/J8M6CN4ZokAJa55ppr9Prrr7uu6/np6baLL75YBw8eVH19vY4fP64vv/xSkvTNN98oPDxcEyZM0A033KAvv/xSgYGBqqqqcm07cOBA5eXluUrWgQMHmp2rpqZGPj4+CgkJUXl5ud5+++0G77/55puuf/br10+9evXS119/7cpXWVmpY8eOndvBANDiMJMEwDJ9+/bV2LFjNXbsWPn6+mr48OG65557XO93795dAwcO1MiRI3XZZZcpOjpa0g+n6Z566im1a9dOISEhWrZsmcLCwhQVFaWRI0dq1KhR+uMf/6hjx45p1KhRcjgcuuyyy7R06dJm5bLb7frd736nm2++Wd27d3edZvvRyZMndeuttyosLExPPPGE2rdvr6VLl+rPf/6zqqurZbPZdP/99zea3QLQutic7ua7AQAALkCcbgMAADCgJAEAABhQkgAAAAwoSQAAAAaUJAAAAANKEgAAgAElCQAAwICSBAAAYPD/AfebAZev+Rs7AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x504 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 7), constrained_layout=True)\n",
    "\n",
    "sns.boxplot(x=\"cluster\", y=\"imd\", hue='cluster', dodge=False, data=imd_england, showfliers=False, linewidth=.75, ax=ax)\n",
    "plt.legend([],[], frameon=False)\n",
    "sns.despine()\n",
    "\n",
    "ax.set_ylabel(\"index of multiple deprivation (rank)\")\n",
    "ax.set_xlabel(\"cluster label\")\n",
    "\n",
    "plt.savefig(f\"figs/imd_rank.pdf\", dpi=150)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}