<jupyter_start><jupyter_text>UK rivers catchment mapInspired by the really cool catchments maps they make at Grasshopper Geography:  https://www.grasshoppergeography.com/River-Maps/  I wanted to see if I could come up with something along the same lines. In the UK the Ordnance Survey has produced and made publicly available a shapefile for all watercourses of the UK except Northern Ireland. It can be found here:  https://www.ordnancesurvey.co.uk/business-government/products/open-map-rivers  Import modules<jupyter_code># Import required modules
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
# Read in watercourse shp file
gdf = gpd.read_file(r'C:\Users\sgibs\Documents\Datasets\os_open_rivers\WatercourseLink.shp')
# What do we have?
print(gdf.shape)
print(gdf.crs)
gdf.head(2)<jupyter_output>(187356, 10)
epsg:27700<jupyter_text>* Gdf contains many ~190,000 individual linestring geometries* Watercourses are made up of a series of linestrings* _Some_ linestrings are labelled with river name, but not many* Majority are of form 'inlandRiver', but also includes tidal river, canal and lake<jupyter_code># Interested in catchment watercourses so drop canals
gdf = gdf.drop(gdf[gdf['form']=='canal'].index)
print(gdf.shape)<jupyter_output>(186100, 10)<jupyter_text>Including tidal rivers has the impact of merging many different catchments into a much smaller number. For instance combining the Trent and the Ouse together at the Humber. Although removing the tidal rivers can also obscure things in a few locations by preventing a number of catchments from joining<jupyter_code># After testing with and without, drop tidal rivers
gdf = gdf.drop(gdf[gdf['form']=='tidalRiver'].index)
print(gdf.shape)<jupyter_output>(175877, 10)<jupyter_text>One of the main difficulties I found with my process was the long time needed to run a number of my functions due to the very large dataset involved. It was much easier to test using a subset<jupyter_code># Take a subset of dataset using geo info
wales = gdf.cx[160000:360000,160000:390000]
# How much smaller is our subset?
print(wales.shape)<jupyter_output>(27010, 10)<jupyter_text>Define processing functions Initially tried to merge touching linestrings to form overall catchments. However there are gaps or missing geometry, which resulted in very few connected geometries. My solution was to buffer the linestrings and merge the resulting polygons. The increased thickness produced overlapping geometries which could be merged.<jupyter_code>def bufferUnion(gdf, size=50):
    """
    buffers gdf geometries and merges all geometries
    prints original number of linestrings and number of new merged polygons
    
    Parameters
    ----------
    gdf : geodataframe
        gdf containing catchment geometries
    
    size : int
        buffer thickness
        
    Returns
    -------
    union : 
    """
    gdfBuffer = gdf.buffer(size)
    union = gdfBuffer.unary_union
    print('{} elements reduced to {}'.format(gdf.shape[0], len(union)))
    return union    
def rebuildGeo(union, crs=gdf.crs):
    """
    creates a geodataframe from single geometry object
    splits individual geometries (multi polygons) into individual catchment objects (merged polygons)
    reassigns a crs to geometries
    
    Parameters
    ----------
    union : single geometry object
        object contains multiple merged overlapping polygons
        
    crs : coordinate reference system string identifier
        default = gdf.crs
        crs to rebuild gdf geometry
        
    Returns
    -------
    gdf : geodataframe
        contains rebuilt geometries per row
    """
    series = pd.Series(union)
    gdf = gpd.GeoDataFrame(series)
    gdf = gdf.rename({0:'geometry'}, axis =1)
    gdf = gdf.set_geometry(col='geometry')
    gdf = gdf.reset_index(drop=True)
    gdf = gdf.set_crs(crs)
    return gdf    
def bufferPlot(gdf, size=200):
    """
    buffers gdf geometry by size for use when plotting
    
    Parameters
    ----------
    gdf : geodataframe
        gdf containing geometry to buffer
        
    size : int
        default = 200
        size to buffer geometry by        
        
    Returns
    -------
    gdf : geodataframe
        contains buffered geometries
    """
    gdf = gdf.buffer(size)
    return gdf
