consolidate_catchments.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "instructional-chapter",
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\n",
    "import os\n",
    "import logging\n",
    "from tqdm.notebook import tqdm\n",
    "import shapely as sh\n",
    "import shapely.geometry\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "import collections\n",
    "import hashlib\n",
    "\n",
    "mpl.rcParams['figure.dpi'] = 144\n",
    "mpl.style.use('scrartcl.mplstyle')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "utility-melbourne",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define helper functions for processing the data.\n",
    "\n",
    "def assert_in(*values, key=None):\n",
    "    \"\"\"\n",
    "    Assert that the value belongs to a set of values.\n",
    "    \"\"\"\n",
    "    def _wrapper(properties, value):\n",
    "        assert value in values, f'{value} is not one of {values}'\n",
    "        if key:\n",
    "            properties[key] = value\n",
    "    return _wrapper\n",
    "\n",
    "\n",
    "def assert_unique(key=None, *, on_error='raise', exceptions=None):\n",
    "    \"\"\"\n",
    "    Assert that the value is unique.\n",
    "    \"\"\"\n",
    "    values = collections.Counter()\n",
    "    exceptions = exceptions or set()\n",
    "    def _wrapper(properties, value):\n",
    "        if value in values and value not in exceptions:\n",
    "            message = f'`{value}` is not a unique value and has occurred {values[value]} times'\n",
    "            if on_error == 'raise':\n",
    "                raise ValueError(message)\n",
    "            elif on_error == 'warn':\n",
    "                logging.warning(message)\n",
    "            else:\n",
    "                raise NotImplementedError\n",
    "            \n",
    "        values[value] += 1\n",
    "        if key:\n",
    "            properties[key] = value\n",
    "    return _wrapper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "electrical-oregon",
   "metadata": {},
   "outputs": [],
   "source": [
    "root = 'zip://data/raw_catchments'\n",
    "crs = 'epsg:27700'\n",
    "\n",
    "# List of recognised catchment owners as supplied by Severn Trent (to make sure there aren't any \n",
    "# unexpected values).\n",
    "severn_trent_allowed_owners = [\n",
    "    'STW', 'Hafren Dyfrdwy', 'Yorkshire Water', 'Anglian Water', 'Private', 'Welsh Water', \n",
    "    'United Utilities', 'Thames Water', None,\n",
    "]\n",
    "\n",
    "def united_utilities_aggregation(catchments):\n",
    "    \"\"\"\n",
    "    Function to aggregate United Utilities subcatchments into catchments at the treatment work \n",
    "    level.\n",
    "    \"\"\"\n",
    "    catchments['reference'] = catchments.reference.str.extract(r'(?P<prefix>UU-\\d+-SC\\d+-\\w+)')\n",
    "    assert catchments.reference.nunique() < len(catchments)\n",
    "    return catchments.dissolve(by='reference').reset_index()\n",
    "\n",
    "\n",
    "# Declare the schema we expect to see for all the different shapefiles. `None` means we acknowledge\n",
    "# the field exists but ignore it in processing. We use the `postprocessor` function to apply \n",
    "# transformations after loading the data, e.g. to aggregate subcatchments.\n",
    "datasets = [\n",
    "    {\n",
    "        'company': 'scottish_water',\n",
    "        'filename': r'scottish_water.zip/DOAs and WWTWs',\n",
    "        'properties': {\n",
    "            'doa_name': assert_unique('name', exceptions=[\n",
    "                # This occurs twice, once for DOA000038 (which should probably be called Maxton)\n",
    "                # and also for DOA002612 which actually has a road called 'Wellrig'.\n",
    "                'Wellrig DOA'\n",
    "            ]),\n",
    "            # See postprocessor for exception.\n",
    "            'doa_refere': assert_unique('reference', exceptions=['DOA000917']),\n",
    "            'SHAPE_STAr': None,\n",
    "            'wams_stw_n': None,\n",
    "        },\n",
    "        # This deals with the one repeated reference in the dataset due to a new development near Dunbar.\n",
    "        'postprocessor': lambda x: x.dissolve(by='reference').reset_index(),\n",
    "    },\n",
    "    {\n",
    "        'company': 'thames_water',\n",
    "        'filename': 'thames_water.zip',\n",
    "        'properties': {\n",
    "            # We do not assert uniqueness for Thames Water Catchments because they have\n",
    "            # provided us with detailed subcatchments (which belong to different STWs).\n",
    "            'STWCODE': 'reference',\n",
    "            'STWNAME': 'name',\n",
    "            'SDACFOUL_1': None,\n",
    "            'RECIEVINGS': None,\n",
    "            'OBJECTID': None, \n",
    "            'CREATIONUS': None,\n",
    "            'DATECREATE': None,\n",
    "            'DATEMODIFI': None,\n",
    "            'LASTUSER': None,\n",
    "            'SHAPEAREA': None,\n",
    "            'SDACFOULBN': None,\n",
    "            'DATEPOSTED': None,\n",
    "            'GLOBALID': None,\n",
    "            'SHAPE_AREA': None, \n",
    "            'SHAPE_LEN': None,\n",
    "        },\n",
    "        'postprocessor': lambda x: x.dissolve(by='reference').reset_index(),\n",
    "    },\n",
    "    {\n",
    "        'company': 'anglian_water',\n",
    "        'filename': 'anglian_water.zip',\n",
    "        'properties': {\n",
    "            'AREANAME': assert_unique('name'),\n",
    "            'AREASHORTC': assert_unique('reference'),\n",
    "            'LIQUIDTYPE': None,\n",
    "            'COLLECTION': None,\n",
    "            'COLLECTIO1': None, \n",
    "            'TREATMENTM': None, \n",
    "            'TREATMENT1': None, \n",
    "            'NOSINLINES': None, \n",
    "            'NOSTERMINA': None, \n",
    "            'NOSPROPS': None, \n",
    "            'NOSDOMESTI': None, \n",
    "            'NOSMIXEDPR': None, \n",
    "            'NOSCOMMPRO': None, \n",
    "            'HOMEPOPN': None, \n",
    "            'HOUSEHOLDP': None, \n",
    "            'DOMESTICPE': None, \n",
    "            'HOLIDAYPE': None, \n",
    "            'INSTITUTIO': None, \n",
    "            'TRADEPE': None, \n",
    "            'IMPORTPE': None, \n",
    "            'EXPORTPE': None, \n",
    "            'NOSSURFACE': None, \n",
    "            'STWS_TOTNU': None, \n",
    "            'SEWOUTFALL': None, \n",
    "            'DATA_CONFI': None,\n",
    "            'ID': None,\n",
    "        }\n",
    "    },\n",
    "    {\n",
    "        'company': 'northumbrian_water',\n",
    "        'filename': 'northumbrian_water.zip/Export',\n",
    "        'properties': {\n",
    "            'NAME': assert_unique('name'),\n",
    "            'ID': assert_unique('reference'),\n",
    "            'GID': None,\n",
    "            'APIC_STYLE': None, \n",
    "            'APIC_SPACE': None, \n",
    "            'APIC_STATE': None, \n",
    "            'APIC_CDATE': None, \n",
    "            'APIC_MDATE': None,\n",
    "        }\n",
    "    },\n",
    "    {\n",
    "        'company': 'severn_trent_water',\n",
    "        'filename': 'severn_trent_water.zip',\n",
    "        'properties': {\n",
    "            # We allow non-unique names for private treatment works because they're not easily \n",
    "            # identifiable.\n",
    "            'CATCHMENT1': assert_unique('name', exceptions={\n",
    "                'private SPS', 'private water works', 'private', 'private STW',\n",
    "                'Surface water only', None,\n",
    "            }), \n",
    "            # Owner was useful in the original dataset because it allowed us to identify records\n",
    "            # that weren't actually from Severn Trent Water. We'll have to do that manually now.\n",
    "            # 'Owner': assert_in(*severn_trent_allowed_owners, key='owner'), \n",
    "            'SU_REFEREN': assert_in('WwTw', 'WwTW', 'ST', None, 'Cess Pit'), \n",
    "            'SAP_FLOC_I': 'reference',\n",
    "            'OBJECTID': None,\n",
    "            'CATCHMENT_': None,\n",
    "            'CACI_2001': None,\n",
    "            'CACI_20010': None,\n",
    "            'CACI_2006': None,\n",
    "            'CACI_2007': None,\n",
    "            'CACI_2030': None,\n",
    "            'PROPERTIES': None,\n",
    "            'STATUS': None,\n",
    "            'UPDATED': None,\n",
    "            'SUGRAN': None,\n",
    "            'COMMENT': None,\n",
    "            'SAP_FLOC_D': None,\n",
    "            'SHAPE_Leng': None,\n",
    "            'STAR_ID': None,\n",
    "            'WORKS_CODE': None,\n",
    "            'WORKS_ID': None,\n",
    "            'CREATED_DA': None,\n",
    "            'MODIFIED_D': None,\n",
    "            'MODIFIED_U': None,\n",
    "            'CREATED_US': None,\n",
    "            'JOB_ID': None,\n",
    "            'QA_STATUS': None,\n",
    "            'SAP_USER_S': None,\n",
    "            'SAP_CLASS': None,\n",
    "            'Shape_STAr': None,\n",
    "            'Shape_STLe': None,\n",
    "            'Shape_ST_1': None,\n",
    "            'Shape_ST_2': None,\n",
    "        },\n",
    "    },\n",
    "    {\n",
    "        'company': 'southern_water',\n",
    "        'filename': 'southern_water.zip',\n",
    "        'properties': {\n",
    "            'Site_Unit_': assert_unique('name'),\n",
    "        }\n",
    "    },\n",
    "    {\n",
    "        'company': 'united_utilities',\n",
    "        'filename': 'united_utilities.zip',\n",
    "        'properties': {\n",
    "            'DA_CODE': assert_unique('reference'),\n",
    "        },\n",
    "        'postprocessor': united_utilities_aggregation,\n",
    "    },\n",
    "    {\n",
    "        'company': 'welsh_water',\n",
    "        'filename': 'welsh_water.zip',\n",
    "        'properties': {\n",
    "            'CATCHMENT_': assert_unique('name'),\n",
    "            'TERMINAL_A': None,\n",
    "            'Area_M': None,\n",
    "        }\n",
    "    },\n",
    "    {\n",
    "        'company': 'wessex_water',\n",
    "        'filename': 'wessex_water.zip',\n",
    "        'properties': {\n",
    "            # There is a shared code because the Corsley Heath catchment comprises\n",
    "            # two disjoint areas.\n",
    "            'SITEID': assert_unique('reference', exceptions=[23390]), \n",
    "            'NAME': 'name',\n",
    "            'Comment': None,\n",
    "        },\n",
    "        # Wessex Water provide subcatchments, and we need to aggregate.\n",
    "        'postprocessor': lambda x: x.dissolve(by='name').reset_index(),\n",
    "    },\n",
    "    {\n",
    "        'company': 'yorkshire_water',\n",
    "        'filename': 'yorkshire_water.zip/EIR - Wastewater Catchments',\n",
    "        'properties': \n",
    "        {\n",
    "            'Company': assert_in('Yorkshire'), \n",
    "            'Name': assert_unique('name'),\n",
    "        }\n",
    "    }\n",
    "]\n",
    "assert len(datasets) == 10\n",
    "\n",
    "# Iterate over all datasets, then over all features in each dataset, and validate the records.\n",
    "parts = []\n",
    "total = 0\n",
    "for dataset in datasets:\n",
    "    try:\n",
    "        # Load the data. Probably should've used geopandas rather than the lower-level fiona. But \n",
    "        # this implementation works for the time being.\n",
    "        company_catchments = []\n",
    "        with fiona.open(os.path.join(root, dataset['filename'])) as fp:\n",
    "            assert fp.crs['init'] == crs, 'dataset is not in the British National Grid projection'\n",
    "            for feature in tqdm(fp, desc=dataset['company']):\n",
    "                # Get the properties that we haven't declared and are thus unexpected. Then complain \n",
    "                # if necessary.\n",
    "                missing = {key: value for key, value in feature['properties'].items() \n",
    "                            if key not in dataset['properties']}\n",
    "                unexpected = set(feature['properties']) - set(dataset['properties'])\n",
    "                if unexpected:\n",
    "                    raise KeyError(f'could not remap unexpected fields {unexpected}')\n",
    "                missing = set(dataset['properties']) - set(feature['properties'])\n",
    "                if missing:\n",
    "                    raise KeyError(f'missing expected fields {missing}')\n",
    "\n",
    "                # Create a record for the consolidated dataset.\n",
    "                properties = {\n",
    "                    'company': dataset['company'],\n",
    "                    'geometry': sh.geometry.shape(feature['geometry']),\n",
    "                }\n",
    "\n",
    "                # Process each feature and either store it if desired or apply a function to \n",
    "                # validate the data.\n",
    "                for key, value in feature['properties'].items():\n",
    "                    key = dataset['properties'][key]\n",
    "                    if callable(key):\n",
    "                        key(properties, value)\n",
    "                    elif key:\n",
    "                        properties[key] = value\n",
    "                        \n",
    "                # Check that we have *some* identifier.\n",
    "                if not any([properties.get('reference'), properties.get('name')]):\n",
    "                    logging.warning(f'{dataset[\"company\"]} has feature without identifiers')\n",
    "\n",
    "                # Store the record for future processing.\n",
    "                company_catchments.append(properties)\n",
    "\n",
    "        # Package the features in a geopandas frame and apply postprocessing steps if necessary.\n",
    "        company_catchments = gpd.GeoDataFrame(company_catchments, crs=crs)\n",
    "        total += len(company_catchments)\n",
    "        postprocessor = dataset.get('postprocessor')\n",
    "        if postprocessor:\n",
    "            company_catchments = postprocessor(company_catchments)\n",
    "            \n",
    "        parts.append(company_catchments)\n",
    "    except Exception as ex:\n",
    "        raise RuntimeError(f'failed to process {dataset[\"filename\"]}') from ex\n",
    "        \n",
    "catchments_raw = pd.concat(parts, ignore_index=True)\n",
    "print(f'loaded {len(catchments_raw)} catchments with {total} shapes')\n",
    "\n",
    "# Filter out catchments for Severn Trent that are odd.\n",
    "fltr = (catchments_raw.company == 'severn_trent_water') & (\n",
    "    (catchments_raw.reference == '0')  # These are private treatment works.\n",
    "    | catchments_raw.reference.str.startswith('TW')  # These are from Thames Water ...\n",
    "    | catchments_raw.reference.str.startswith('AW')  # ... Anglian Water ...\n",
    "    | catchments_raw.reference.str.startswith('UU')  # ... United Utilities ...\n",
    "    | catchments_raw.reference.str.startswith('YW')  # ... Yorkshire Water ...\n",
    "    | catchments_raw.reference.str.startswith('WW')  # ... Welsh Water ...\n",
    ")\n",
    "catchments_raw = catchments_raw[~fltr]\n",
    "print(f'dropped {fltr.sum()} treatment works not owned by Severn Trent Water')\n",
    "\n",
    "# Backfill the name if only the ref is available and vice-versa.\n",
    "catchments_raw.name.fillna(catchments_raw.reference, inplace=True)\n",
    "catchments_raw.reference.fillna(catchments_raw.name, inplace=True)\n",
    "\n",
    "# Buffer to make the shapes valid.\n",
    "print(f'{100 * catchments_raw.is_valid.mean():.3f}% of catchments are valid')\n",
    "catchments_raw['geometry'] = catchments_raw.geometry.buffer(0)\n",
    "\n",
    "# Show the breakdown by company.\n",
    "counts = catchments_raw.groupby('company').count()\n",
    "print(counts.sum())\n",
    "counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bottom-roads",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply manual fixes to the data.\n",
    "update = [\n",
    "    [('thames_water', 'EAST HYDE STW'),\n",
    "     {'comment': 'Overlap with the Chalton and Offley catchments from Anglian Water.'}],\n",
    "    [('thames_water', 'NAGS HEAD LANE STW'),\n",
    "     {'comment': 'Overlap with the (Shenfield and Hutton) and Upminster catchments from Anglian Water.'}],\n",
    "    [('thames_water', 'BARKWAY STW'),\n",
    "     {'comment': 'The catchment as substantial overlap with the Barley catchment from anglian_water.'}],\n",
    "    [('anglian_water', 'Chalton'),\n",
    "     {'comment': 'Overlap with the East Hyde catchment from Thames Water.'}],\n",
    "    [('anglian_water', 'Shenfield And Hutton'),\n",
    "     {'comment': 'Overlap with the Nags Head Lane catchment from Thames Water.'}],\n",
    "    [('anglian_water', 'Buckminster'),\n",
    "     {'comment': 'Overlap with the Waltham catchment from Severn Trent Water.'}],\n",
    "    [('anglian_water', 'Barley'),\n",
    "     {'comment': 'Overlap with the BARKWAY STW catchment from thames_water.'}],\n",
    "    [('severn_trent_water', 'CLAY CROSS (WRW)'),\n",
    "     {'comment': 'Overlap with the Danesmoor catchment from Yorkshire Water.'}],\n",
    "    [('severn_trent_water', 'WALTHAM (WRW)'),\n",
    "     {'comment': 'Overlap with the Buckminster catchment from Anglian Water.'}],\n",
    "    [('severn_trent_water', 'ABBEY LATHE - MALTBY (WRW)'),\n",
    "     {'comment': 'Overlap with the Aldwarke catchment from Yorkshire Water.'}],\n",
    "    [('yorkshire_water', 'Danesmoor'),\n",
    "     {'comment': 'Overlap with the Clay Cross catchment from Severn Trent Water.'}],\n",
    "    [('yorkshire_water', 'Aldwarke'),\n",
    "     {'comment': 'Overlap with the Abbey Lathe - Maltby catchment from Severn Trent Water.'}],\n",
    "    [('severn_trent_water', 'WIGMORE (WRW)'), \n",
    "     {'comment': 'Overlap with the The Orchards SPS catchment from welsh_water.'}],\n",
    "    [('welsh_water', 'The Orchards SPS'), \n",
    "     {'comment': 'Overlap with the WIGMORE (WRW) catchment from severn_trent_water.'}],\n",
    "    [('united_utilities', 'UU-08-SC42-BROMB'), \n",
    "     {'comment': 'Overlap with the NESTON catchment from welsh_water.'}],\n",
    "    [('welsh_water', 'NESTON'), \n",
    "     {'comment': 'Overlap with the UU-08-SC42-BROMB catchment from united_utilities.'}],\n",
    "    [('severn_trent_water', 'SCARCLIFFE (WRW)'), \n",
    "     {'comment': 'Overlap with the Bolsover catchment from yorkshire_water.'}],\n",
    "    [('yorkshire_water', 'Bolsover'), \n",
    "     {'comment': 'Overlap with the SCARCLIFFE (WRW) catchment from severn_trent_water.'}],\n",
    "    [('severn_trent_water', 'CLOWNE (WRW)'), \n",
    "     {'comment': 'Overlap with the Staveley catchment from yorkshire_water.'}],\n",
    "    [('yorkshire_water', 'Staveley'), \n",
    "     {'comment': 'Overlap with the CLOWNE (WRW) catchment from severn_trent_water.'}],\n",
    "    [('southern_water', 'LONGFIELD'), \n",
    "     {'comment': 'Overlap with the LONG REACH STW catchment from thames_water.'}],\n",
    "    [('united_utilities', 'UU-05-SC26-COLNE'), \n",
    "     {'comment': 'Overlap with the Foulridge catchment from yorkshire_water.'}],\n",
    "    [('yorkshire_water', 'Foulridge'), \n",
    "     {'comment': 'Overlap with the UU-05-SC26-COLNE catchment from united_utilities.'}],\n",
    "    [('anglian_water', 'Offley'), \n",
    "     {'comment': 'Overlap with the EAST HYDE STW catchment from thames_water.'}],\n",
    "    [('thames_water', 'RIVERSIDE STW'), \n",
    "     {'comment': 'Overlap with the Upminster catchment from anglian_water.'}],\n",
    "    [('anglian_water', 'Upminster'), \n",
    "     {'comment': 'Overlap with the RIVERSIDE STW and NAGS HEAD LANE STW catchments from thames_water.'}],\n",
    "    [('southern_water', 'REDLYNCH'), \n",
    "     {'comment': 'Overlap with the DOWNTON STW CATCHMENT catchment from wessex_water.'}],\n",
    "    [('wessex_water', 'DOWNTON STW CATCHMENT'), \n",
    "     {'comment': 'Overlap with the REDLYNCH catchment from southern_water.'}],\n",
    "    [('severn_trent_water', 'COALEY (WRW)'), \n",
    "     {'comment': 'Overlap with the NORTH NIBLEY STW CATCHMENT catchment from wessex_water.'}],\n",
    "    [('wessex_water', 'NORTH NIBLEY STW CATCHMENT'), \n",
    "     {'comment': 'Overlap with the COALEY (WRW) catchment from severn_trent_water.'}],\n",
    "    [('united_utilities', 'UU-11-SC58-ELLES'), \n",
    "     {'comment': 'Overlap with the CHESTER catchment from welsh_water.'}],\n",
    "    [('welsh_water', 'CHESTER'), \n",
    "     {'comment': 'Overlap with the UU-11-SC58-ELLES and UU-11-SC58-WAVER catchments from united_utilities.'}],\n",
    "    [('severn_trent_water', 'LYDNEY (WRW)'), \n",
    "     {'comment': 'Overlap with the NEWLAND catchment from welsh_water.'}],\n",
    "    [('welsh_water', 'NEWLAND'), \n",
    "     {'comment': 'Overlap with the LYDNEY (WRW) catchment from severn_trent_water.'}],\n",
    "    [('thames_water', 'GUILDFORD STW'), \n",
    "     {'comment': 'Overlap with the GUILDFORD catchment from southern_water.'}],\n",
    "    [('southern_water', 'GUILDFORD'), \n",
    "     {'comment': 'Overlap with the GUILDFORD STW catchment from thames_water.'}],\n",
    "    [('severn_trent_water', 'STRONGFORD (WRW)'), \n",
    "     {'comment': 'Overlap with the UU-11-SC59-MADEL and UU-11-SC60-KIDSG catchments from united_utilities.'}],\n",
    "    [('united_utilities', 'UU-11-SC59-MADEL'), \n",
    "     {'comment': 'Overlap with the STRONGFORD (WRW) catchment from severn_trent_water.'}],\n",
    "    [('southern_water', 'PENNINGTON'), \n",
    "     {'comment': 'Overlap with the CHRISTCHURCH STW CATCHMENT catchment from wessex_water.'}],\n",
    "    [('wessex_water', 'CHRISTCHURCH STW CATCHMENT'), \n",
    "     {'comment': 'Overlap with the PENNINGTON catchment from southern_water.'}],\n",
    "    [('thames_water', 'LONG REACH STW'), \n",
    "     {'comment': 'Overlap with the IDE HILL TO THAMES, LONG HILL TO THAMES, LONGFIELD, and NORTHFLEET catchments from southern_water.'}],\n",
    "    [('southern_water', 'IDE HILL TO THAMES'), \n",
    "     {'comment': 'Overlap with the LONG REACH STW catchment from thames_water.'}],\n",
    "    [('united_utilities', 'UU-08-SC41-BIRKE'), \n",
    "     {'comment': 'Overlap with the HESWALL catchment from welsh_water.'}],\n",
    "    [('welsh_water', 'HESWALL'), \n",
    "     {'comment': 'Overlap with the UU-08-SC41-BIRKE catchment from united_utilities.'}],\n",
    "    [('severn_trent_water', 'BRANTON (WRW)'), \n",
    "     {'comment': 'Overlap with the Sandall catchment from yorkshire_water.'}],\n",
    "    [('yorkshire_water', 'Sandall'), \n",
    "     {'comment': 'Overlap with the BRANTON (WRW) catchment from severn_trent_water.'}],\n",
    "    [('united_utilities', 'UU-11-SC60-KIDSG'), \n",
    "     {'comment': 'Overlap with the STRONGFORD (WRW) catchment from severn_trent_water.'}],\n",
    "    [('united_utilities', 'UU-11-SC58-WAVER'), \n",
    "     {'comment': 'Overlap with the CHESTER catchment from welsh_water.'}],\n",
    "    [('severn_trent_water', 'DINNINGTON (WRW)'), \n",
    "     {'comment': 'Overlap with the Woodhouse Mill catchment from yorkshire_water.'}],\n",
    "    [('yorkshire_water', 'Woodhouse Mill'), \n",
    "     {'comment': 'Overlap with the DINNINGTON (WRW) catchment from severn_trent_water.'}],\n",
    "    [('southern_water', 'LONGFIELD HILL TO THAMES'), \n",
    "     {'comment': 'Overlap with the LONG REACH STW catchment from thames_water.'}],\n",
    "    [('southern_water', 'NORTHFLEET'), \n",
    "     {'comment': 'Overlap with the LONG REACH STW catchment from thames_water.'}],\n",
    "]\n",
    "\n",
    "# Treatment works to remove because they are duplicates. We can check which is the correct company \n",
    "# using https://www.water.org.uk/advice-for-customers/find-your-supplier/.\n",
    "drop = {\n",
    "    'northumbrian_water': [\n",
    "        # Great Smeaton, DL6 2ET, supplied Yorkshire Water.\n",
    "        'GREAT SMEATON STW NZ34043201',\n",
    "        # Brampton, CA8 2Qg, United Utilities.\n",
    "        'TINDALE STW NY61599206',\n",
    "    ],\n",
    "    'thames_water': [\n",
    "        # Wingrave, HP22 4PS, supplied by Anglian Water.\n",
    "        'WINGRAVE STW',\n",
    "        # Stewkley, LU7 0EL, supplied by Anglian Water.\n",
    "        'STEWKLEY STW',\n",
    "        # RH12 4BB, supplied by Southern Water.\n",
    "        'COLGATE STW',\n",
    "        # Anglian.\n",
    "        'DUNSTABLE (AW) STW',\n",
    "    ],\n",
    "    'united_utilities': [\n",
    "        # Doveholes, SK17 8BL, supplied by Severn Trent.\n",
    "        'UU-10-SC54-DOVEH',\n",
    "        # Various catchments supplied by Severn Trent.\n",
    "        'UU-11-SC59-SVNTR',\n",
    "    ],\n",
    "    'severn_trent_water': [\n",
    "        # Acton Green, SWR6 5AA, supplied by Welsh Water.\n",
    "        'ACTON GREEN (WRW)',\n",
    "        # Severn Trent are missing Thealby and Normanby \n",
    "        # (but the Anglian Water catchment is rather coarse).\n",
    "        'BURTON STATHER (WRW)',\n",
    "    ],\n",
    "    'anglian_water': [\n",
    "        # SG4 7AA, Thames Water.\n",
    "        'Weston',\n",
    "        # CM24 8BE, Thames Water.\n",
    "        'Henham',\n",
    "        # DN9 3ED, Severn Trent.\n",
    "        'Ex Raf Finningley',\n",
    "        # Charndon, Thames.\n",
    "        'Chardon',\n",
    "        # This is technically Anglian, but Thames have better data.\n",
    "        # Specifically, Anglian are missing Marston St Lawrence.\n",
    "        'Halse',\n",
    "        # Severn Trent have better data. Anglian are missing Bothamsall.\n",
    "        'Elkesley',\n",
    "        # LE15 7PL, Severn Trent.\n",
    "        'Market Overton',\n",
    "        # The Severn Trent data is more refined and actually covers Alkborough.\n",
    "        'Alkborough',\n",
    "        # Covered by Thames Water.\n",
    "        'Graveley',\n",
    "    ],\n",
    "    'welsh_water': [\n",
    "        # SY12 9LD, Severn Trent Water.\n",
    "        'DUDLESTON HEATH',\n",
    "        # sy13 2af, Severn Trent.\n",
    "        'ASH MAGNA',\n",
    "    ],\n",
    "    'southern_water': [\n",
    "        # Thames Water provide the treatment in Crawley.\n",
    "        'COPTHORNE',\n",
    "        # Covered by Thames.\n",
    "        'SMALLFIELD',\n",
    "        'HASLEMERE',\n",
    "        'GREENHYTHE',\n",
    "    ]\n",
    "}\n",
    "\n",
    "catchments = catchments_raw.copy()\n",
    "\n",
    "# Drop each catchment individually, ensuring that we drop exactly one per name/company.\n",
    "for company, names in drop.items():\n",
    "    for name in names:\n",
    "        f = (catchments.name == name) & (catchments.company == company)\n",
    "        assert f.sum() == 1, (company, name, f.sum())\n",
    "        catchments = catchments[~f]\n",
    "print(f'dropped {sum(len(names) for names in drop.values())} catchments')\n",
    "\n",
    "# Update each record in turn, ensuring we only update exactly one per name/company(/reference).\n",
    "for key, item in update:\n",
    "    if len(key) == 2:\n",
    "        (company, name) = key\n",
    "        f = (catchments.name == name) \\\n",
    "            & (catchments.company == company)\n",
    "    else:\n",
    "        (company, name, reference) = key\n",
    "        f = (catchments.name == name) \\\n",
    "            & (catchments.company == company) \\\n",
    "            & (catchments.reference == reference)\n",
    "    assert f.sum() == 1, (company, name, f.sum())\n",
    "\n",
    "    for key, value in item.items():\n",
    "        catchments.loc[f, key] = value\n",
    "\n",
    "print(f'updated {len(update)} catchments')\n",
    "\n",
    "# Show the counts by company after updates.\n",
    "counts = catchments.groupby('company').count()\n",
    "print(counts.sum())\n",
    "counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "outside-arabic",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show an overview of the data we are now left with.\n",
    "ax = catchments.plot(\n",
    "    column='company', \n",
    "    categorical=True, \n",
    "    legend=True,\n",
    "    legend_kwds={\n",
    "        'fontsize': 'x-small',\n",
    "        'ncol': 2,\n",
    "    }\n",
    ")\n",
    "ax.ticklabel_format(scilimits=(0, 0), useMathText=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "broke-highway",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the intersections of catchments.\n",
    "sindex = catchments.geometry.sindex\n",
    "indices = sindex.query_bulk(catchments.geometry).T\n",
    "print(f'found {len(indices)} intersections')\n",
    "\n",
    "# Remove self-intersections and only report each intersection once.\n",
    "indices = indices[indices[:, 0] < indices[:, 1]]\n",
    "print(f'found {len(indices)} intersections without self intersections')\n",
    "\n",
    "# Calculate the actual intersection areas and remove zero-area intersections.\n",
    "i, j = indices.T\n",
    "x = catchments.iloc[i].geometry.reset_index(drop=True)\n",
    "y = catchments.iloc[j].geometry.reset_index(drop=True)\n",
    "intersection_areas = x.intersection(y).area.values\n",
    "f = intersection_areas > 0\n",
    "intersection_areas = intersection_areas[f]\n",
    "indices = indices[f]\n",
    "\n",
    "print(f'found {len(indices)} intersections with non-zero intersection area')\n",
    "\n",
    "# Evaluate the unions.\n",
    "i, j = indices.T\n",
    "x = catchments.iloc[i].geometry.reset_index(drop=True)\n",
    "y = catchments.iloc[j].geometry.reset_index(drop=True)\n",
    "union_areas = x.union(y).area.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "boolean-front",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a data frame that has information about both pairs of an intersection as well as \n",
    "# intersection information, such as the area.\n",
    "parts = [catchments.iloc[i].copy().reset_index(drop=True), catchments.iloc[j].copy().reset_index(drop=True)]\n",
    "for suffix, part in zip(['_x', '_y'], parts):\n",
    "    part.columns = [col + suffix for col in part]\n",
    "overlapping = gpd.GeoDataFrame(pd.concat(parts, axis=1))\n",
    "overlapping['intersection_areas'] = intersection_areas\n",
    "overlapping['union_areas'] = union_areas\n",
    "overlapping['iou'] = overlapping.intersection_areas / overlapping.union_areas\n",
    "overlapping = overlapping.sort_values('iou', ascending=False)\n",
    "\n",
    "# Restrict to settings where the water companies that provided the data differ (e.g. to avoid self \n",
    "# overlap within subcatchments).\n",
    "f = (overlapping.company_y != overlapping.company_x)\n",
    "print(f'{f.sum()} overlapping catchments from different companies')\n",
    "\n",
    "# Filter out the pairs where one or more catchments have an annotation.\n",
    "f &= (overlapping.comment_x.isnull() | overlapping.comment_y.isnull())\n",
    "overlapping = overlapping.loc[f]\n",
    "\n",
    "# Manual exclusion of items that have been investigated and considered insignificant.\n",
    "ignore = [\n",
    "    {\n",
    "        'company_x': 'severn_trent_water',\n",
    "        'name_x': 'private',\n",
    "        'company_y': 'welsh_water',\n",
    "        'name_y': 'LYDBROOK',\n",
    "    }\n",
    "]\n",
    "for item in ignore:\n",
    "    f = True\n",
    "    for key, value in item.items():\n",
    "        f = f & (overlapping[key] == value)\n",
    "    overlapping = overlapping.loc[~f]\n",
    "\n",
    "assert len(overlapping) == 0, f'{len(overlapping)} unaccounted overlaps found'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "approximate-suicide",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the breakdown of areas covered by each company in square kilometres.\n",
    "areas = catchments.groupby('company').geometry.apply(lambda x: x.area.sum()) / 1e6\n",
    "print(f'total area covered: {areas.sum()}')\n",
    "areas.round()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "shaped-baking",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show an overview figure.\n",
    "fig = plt.figure(figsize=(5.8, 5.8))\n",
    "gs = fig.add_gridspec(2, 2, width_ratios=[3, 2])\n",
    "ax3 = fig.add_subplot(gs[:, 0])\n",
    "tolerance = 1000\n",
    "countries = gpd.read_file('data/geoportal.statistics.gov.uk/countries20_BGC.zip')\n",
    "countries.simplify(tolerance=tolerance).plot(\n",
    "    ax=ax3, facecolor='none', edgecolor='k', lw=.5)\n",
    "simplified = catchments.copy()\n",
    "simplified['geometry'] = simplified.simplify(tolerance=tolerance)\n",
    "simplified.plot(\n",
    "    ax=ax3, column='company', zorder=9, edgecolor='face', lw=.5)\n",
    "ax3.set_ylim(top=1.5e6)\n",
    "\n",
    "ax1 = fig.add_subplot(gs[0, 1])\n",
    "ax2 = fig.add_subplot(gs[1, 1])\n",
    "\n",
    "alpha = 0.5\n",
    "\n",
    "# Easy to tell that they're duplicates.\n",
    "colors = ['C0', 'C3']\n",
    "subset = catchments_raw.loc[catchments_raw.name.str.lower().str.contains('market overton').fillna(False)]\n",
    "subset.plot(ax=ax1, facecolor=colors, alpha=alpha, edgecolor='face')\n",
    "handles = [mpl.patches.Rectangle((0, 0), 1, 1, color=color, alpha=alpha) \n",
    "           for color in colors]\n",
    "\n",
    "# Hard to tell what's going on.\n",
    "colors = ['C5', 'C0']\n",
    "subset = catchments_raw.loc[np.in1d(catchments_raw.name, ['EAST HYDE STW', 'Chalton'])]\n",
    "subset.plot(ax=ax2, facecolor=colors, alpha=alpha, edgecolor='face')\n",
    "\n",
    "for ax, label in zip([ax3, ax1, ax2], ['(a)', '(b)', '(c)']):\n",
    "    ax.set_axis_off()\n",
    "    ax.text(0.05, 0.05, label, transform=ax.transAxes)\n",
    "    \n",
    "labels = names = [x.replace('_', ' ').title() for x \n",
    "         in sorted(catchments.company.unique())]\n",
    "handles = [mpl.patches.Rectangle((0, 0), 1, 1, color=f'C{i}') \n",
    "           for i in range(len(labels))]\n",
    "fig.legend(handles, labels, ncol=2, loc='upper left')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('overview.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hydraulic-contract",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate unique identifiers for the catchments based on the first ten characters of the sha1 hash\n",
    "# of the geometry rounded to the nearest metre. Then save the file.\n",
    "def hash_function(geometry):\n",
    "    coords = geometry.envelope.exterior.coords\n",
    "    return hashlib.sha1(np.round(coords)).hexdigest()[:10]\n",
    "\n",
    "identifiers = catchments.geometry.apply(hash_function)\n",
    "assert identifiers.nunique() == len(identifiers)\n",
    "\n",
    "catchments['identifier'] = identifiers\n",
    "catchments.to_file('data/wastewater_catchment_areas_public/catchments_consolidated.shp')"
   ]
  }
 ],
 "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
}