match_waterbase_and_catchments.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "serious-hospital",
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona  # Needs to be imported first for geopandas to work.\n",
    "import geopandas as gpd\n",
    "import matplotlib as mpl\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pathlib\n",
    "import re\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "mpl.rcParams['figure.dpi'] = 144\n",
    "mpl.style.use('scrartcl.mplstyle')\n",
    "\n",
    "ROOT = pathlib.Path('data/wastewater_catchment_areas_public')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adverse-knife",
   "metadata": {},
   "outputs": [],
   "source": [
    "catchments = gpd.read_file(ROOT / 'catchments_consolidated.shp')\n",
    "catchments.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fifth-inspector",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the waterbase data.\n",
    "uwwtps = gpd.GeoDataFrame(pd.read_csv(ROOT / 'waterbase_consolidated.csv'))\n",
    "\n",
    "# Identify the treatment plants that were inactive in the most recent report.\n",
    "inactive = []\n",
    "for uwwCode, subset in uwwtps.groupby('uwwCode'):\n",
    "    subset = subset.sort_values('year')\n",
    "    item = subset.iloc[-1]\n",
    "    if item.uwwState == 'inactive':\n",
    "        inactive.append(item)\n",
    "inactive = pd.DataFrame(inactive)\n",
    "\n",
    "# Drop treatment works that are inactive.\n",
    "print(f'starting with {uwwtps.uwwCode.nunique()} treatment plants')\n",
    "uwwtps = uwwtps.loc[~np.in1d(uwwtps.uwwCode, inactive.uwwCode)]\n",
    "print(f'retained {uwwtps.uwwCode.nunique()} treatment plants after removing inactive plants')\n",
    "\n",
    "# Drop all treatment works from South West Water, Northern Ireland Water, and Gibraltar.\n",
    "# Remove all treatment works from South West Water and Northern Ireland Water.\n",
    "uwwtps = uwwtps[~(\n",
    "    uwwtps.uwwCode.str.startswith('UKENSW_SWS') | \n",
    "    uwwtps.uwwCode.str.startswith('UKNI') | \n",
    "    uwwtps.uwwCode.str.startswith('UKGIB')\n",
    ")]\n",
    "\n",
    "print(f'retained {uwwtps.uwwCode.nunique()} treatment plants after removing treatment plants owned '\n",
    "      'by water companies that did not supply data')\n",
    "\n",
    "# Ensure that all retained treatment plants have reported at least one record for capacity and load.\n",
    "grouped = uwwtps.groupby('uwwCode')\n",
    "np.testing.assert_array_less(1, grouped.uwwCapacity.max())\n",
    "np.testing.assert_array_less(1, grouped.uwwLoadEnteringUWWTP.max())\n",
    "\n",
    "# Reproject into the british national grid.\n",
    "uwwtps['geometry'] = gpd.points_from_xy(uwwtps.uwwLongitude, uwwtps.uwwLatitude)\n",
    "uwwtps = gpd.GeoDataFrame(uwwtps).set_crs('epsg:4326').to_crs('epsg:27700')\n",
    "\n",
    "# Group by treatment code and take the most recent data.\n",
    "uwwtps = uwwtps.sort_values('year').groupby('uwwCode').last().reset_index()\n",
    "uwwtps.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "invisible-chapel",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate the distance matrix between catchments and treatment works.\n",
    "distances = []\n",
    "for _, catchment in tqdm(catchments.iterrows(), total=len(catchments)):\n",
    "    distances.append(uwwtps.distance(catchment.geometry, align=False))\n",
    "# Cast to a dataframe.\n",
    "distances = np.asarray(distances)\n",
    "distances = pd.DataFrame(distances, index=catchments.identifier, columns=uwwtps.uwwCode)\n",
    "distances.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "later-wonder",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sweep over different distance thresholds and evaluate how many entities are close to each other\n",
    "# at that radius.\n",
    "thresholds = np.logspace(1, 3, 50)\n",
    "num_uwwtps_near_catchment = []\n",
    "num_catchments_near_uwwtp = []\n",
    "for threshold in thresholds:\n",
    "    near = distances < threshold\n",
    "    num_uwwtps_near_catchment.append(near.sum(axis=1))\n",
    "    num_catchments_near_uwwtp.append(near.sum(axis=0))\n",
    "    \n",
    "num_uwwtps_near_catchment = np.asarray(num_uwwtps_near_catchment)\n",
    "num_catchments_near_uwwtp = np.asarray(num_catchments_near_uwwtp)\n",
    "assert num_uwwtps_near_catchment.shape == (len(thresholds), len(catchments))\n",
    "assert num_catchments_near_uwwtp.shape == (len(thresholds), len(uwwtps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "centered-accuracy",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the fraction of UWWTPs that have exactly one catchment within a distance given by the \n",
    "# threshold. As expected, there is an optimal distance: small thresholds are too restrictive whereas\n",
    "# large ones \"confuse\" different catchments.\n",
    "fig, ax = plt.subplots()\n",
    "fraction_single_catchment_near_uwwtp = np.mean(num_catchments_near_uwwtp == 1, axis=1)\n",
    "line, = ax.plot(thresholds, fraction_single_catchment_near_uwwtp)\n",
    "i = np.argmax(fraction_single_catchment_near_uwwtp)\n",
    "ax.scatter(thresholds[i], fraction_single_catchment_near_uwwtp[i], color=line.get_color())\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlabel('Thresholds (metres)')\n",
    "ax.set_ylabel('Fraction of UWWTPs with unique\\ncatchment @ threshold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "painted-whale",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find all combinations of catchments and treatment plants where both the treatment plant \n",
    "# and the catchment are associated with exactly one other treatment plant/catchment at a given \n",
    "# distance.\n",
    "thresholds = np.linspace(0, 500, 50)\n",
    "num = []\n",
    "for threshold in tqdm(thresholds):\n",
    "    binary = (distances <= threshold).values\n",
    "    f = (binary.sum(axis=1) == 1)[:, None] & ((binary).sum(axis=0) == 1) & binary\n",
    "    catchment_idx, stp_idx = np.nonzero(f)\n",
    "    num.append(len(catchment_idx))\n",
    "    \n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(thresholds, num)\n",
    "ax.set_xlabel('Thresholds (metres)')\n",
    "ax.set_ylabel('# unique pairs @ threshold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "formed-capital",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: distance-based matching.\n",
    "\n",
    "# Find all combinations of catchments and treatment plants where both the treatment plant and the \n",
    "# catchment are associated with exactly one other treatment plant/catchment below a given distance.\n",
    "threshold = 100\n",
    "binary = (distances <= threshold).values\n",
    "f_distance = (binary.sum(axis=1) == 1)[:, None] & ((binary).sum(axis=0) == 1) & binary\n",
    "print('distance matched', f_distance.sum())\n",
    "\n",
    "# Validate uniqueness to ensure the logic is consistent.\n",
    "np.testing.assert_array_less(f.sum(axis=0), 2)\n",
    "np.testing.assert_array_less(f.sum(axis=1), 2)\n",
    "\n",
    "def normalise_name(value):\n",
    "    \"\"\"\n",
    "    Function to normalise treatment names by removing special characters and common abbreviations.\n",
    "    \"\"\"\n",
    "    if pd.isnull(value):\n",
    "        return value\n",
    "    # Drop all special characters.\n",
    "    value, _ = re.subn(r'[\\(\\)\\[\\]\\{\\}\\.\\&-]', '', value.lower())\n",
    "    # Drop stw and wwtw (often suffixes).\n",
    "    value, _ = re.subn(r'\\b(stw|wwtw|doa|wrw|bucks)\\b', '', value)\n",
    "    # Drop all whitespace.\n",
    "    value, _ = re.subn(r'\\s', '', value)\n",
    "    return value\n",
    "\n",
    "# Step 2: name-based matching.\n",
    "\n",
    "# Add all treatment works where the normalised names are identical and the distance is less than \n",
    "# another threshold.\n",
    "f_name = (catchments['name'].apply(normalise_name).values[:, None] == uwwtps.uwwName.apply(normalise_name).values) \\\n",
    "    & (distances.values <= 2500)\n",
    "\n",
    "print('name matched', f_name.sum())\n",
    "\n",
    "# Create the auto-matched lookup table.\n",
    "catchment_idx, stp_idx = np.nonzero(f_distance | f_name)\n",
    "auto_matched = pd.concat([\n",
    "    uwwtps.iloc[stp_idx].reset_index(drop=True), \n",
    "    catchments.iloc[catchment_idx].reset_index(drop=True).drop('geometry', axis=1)\n",
    "], axis=1)[['identifier', 'name', 'uwwCode', 'uwwName']]\n",
    "print('total automatically matched', auto_matched.uwwCode.nunique())\n",
    "auto_matched.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "provincial-stomach",
   "metadata": {},
   "outputs": [],
   "source": [
    "manual_matched = pd.DataFrame([\n",
    "    ('214873b5cf', 'UU-09-SC47-DAVYH', 'UKENNW_UU_TP000042', 'Manchester and Salford (Davyhulme)STW'),\n",
    "    ('d793a815d9', 'UU-10-SC50-ROYTO', 'UKENNW_UU_TP000105', 'ROYTON STW'),\n",
    "    ('b65d2eebc7', 'UU-10-SC51-FAILS', 'UKENNW_UU_TP000048', 'FAILSWORTH STW'),\n",
    "    ('bb6ca16304', 'UU-10-SC54-HAZEL', 'UKENNW_UU_TP000058', 'HAZEL GROVE STW'),\n",
    "    ('59bd836fc3', 'UU-10-SC54-STOCK', 'UKENNW_UU_TP000116', 'STOCKPORT STW'),\n",
    "    ('5900c96825', 'Shenfield and Hutton', 'UKENAN_AW_TP000230', 'SHENFIELD STW'),\n",
    "    ('4647e75f9e', 'CARDIFF BAY', 'UKWAWA_WW_TP000026', 'CARDIFF STW'),\n",
    "    # These two are the same treatment work.\n",
    "    ('cd33d70a80', 'Catterick Village', 'UKENNE_YW_TP000007', 'CATTERICK STW'),\n",
    "    ('cd33d70a80', 'Catterick Village', 'UKENNE_YW_TP000008', 'COLBURN STW'),\n",
    "    ('27ec2163de', 'UU-11-SC56-NTHWI', 'UKENNW_UU_TP000097', 'NORTHWICH STW'),\n",
    "    # Cuddington has been decommissioned and now feeds into Northwich, but it's all a bit tricky\n",
    "    # because there can be pumping going on between the different treatment plants\n",
    "    # (https://www.unitedutilities.com/about-us/cheshire/).\n",
    "    (None, None, 'UKENNW_UU_TP000040', 'CUDDINGTON STW'),\n",
    "    # Cotgrave feeds into Radcliffe.\n",
    "    # https://www.nottinghamshire.gov.uk/planningsearch/DisplayImage.aspx?doc=cmVjb3JkX251bWJlcj02OTU4JmZpbGVuYW1lPVxcbnMwMS0wMDI5XGZpbGVkYXRhMiRcREIwMy0wMDMwXFNoYXJlZEFwcHNcRExHU1xQbGFuc1xQTEFOTklOR1xTY3ItMzY0OFxSYWRjbGlmZmUgRUlBIFNjcmVlbmluZyBPcGluaW9uIFJlcXVlc3QucGRmJmltYWdlX251bWJlcj01JmltYWdlX3R5cGU9cGxhbm5pbmcmbGFzdF9tb2RpZmllZF9mcm9tX2Rpc2s9MTIvMDQvMjAxNyAxMDozMzo0MQ==\n",
    "    (None, None, 'UKENMI_ST_TP000064', 'COTGRAVE STW'),\n",
    "    ('fb13da240a', 'MANSFIELD-BATH LANE (WRW)', 'UKENMI_ST_TP000143', 'MANSFIELD STW'),\n",
    "    # https://goo.gl/maps/pFhJy3ZAxyABKwWc6\n",
    "    ('9fcbac5ecc', 'STOKE BARDOLPH (WRW)', 'UKENMI_ST_TP000163', 'NOTTINGHAM STW'),\n",
    "    ('7be11f772b', 'BEESTON -LILAC GROVE (WRW)', 'UKENMI_ST_TP000025', 'BEESTON STW'),\n",
    "    ('c861b125a1', 'STRATFORD-MILCOTE (WRW)', 'UKENMI_ST_TP000206', 'STRATFORD  STW'),\n",
    "    # Long Marston has been decommissioned.\n",
    "    # https://waterprojectsonline.com/custom_case_study/pebworth-long-marston-stratford-upon-avon-transfer/\n",
    "    (None, None, 'UKENMI_ST_TP000135', 'LONG MARSTON STW'),\n",
    "    ('01bf76d616', 'MAPLE LODGE STW', 'UKENTH_TWU_TP000106', 'MAPLE LODGE, BUCKS STW\"'),\n",
    "    ('034de6cd3b', 'Eastwood', 'UKENNE_YW_TP000079', 'TODMORDEN   STW'),\n",
    "    ('c64f685e8b', 'Neiley', 'UKENNE_YW_TP000109', 'HOLMFIRTH   STW'),\n",
    "    ('7e5b0d4c85', 'UU-05-SC28-BLACK', 'UKENNW_UU_TP000018', 'BLACKBURN STW'),\n",
    "    ('ccad6a0777', 'UU-05-SC28-DARWE', 'UKENNW_UU_TP000041', 'DARWEN STW'),\n",
    "    ('31ec9a8fc2', 'UU-06-SC35-HUYTO', 'UKENNW_UU_TP000066', 'HUYTON STW'),\n",
    "    ('a83266015d', 'UU-06-SC35-LIVER', 'UKENNW_UU_TP000080', 'LIVERPOOL SOUTH [WOOLTON]) STW'),\n",
    "    ('0d2d457ce0', 'MARGATE AND BROADSTAIRS', 'UKENSO_SW_TP000019', 'BROADSTAIRS/MARGATE OUT CSM (Weatherlees B)'),\n",
    "    ('bbd09a27de', 'WEATHERLEES HILL', 'UKENSO_SW_TP000022', 'RAMSGATE, SANDWICH, DEAL   STW\"'),\n",
    "    ('00e1c4c53b', 'UU-09-SC46-ECCLE', 'UKENNW_UU_TP000046', 'ECCLES STW'),\n",
    "    ('feab30c503', 'UU-07-SC40-WORSL', 'UKENNW_UU_TP000140', 'WORSLEY STW'),\n",
    "    ('a3aea16d4a', 'ABINGDON STW', 'UKENTH_TWU_TP000001', 'ABINGDON (OXON STW)'),\n",
    "    ('b428837dda', 'NAGS HEAD LANE STW', 'UKENTH_TWU_TP000115', 'NAGS HEAD LANE ( BRENTWOOD STW'),\n",
    "    ('3d27aa0af3', 'RIVERSIDE STW', 'UKENTH_TWU_TP000125', 'LONDON (Riverside STW)'),\n",
    "    ('cee1fb4e33', 'UU-07-SC40-GLAZE', 'UKENNW_UU_TP000053', 'GLAZEBURY STW'),\n",
    "    ('feec09aaaf', 'UU-07-SC40-LEIGH', 'UKENNW_UU_TP000078', 'LEIGH STW'),\n",
    "    ('3b82554080', 'UU-07-SC39-BOLTO', 'UKENNW_UU_TP000019', 'BOLTON STW'),\n",
    "    ('51bdbc0969', 'UU-07-SC37-BURYZ', 'UKENNW_UU_TP000026', 'BURY STW'),\n",
    "    ('f08942bc5a', 'UU-10-SC52-ASHUL', 'UKENNW_UU_TP000007', 'ASHTON-UNDER-LYNE STW'),\n",
    "    ('5d269b4c58', 'UU-10-SC52-DUKIN', 'UKENNW_UU_TP000044', 'DUKINFIELD STW'),\n",
    "    ('b24ca82746', 'RHIWSAESON (NEW)', 'UKWAWA_WW_TP000024', 'RHIWSAESON STW RHIWSAESON LLANTRI STW'),\n",
    "    ('a7ecedaef5', 'COSLECH', 'UKWAWA_WW_TP000020', 'SOUTH ELY VALLEY   STW'),\n",
    "    ('2e2e1b5104', 'Dalscone DOA', 'UKSC_TP00053', 'DALSCONE S.T.W. (NEW)'),\n",
    "    ('2e2e1b5104', 'Dalscone DOA', 'UKSC_TP00054', 'DALSCONE S.T.W.  (OLD)'),\n",
    "    ('5f66794c55', 'Whilton', 'UKENAN_AW_TP000286', 'DAVENTRY   STW'),\n",
    "    ('3b3e3efde8', 'BAKEWELL, PICKORY CORNER (WRW)', 'UKENMI_ST_TP000013', 'BAKEWELL STW'),\n",
    "    ('2eb71ffcfc', 'BOTTESFORD-STW', 'UKENMI_ST_TP000033', 'BOTTESFORD STW'),\n",
    "    ('e4a961fba1', 'KEMPSEY WORKS (WRW)', 'UKENMI_ST_TP000119', 'KEMPSEY STW'),\n",
    "    ('53fe082858', 'CLAYMILLS (WRW)', 'UKENMI_ST_TP000056', 'BURTON ON TRENT   STW'),\n",
    "    ('54e25cb6a9', 'MELTON (WRW)', 'UKENMI_ST_TP000152', 'MELTON MOWBRAY STW'),\n",
    "    ('b0753764b4', 'STANLEY DOWNTON (WRW)', 'UKENMI_ST_TP000208', 'STROUD STW'),\n",
    "    ('b11e60b850', 'ROUNDHILL (WRW)', 'UKENMI_ST_TP000180', 'STOURBRIDGE & HALESOWEN   STW'),\n",
    "    ('ffe8483a02', 'ASH VALE STW', 'UKENTH_TWU_TP000009', 'ASH VALE, STRATFORD ROAD, NORTH STW\"'),\n",
    "    ('e526b8f3b6', 'BRACKNELL STW', 'UKENTH_TWU_TP000025', 'BRACKNELL, HAZELWOOD LANE, BINF STW\"'),\n",
    "    ('29d52395a3', 'UU-04-SC21-PREST', 'UKENNW_UU_TP000102', 'PRESTON (CLIFTON MARSH) STW'),\n",
    "    ('c1e0195f63', 'WINDSOR STW', 'UKENTH_TWU_TP000152', 'WINDSOR, HAM ISLAND, OLD WINDSO STW\"'),\n",
    "    ('743841f528', 'CAMBERLEY STW', 'UKENTH_TWU_TP000033', 'CAMBERLEY, CAMBERLEY, SURREY STW\"'),\n",
    "    ('80c85e78c6', 'DORCHESTER STW', 'UKENTH_TWU_TP000056', 'DORCHESTER STW (OXON)'),\n",
    "    ('0f5561c4b9', 'STANFORD RIVERS STW', 'UKENTH_TWU_TP000136', 'STANFORD RIVERS, ONGAR, ESSEX STW\"'),\n",
    "    ('45b2ebfd2f', 'UU-03-SC17-KRKLO', 'UKENNW_UU_TP000151', 'KIRKBY LONSDALE STW HUMUS TANK EFFLUENT'),\n",
    "    ('935c2734be', 'MARLBOROUGH STW', 'UKENTH_TWU_TP000108', 'MARLBOROUGH, MARLBOROUGH, WILTS STW\"'),\n",
    "    ('ba3e9907f2', 'EARLSWOOD STW', 'UKENTH_TWU_TP000123', 'REIGATE STW'),\n",
    "    ('98a14c9923', 'STANTON - DERBYSHIRE (WRW)', 'UKENMI_ST_TP000201', 'STANTON   STW'),\n",
    "    ('1d9dcc765f', 'WORMINGHALL STW', 'UKENTH_TWU_TP000158', 'WORMINGHALL, WORMINGHALL, BUCKS STW\"'),\n",
    "    ('a087889829', 'SUTTON BONNINGTON (WRW)', 'UKENMI_ST_TP000278', 'SUTTON BONINGTON STW, FE\"'),\n",
    "    ('fd7e5039c1', 'STANSTED MOUNTFITCHET STW', 'UKENTH_TWU_TP000137', 'STANSTED MOUNTFITCHET, STANSTED STW\"'),\n",
    "    ('0d09ae579a', 'RAMSBURY STW', 'UKENTH_TWU_TP000121', 'RAMSBURY, RAMSBURY, MARLBOROUGH STW\"'),\n",
    "    ('0de7c65bea', 'UU-06-SC31-HSKBK', 'UKENNW_UU_TP000060', 'HESKETH BANK STW'),\n",
    "    ('a8d4ccb494', 'CHICKENHALL EASTLEIGH', 'UKENSO_SW_TP000013', 'EASTLEIGH   STW'),\n",
    "    ('0340b3c3c4', 'WALSALL WOOD (WRW)', 'UKENMI_ST_TP000221', 'WALSALL NORTH   STW'),\n",
    "    ('47aebf7f1e', 'UU-03-SC18-LANCA', 'UKENNW_UU_TP000076', 'LANCASTER (STODDAY) STW'),\n",
    "    ('95fa57e680', 'MILL GREEN STW', 'UKENTH_TWU_TP000111', 'MILL GREEN, HATFIELD, HERTS STW\"'),\n",
    "    ('f01161329e', 'SLOUGH STW', 'UKENTH_TWU_TP000133', 'SLOUGH, WOOD STW\"'),\n",
    "    ('d2e0ce618b', 'FINSTOCK STW', 'UKENTH_TWU_TP000067', 'FINSTOCK, FINSTOCK, OXON STW\"'),\n",
    "    ('15831f78a8', 'WOTTON UNDER EDGE STW CATCHMENT', 'UKENSW_WXW_TP000109', 'WOTTON UNDER EDGE STW'),\n",
    "    ('8237f30715', 'STANFORD IN THE VALE STW', 'UKENTH_TWU_TP000169', 'Standford in the Vale STW'),\n",
    "    ('a1e551a2fd', 'MORESTEAD ROAD WINCHESTER', 'UKENSO_SW_TP000003', 'WINCHESTER CENTRAL AND SOUTH (MORESTEAD) STW'),\n",
    "    ('3b917f6c5f', 'BASINGSTOKE STW', 'UKENTH_TWU_TP000013', 'BASINGSTOKE, WILDMOOR, BASINGST STW\"'),\n",
    "    ('1903e1316b', 'FULLERTON', 'UKENSO_SW_TP000006', 'ANDOVER STW'),\n",
    "    ('e9216a235f', 'Fochabers DOA', 'UKSC_TP00077', 'FOCHABERS WWTP'),\n",
    "    ('965a041c6a', 'THAME STW', 'UKENTH_TWU_TP000140', 'THAME, THAME, OXON STW\"'),\n",
    "    ('c836b06a6e', 'GODALMING STW', 'UKENTH_TWU_TP000071', 'GODALMING, UNSTEAD, GODALMING, STW\"'),\n",
    "    ('bb40f7f215', 'UU-03-SC16-GRNGS', 'UKENNW_UU_TP000055', 'GRANGE-OVER-SANDS STW'),\n",
    "    ('60cd9fef4c', 'DRENEWYDD - OSWESTRY (WRW)', 'UKENMI_ST_TP000166', 'OSWESTRY DRENEWYDD STW'),\n",
    "    ('9edf9e05e3', 'HIGHER HEATH-PREES (WRW)', 'UKENMI_ST_TP000268', 'PREES HIGHER HEATH STW'),\n",
    "    ('5b6ed3d068', 'RUSHMOOR (WRW)', 'UKENMI_ST_TP000184', 'TELFORD   STW'),\n",
    "    ('d8b6f988ea', 'CHESHAM STW', 'UKENTH_TWU_TP000039', 'CHESHAM, BUCKS STW\"'),\n",
    "    ('04a2d621cd', 'UU-07-SC40-TYLDE', 'UKENNW_UU_TP000120', 'TYLDESLEY STW'),\n",
    "    ('5353ff23aa', 'APPLETON STW', 'UKENTH_TWU_TP000006', 'APPLETON, ABINGDON, OXON STW\"'),\n",
    "    ('8fccf79af1', 'WHITE WALTHAM STW', 'UKENTH_TWU_TP000150', 'WHITE WALTHAM, WHITE WALTHAM, B STW\"'),\n",
    "    ('a2a6cca4f6', 'POWICK (WRW)', 'UKENMI_ST_TP000173', 'POWICK NEW STW'),\n",
    "    ('c33bed15a5', 'Knostrop Merge High + Low', 'UKENNE_YW_TP000098', 'LEEDS (KNOSTROP) STW'),\n",
    "    ('11d6a4c3d9', 'HAYDEN (WRW)', 'UKENMI_ST_TP000256', 'CHELTENHAM STW'),\n",
    "    ('5b3b7769c8', 'Withernsea No. 2', 'UKENNE_YW_TP000139', 'WITHERNSEA OUTFALL STW'),\n",
    "    ('864e5eaa8d', 'TRIMDON VILLAGE STW NZ37346301', 'UKENNE_NU_TP000052', 'TRIMDON VILLAGE STW'),\n",
    "    ('d31b19b35b', 'NORTH TIDWORTH STW CATCHMENT', 'UKENSW_VE_TP000001', 'TIDWORTH GARRISON STW FE'),\n",
    "    ('f00ad841f9', 'PEACEHAVEN', 'UKENSO_SW_TP000126', 'PEACEHAVEN WASTEWATER TREATMENT WKS'),\n",
    "    # Morecambe is auto-matched to Middleton. But we also need the Morecambe catchment area itself.\n",
    "    # (cf. page 3 of https://www.unitedutilities.com/globalassets/documents/pdf/7013b-morcambe-a5-6pp-flyer-v10_acc17.pdf).\n",
    "    ('15d6d1fe75', 'UU-03-SC18-MOREC', 'UKENNW_UU_TP000092', 'MORECAMBE STW'),\n",
    "    # No good matches.\n",
    "    (None, None, 'UKSC_TP00011', 'ANNANDALE WATER MSA S.T.W.'),\n",
    "    (None, None, 'UKENTH_TWU_TP000162', 'ALDERSHOT MILITARY STW'),\n",
    "    (None, None, 'UKSC_TP00201', 'FASLANE STW'),\n",
    "    (None, None, 'UKSC_TP00085', 'GLENEAGLES STW'),\n",
    "    (None, None, 'UKWAWA_WW_TP000128', 'BLUESTONE LEISURE CANASTON'),\n",
    "    (None, None, 'UKENMI_ST_TP000227', 'WELSHPOOL STW'),\n",
    "    (None, None, 'UKSC_TP00171', 'SOUTHERNESS STW'),\n",
    "    (None, None, 'UKENMI_ST_TP000127', 'KNIGHTON STW'),\n",
    "    (None, None, 'UKENMI_ST_TP000133', 'LLANIDLOES STW'),\n",
    "    (None, None, 'UKENMI_ST_TP000161', 'NEWTOWN STW'),\n",
    "    (None, None, 'UKWAWA_WW_TP000011', 'MERTHYR TYDFIL STW'),\n",
    "], columns=['identifier', 'name', 'uwwCode', 'uwwName'])\n",
    "\n",
    "matched = pd.concat([auto_matched, manual_matched])\n",
    "# Ensure that the matching accounts for exactly the right set of treatment works.\n",
    "assert set(matched.uwwCode) == set(uwwtps.uwwCode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67b0e472",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure there are no invalid identifiers amongst the manual matches. That might happen\n",
    "# if the underlying data change, e.g. as part of the request via whatdotheyknow.com.\n",
    "invalid_identifiers = set(manual_matched.identifier) - set(catchments.identifier) - {None}\n",
    "if invalid_identifiers:\n",
    "    raise ValueError(f'found {len(invalid_identifiers)} invalid identifiers: {invalid_identifiers}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "optimum-finnish",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate the distance between the catchment and the stp to give \"a feel\" for how good the matching\n",
    "# is.\n",
    "catchments_indexed = catchments.set_index('identifier')\n",
    "uwwtps_indexed = uwwtps.set_index('uwwCode')\n",
    "matched['distance'] = matched.apply(\n",
    "    lambda x: None if pd.isnull(x.identifier) else \n",
    "    catchments_indexed.geometry.loc[x.identifier].distance(uwwtps_indexed.geometry.loc[x.uwwCode]), axis=1)\n",
    "matched.sort_values('uwwCode').to_csv(ROOT / 'waterbase_catchment_lookup.csv', index=False)\n",
    "# Show the 99% distance between matched catchments and stps.\n",
    "matched.distance.quantile([.5, .99])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adapted-jersey",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show catchments that are matched to more than one treatment work. This should NEVER happen, unless\n",
    "# the same treatment work has been assigned different identifiers by UWWTD.\n",
    "\n",
    "exceptions = [\n",
    "    '2e2e1b5104',  # Dalscone New and Old.\n",
    "    'a3aea16d4a',  # Abingdon appears twice.\n",
    "    '8a5a88fe31',  # Balby appears twice.\n",
    "    'd8270e9ad8',  # Catterick village and Colburn STWs are in the same location...\n",
    "    'cd33d70a80',  # ... and the catchment needs to be linked to both treatment works.\n",
    "]\n",
    "\n",
    "counts = 0\n",
    "for identifier, subset in matched.groupby('identifier'):\n",
    "    if identifier in exceptions:\n",
    "        continue\n",
    "    if len(subset) > 1:\n",
    "        print(subset)\n",
    "        counts += 1\n",
    "        \n",
    "assert counts == 0, 'data doesn\\'t make sense'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "parental-yemen",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show treatment works broken down by company.\n",
    "counts = pd.merge(matched, catchments, on='identifier', how='left').groupby('company').uwwCode.nunique()\n",
    "print('total number of treatment works', counts.sum())\n",
    "counts"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}