def addHullsAreaCentrGeo(gdf):
    """
    adds convex hull and centroid geometries to gdf
    calculates convex hull area
    
    Parameters
    ----------
    gdf : geodataframe
        gdf to add columns to
        
    Returns
    -------
    gdf : geodataframe
        gdf containing additional columns
    """
    gdf['hullGeometry'] = gdf['geometry'].convex_hull
    gdf['hullArea'] = gdf.set_geometry(col='hullGeometry').area
    gdf['centroidGeometry'] = gdf.set_geometry(col='hullGeometry').centroid
    return gdf
def areaFilter(gdf, areaThreshold = 1e7):
    """
    removes rows from gdf with convex hull geometries of area less than the area threshold
    
    Parameters
    ----------
    gdf : geodataframe
        gdf from which to remove geometries
    
    areaThreshold : int
        default = 1e7
        minimum area for convex hull geometries in returned gdf
        
    Returns
    -------
    gdf : geodataframe
        original gdf with convex hull geometry areas above areaThreshold
    """
    indexDrop = gdf[gdf['hullArea'] < areaThreshold].index
    print('For area threshold of {}m^2:\n{} rows remaining from {}'.format(areaThreshold,len(gdf)-len(indexDrop),len(gdf)))
    return gdf.drop(indexDrop)<jupyter_output><empty_output><jupyter_text>Processing Care - long time to execute depending on size of input<jupyter_code># Create union object
union = bufferUnion(wales)<jupyter_output>27010 elements reduced to 765<jupyter_text>Care - long time to execute:<jupyter_code># Rebuild geometry
gdfUnion = rebuildGeo(union)
# Save merged and rebuilt geometry to file
#gdfUnion.to_file('OSOpenRivers/catchments.shp', driver = 'ESRI Shapefile', encoding='utf-8')
# Add convex hull and centroid geometries and calculate convex hull area variable
gdfPlot = addHullsAreaCentrGeo(gdfUnion)
# Use convex hull areas to remove small catchments
gdfPlot = areaFilter(gdfPlot, areaThreshold = 1e6)<jupyter_output>For area threshold of 1000000.0m^2:
301 rows remaining from 765<jupyter_text>Care - long time to execute:<jupyter_code># Buffer catchment geometry for plotting
gdfPlot['geometryBuff'] = bufferPlot(gdfPlot)<jupyter_output><empty_output><jupyter_text>A key difficulty when plotting the final catchments was the issue of adjacent geometries plotting in the same colour. To try and address this I labelled each catchment with random integer 0-9 and plotted this as a categorical variable. There is still a significant chance of adjacent catchments sharing the same random label.<jupyter_code># Assign random integer label to each row
gdfPlot['rand'] = np.random.randint(0,10,size=len(gdfPlot))
# Figure for plotting catchment map
fig, ax = plt.subplots(figsize=(8,8))

#gdfPlot.set_geometry('centroidGeometry').plot(ax=ax,
#gdfPlot.set_geometry('hullGeometry').plot(ax=ax,
gdfPlot.set_geometry('geometryBuff').plot(ax=ax,
                  column='rand',
                  #categorical=True,
                  cmap='tab10')

# Remove map axes
ax.set_axis_off()

# Save file, high dpi for printing
plt.savefig('OSOpenRivers/catchments.png', dpi=500)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08440860215053764, 'fraction_numerical': 0.016263440860215053, 'mean_word_length': 3.659361302442079, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 31, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '274855', 'license': 'MIT', 'n_tokens_mistral': 2352, 'n_tokens_neox': 2159, 'n_words': 744, 'path': 'notebooks/OSOpenRivers.ipynb', 'repo_name': 'phtevegibson/data-projects', 'size': '7440'}