Create maps#

Using IVT catchment and ambulance trust shapefiles, create maps - for example, population weighted mean IMD score per IVT catchment, or mean transfer time by ambulance trust.

Notes:

  • When using the 2018 age CSV file, there were two Welsh hospital catchments (CF311RQ Bridgend, CF144XW Cardiff) that each contained two LSOAs that are not present in that age file (Rhondda Cynon Taf 023F, Rhondda Cynon Taf 027F, Rhondda Cynon Taf 027G, Rhondda Cynon Taf 023G). Since we use the population count from this age csv file to calculate a weighted sum for the IVT catchment, these LSOAs were previously (1) excluded from the weighted mean, and (2) excluded from the map. However, in the 2020 age CSV file now used, these areas have data, so these steps are no longer required and have been removed.

  • To help with geopandas use MichaelAllen1966/2010_geopandas

  • Uses MatPlotLib for mapping the data since that gives us more power than the GeoPandas DataFrame plot method.

1 Set up#

1.1 Import libraries and define file paths#

import contextily as ctx
from dataclasses import dataclass
import geopandas as gpd
from mapclassify import greedy
import matplotlib.pyplot as plt
from matplotlib.transforms import Bbox
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
import numpy as np
import os
import pandas as pd
from pandas.api.types import is_numeric_dtype

# Linting
%load_ext pycodestyle_magic
%pycodestyle_on
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Input In [1], in <cell line: 1>()
----> 1 import contextily as ctx
      2 from dataclasses import dataclass
      3 import geopandas as gpd

ModuleNotFoundError: No module named 'contextily'
# Define file paths
@dataclass(frozen=True)
class Paths:
    '''Singleton object for storing paths to data and database.'''

    data = './data'
    collated = 'collated_data_amb.csv'

    hospitals = 'hospitals'
    stroke_hospitals = 'stroke_hospitals_2022.csv'

    shapefiles = 'shapefiles'
    ivt_catchment = 'ivt_catchment.shp'
    amb_catchment = 'amb_catchment.shp'
    countries = 'Countries_(December_2022)_GB_BFC.zip'


paths = Paths()

1.2 Load data#

We convert data to epsg:3857, as that is needed when using base maps. Base maps are a reference map that you can then overlay your data onto - for example, a topographic basemap showing elevation. We can set the crs (Coordinate Reference System) when loading a GeoPandas DataFrame. CRS to consider:

  • EPSG:27700 – BNG (British National Grid Eastings and Northings).

  • EPSG:27700 OSGB 1936 / British National Grid – United Kingdom Ordnance Survey. Co-ordinates are in Eastings (X) and Northings (Y).

  • EPSG:4326 WGS 84 – WGS84 - World Geodetic System 1984. Co-ordinates are in Longitude (X) and Latitude (Y).

  • EPSG:3857 - projection for displaying lat/long as a flat map

IVT catchment shapefile - shape file of the catchments for the IVT units in England and Wales, which we load into a GeoPandas DataFrame. Note that we can set the crs (Coordinate Reference System) when loading a GeoPandas DataFrame. EPSG:27700 is the crs to use when geography is in BNG (British National Grid Eastings and Northings). Note: If this fails to run, you need to run 02_create_ivt_catchment_shapefile.ipynb to generate the IVT catchment shapefile.

# Load IVT catchment shapefile
gdf_ivt_catchment = gpd.read_file(
    os.path.join(paths.data, paths.shapefiles, paths.ivt_catchment),
    crs='EPSG:27700')

# Set index
gdf_ivt_catchment.set_index(['closest_iv'], inplace=True)

# Preview the dataframe
display(gdf_ivt_catchment.head(2))

# Plot the shapefile
gdf_ivt_catchment.plot(figsize=(6, 6))
plt.show()

# Convert crs to epsg 3857
gdf_ivt_catchment = gdf_ivt_catchment.to_crs(epsg=3857)
LSOA11CD LSOA11NM LSOA11NMW geometry
closest_iv
B152TH E01008881 Birmingham 067A Birmingham 067A MULTIPOLYGON (((403011.723 266679.152, 402597....
B714HJ E01008899 Birmingham 037A Birmingham 037A MULTIPOLYGON (((395950.714 272162.470, 395950....
_images/afc10f8220f0c3fd5783b9aa69ec469126718223b45a32ff9625d7e73b2cc26b.png

Country outline - used to produce outline of countries on the map.

# Load country outline
outline = gpd.read_file(os.path.join(
    paths.data, paths.shapefiles, paths.countries))

# Drop Scotland
mask = outline['CTRY22NM'] != 'Scotland'
outline = outline[mask]

# Convert to epsg:3857
outline = outline.to_crs(epsg=3857)

# Preview the dataframe
display(outline.head())

# Plot the shapefile
outline.plot(figsize=(6, 6))
plt.show()
FID CTRY22CD CTRY22NM CTRY22NMW BNG_E BNG_N LONG LAT GlobalID SHAPE_Leng SHAPE_Area geometry
0 1 E92000001 England Lloegr 394883 370883 -2.07811 53.23497 {096C3540-B006-4D2D-BCDF-B3E945B80D51} 1.358939e+07 1.304623e+11 MULTIPOLYGON (((-712384.210 6422975.926, -7123...
2 3 W92000004 Wales Cymru 263405 242881 -3.99417 52.06741 {1DFA050A-C94B-4E58-8210-DE8274C7C33E} 3.273688e+06 2.078260e+10 MULTIPOLYGON (((-347425.136 6688819.680, -3474...
_images/0245cd241975c5c449acf2768ac9dcbd0a12e146fedd68197d6ef0267276be99.png

LSOA information - load collated_data.csv produced by 01_combine_demographic_data.ipynb which has information on each LSOA.

df_lsoa = pd.read_csv(os.path.join(paths.data, paths.collated))
df_lsoa.head(2)
LSOA admissions closest_ivt_unit closest_ivt_unit_time closest_mt_unit closest_mt_unit_time closest_mt_transfer closest_mt_transfer_time total_mt_time ivt_rate ... age_band_males_65 age_band_males_70 age_band_males_75 age_band_males_80 age_band_males_85 age_band_males_90 ambulance_service local_authority_district_22 LAD22NM country
0 Welwyn Hatfield 010F 0.666667 SG14AB 18.7 NW12BU 36.9 CB20QQ 39.1 57.8 6.8 ... 33 28 26 14 5 3 East of England Welwyn Hatfield Welwyn Hatfield England
1 Welwyn Hatfield 012A 4.000000 SG14AB 19.8 NW12BU 36.9 CB20QQ 39.1 58.9 6.8 ... 24 18 21 12 5 4 East of England Welwyn Hatfield Welwyn Hatfield England

2 rows × 115 columns

Hospital data - load data on each hospital, with aim of getting locations of hospitals that provide IVT.

# Read in data on each hospital
gdf_units = gpd.read_file(
    os.path.join(paths.data, paths.hospitals, paths.stroke_hospitals))

# Combine get geometry from Easting and Northing columns
gdf_units["geometry"] = gpd.points_from_xy(
        gdf_units.Easting, gdf_units.Northing)
gdf_units = gdf_units.set_crs(epsg=27700)

# Restrict to the units delivering thrombolysis
mask = gdf_units['Use_IVT'] == '1'
gdf_units = gdf_units[mask]

# Convert crs to epsg 3857
gdf_units = gdf_units.to_crs(epsg=3857)

# Preview dataframe
gdf_units.head(2)
Postcode Hospital_name Use_IVT Use_MT Use_MSU Country Strategic Clinical Network Health Board / Trust Stroke Team SSNAP name ... ivt_rate Easting Northing long lat Neuroscience 30 England Thrombectomy Example hospital_city Notes geometry
0 RM70AG RM70AG 1 1 1 England London SCN Barking Havering and Redbridge University Hospitals N... Queens Hospital Romford HASU ... 11.9 551118 187780 0.179030640661934 51.5686465521504 1 0 Romford POINT (19929.600 6722503.874)
1 E11BB E11BB 1 1 1 England London SCN Barts Health NHS Trust The Royal London Hospital Royal London Hospital HASU ... 13.4 534829 181798 -0.0581329916047372 51.5190178361295 1 1 Royal London POINT (-6471.335 6713620.606)

2 rows × 22 columns

Ambulance catchment shapefile

# Load IVT catchment shapefile
gdf_amb_catchment = gpd.read_file(
    os.path.join(paths.data, paths.shapefiles, paths.amb_catchment),
    crs='EPSG:27700')

# Set index
gdf_amb_catchment.set_index(['ambulance_'], inplace=True)

# Preview the dataframe
display(gdf_amb_catchment.head(2))

# Plot the shapefile
gdf_amb_catchment.plot(figsize=(6, 6))
plt.show()

# Convert crs to epsg 3857
gdf_amb_catchment = gdf_amb_catchment.to_crs(epsg=3857)
LSOA11CD LSOA11NM LSOA11NMW geometry
ambulance_
East Midlands E01013128 North East Lincolnshire 010A North East Lincolnshire 010A MULTIPOLYGON (((433631.972 296489.554, 433631....
East of England E01015589 Peterborough 004A Peterborough 004A MULTIPOLYGON (((503939.688 193345.500, 503939....
_images/4b7183c4b9444d9c75fed037aac0aa03c251f36050999d0b55a3f092130d7d77.png

1.3 Define functions for calculating weighted mean and proportions#

Group data by specified geography (e.g. closest IVT unit) and calculate the weighted average of specified columns - see stackoverflow tutorial here.

def weighted_av(df, group_by, col, new_col, weight='population_all'):
    '''
    Groupby specified column and then find weighted average of another column
    Returns normal mean and weighted mean
    Inputs:
    - df - dataframe, contains data and columns to do calculation
    - group_by - string, column to group by
    - col - string, column to find average of
    - new_col - string, name of column which will contain weighted mean
    - weight - string, name of column to weight by, default 'population_all'
    '''
    res = (df
           .groupby(group_by)
           .apply(lambda x: pd.Series([np.mean(x[col]),
                                       np.average(x[col], weights=x[weight])],
                                      index=['mean', new_col])))
    return (res)

Group data by specific geography (e.g. ambulance trust) and calculate proportion.

def find_proportion(subgroup_col, overall_col, proportion_col, group_col, df):
    '''
    Groupby a particular geography, and find a proportion from two provided
    columns
    Inputs:
    subgroup_col - string, numerator, has population counts for a subgroup
    overall_col - string, denominator, has population count overall
    proportion_col - string, name for the column produced
    group_col - string, column to groupby when calculate proportions
    df - dataframe, with columns used above
    '''
    # Sum the columns for each group
    res = df.groupby(group_col).agg({
        subgroup_col: 'sum',
        overall_col: 'sum'
    })

    # Find the proportions for each group
    res[proportion_col] = res[subgroup_col] / res[overall_col]

    return res

Find counts of people in each category, then calculate proportion.

def find_count_and_proportion(df_lsoa, area, col, count, prop):
    '''
    Finds count of people in each category for each area (based on LSOA info
    and counts), and then calculates proportion in particular category.
    Requires col to contain 'True' 'False' where 'True' is the proportion you
    want to find (so calculates 'True' / 'True'+'False')
    Inputs:
    - df_lsoa - dataframe with information on each LSOA
    - area - string, column to group areas by
    - col - string, column with data of interest
    - count - string, column with population count
    - prop - string, name of new column with proportion
    '''
    # Find count for each area of people in each category
    counts = (df_lsoa.groupby([area, col])[count].sum())

    # Reformat dataframe
    counts = counts.reset_index(name='n')
    counts = counts.pivot(
        index=area,
        columns=col,
        values='n').reset_index().rename_axis(None, axis=1)

    # Set NaN to 0
    counts = counts.replace(np.nan, 0)

    # Find proportion
    counts[prop] = counts['True'] / (counts['True'] + counts['False'])

    # Set index to the area unit
    counts = counts.set_index(area)

    return (counts)

1.4 Extract boundaries#

Extract the bounds of each IVT catchment.

eng_wales_bounds = gdf_ivt_catchment.bounds
eng_wales_bounds.head(2)
minx miny maxx maxy
closest_iv
B152TH -228159.802962 6.840962e+06 -186486.338536 6.891587e+06
B714HJ -230754.673759 6.862190e+06 -174024.968174 6.926919e+06

1.5 Define mapping functions#

Non-overlapping labels of each stroke unit - see stackoverflow explanation here.

def add_nonoverlapping_text_labels(units, ax, col, y_step=0.05, fontsize=7):
    '''
    Creates non-overlapping text labels for stroke units to the input axis
    Inputs:
    - units - dataframe, contains data on each hospital
    - ax - axis object, to add labels to
    - col - string, label for each position
    - y_step - number, amount to adjust position on y axis by
    - fontsize - number, size of font in labels
    Output:
    - ax - axis object with addition of text labels
    '''
    # Create empty array
    text_rectangles = []

    # Sort labels descending in y axis (gives better results)
    units.loc[:, 'sort_by'] = units.geometry.y
    units = units.sort_values(by='sort_by', ascending=False, axis=0)
    del units['sort_by']

    # Add labels to the plot
    # Loop through each stroke unit - in each loop, get x, y and label
    geom_unit = zip(units.geometry.x, units.geometry.y, units[col])
    for x, y, label in geom_unit:

        # Add label to the axis (xy is location of point)
        text = ax.annotate(
            text=label,
            xy=(x, y),  # location of point
            xytext=(8, 8),  # location of text that goes with point
            textcoords='offset points',  # offset text in points from xy value
            fontsize=fontsize,
            bbox=dict(facecolor='w', alpha=0.3, edgecolor='none',
                      boxstyle="round", pad=0.1))

        rect = text.get_window_extent()
        for other_rect in text_rectangles:
            while Bbox.intersection(rect, other_rect):  # overlapping
                x, y = text.get_position()
                text.set_position((x, (y - y_step)))
                rect = text.get_window_extent()
        text_rectangles.append(rect)

    return (ax)
def plot_map(gdf, col, col_readable, ax, include_legend, legend_kwds,
             units, base_map, show_labels, cmap):
    '''
    Create map of catchment areas, coloured by the provided column.
    Inputs:
    - gdf - geopandas dataframe with IVT catchment areas and col
    - col - string, column to colour map by
    - col_readable - string, column used to colour map
    - ax - axis object, to plot map on
    - include_legend - boolean, whether to display legend
    - legend_kwds - dictionary, characteristics of legend
    - units - dataframe, information on each stroke unit inc. location
    - base_map - boolean, whether to add base map
    - show_labels - boolean, whether to add hospital labels
    - cmap - string, colourmap to use
    '''
    # Plot IVT catchment areas
    gdf.plot(
        ax=ax,
        column=col,  # Column to apply colour
        antialiased=False,  # Avoids artifact boundry lines
        edgecolor='face',  # Make LSOA boundry same colour as area
        linewidth=0.0,  # Use linewidth=0 to hide boarder lines
        # vmin=0,  # Manual scale min (remove to make automatic)
        # vmax=70,  # Manual scale max (remove to make automatic)
        cmap=cmap,  # Colour map to use
        legend=include_legend,  # Whether to display legend
        legend_kwds=legend_kwds,  # Characteristics of legend
        alpha=0.7)  # Set transparency (to reveal basemap)

    # Add title if categorical (for numeric, it is done with legend_kwds)
    if not is_numeric_dtype(gdf[col]) and include_legend:
        leg = ax.get_legend()
        leg.set_title(col_readable)

    # Add country border
    outline.plot(ax=ax, edgecolor='k', facecolor='None', linewidth=1.0)

    # Add location of hospitals and - if True - labels for each hospital
    units.plot(ax=ax, edgecolor='k', facecolor='w', markersize=120, marker='*')
    if show_labels:
        ax = add_nonoverlapping_text_labels(
            units, ax, 'hospital_city', y_step=0.05, fontsize=8)

    # If True, add base map, specifying same CRS and using manual zoom to
    # adjust level of detail
    if base_map:
        ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, zoom=10)

    # Remove x and y ticks
    plt.setp(ax.get_xticklabels(), visible=False)
    plt.setp(ax.get_yticklabels(), visible=False)

    # Turn off axis line and numbers
    ax.set_axis_off()
def create_map(gdf, col, col_readable, gdf_units,
               base_map, show_labels, cmap='inferno_r', include_legend=True):
    '''
    Creates map of England and Wales, with an inset map of London. Shows
    IVT catchment areas, coloured by value from specified column. Uses the
    plot_map() function to create each map (main and inset).
    Inputs:
    - gdf - geopandas dataframe with polygons and col
    - col - string, column to colour map by
    - col_readable - string, used for legend and plot title
    - gdf_units - dataframe, information on each stroke unit inc. location
    - base_map - boolean, whether to add base map
    - cmap - string, colourmap to use on maps, default 'inferno_r'
    - include_legend - boolean, whether to include a legend, default True
    '''
    # Set up figure with max dimensions of 10x10 inch
    fig, ax = plt.subplots(figsize=(20, 20))

    # Set legend kwds if is numeric, else blank
    if is_numeric_dtype(gdf[col]):
        legend_kwds = {'shrink': 0.5, 'label': f'{col_readable}'}
    else:
        legend_kwds = {}

    # Plot map for England and Wales
    plot_map(gdf, col, col_readable, ax,
             include_legend=include_legend,
             legend_kwds=legend_kwds,
             units=gdf_units,
             base_map=base_map,
             show_labels=show_labels,
             cmap=cmap)

    # Make space at bottom for inset map of London
    ax.set_ylim(ax.get_ylim()[0] - 100000, ax.get_ylim()[1])

    # Make space on right for hospital name label
    ax.set_xlim(ax.get_xlim()[0], ax.get_xlim()[1] + 50000)

    # Create inset axis in bottom right corner
    axins = ax.inset_axes([0.8, 0.05, 0.18, 0.18])

    # Identify London hospitals to go in the inset map
    london_mask = gdf_units["Strategic Clinical Network"] == "London SCN"
    london_hospitals = gdf_units["Hospital_name"][london_mask].to_list()
    london_units = gdf_units[london_mask]

    # Identify map area to plot in the inset map
    # initialise exteme values
    minx_ = np.inf
    miny_ = np.inf
    maxx_ = 0
    maxy_ = 0
    # Find min and max x and y for the london hospitals
    for h in london_hospitals:
        minx, miny, maxx, maxy = eng_wales_bounds.loc[h]
        minx_ = min(minx_, minx)
        miny_ = min(miny_, miny)
        maxx_ = max(maxx_, maxx)
        maxy_ = max(maxy_, maxy)

    # Set extent of inset map
    axins.set_xlim(minx_, maxx_)
    axins.set_ylim(miny_, maxy_)

    # Define lines connecting inset map to main map
    mark_inset(ax, axins, loc1=2, loc2=3, fc="none", ec="0.5")

    # Plot inset map of London
    plot_map(gdf, col, col_readable, axins,
             include_legend=False,
             legend_kwds={},
             units=london_units,
             base_map=base_map,
             show_labels=show_labels,
             cmap=cmap)

    # Set title
    ax.set_title(f'{col_readable}')

    # Adjust for printing
    ax.margins(0)
    ax.apply_aspect()
    plt.subplots_adjust(left=0.01, right=1.0, bottom=0.0, top=1.0)
    plt.savefig(f'figures/map_{col}.png', dpi=300)
    plt.show()

    # If colouring column is numeric, do some further analysis
    if is_numeric_dtype(gdf[col]):
        # Print information on the range and median
        display(f'Range: {round(gdf[col].min(),3)} to ' +
                f'{round(gdf[col].max(),3)}')
        display(f'Median: {round(gdf[col].median(),3)}')

        # Histogram showing range in that col
        plt.hist(gdf[col])
        plt.xlabel(col_readable)
        plt.ylabel('Number of catchment areas')
        plt.show()

2 IVT catchment statistics#

IVT catchment map#

Information on using different colours for adjacent polygons available here.

# Use greedy to assign colours wherein adjacent polygons are never the same
# min_distance is the minimal distance between colours
# balance=count means we try to balance number of features per colour
gdf_ivt_catchment['adjacent_colours'] = greedy(
    gdf_ivt_catchment, min_distance=1, balance='count').astype(str)
/home/amy/miniconda3/envs/geopandas/lib/python3.9/site-packages/libpysal/weights/util.py:1658: FutureWarning: The `query_bulk()` method is deprecated and will be removed in GeoPandas 1.0. You can use the `query()` method instead.
  inp, res = gdf.sindex.query_bulk(gdf.geometry, predicate=predicate)
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='adjacent_colours',
    col_readable=('Map of catchment areas for each IVT unit\n(based on ' +
                  'closest unit for each LSOA)'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='viridis',
    include_legend=False)
_images/47808118a9dfa476032ba3bb4598dbe2d6d02b4e3b742e9533805c1bcdb9bfdd.png

Admissions to each unit (2021-22)#

# Get admissions to each unit from 2021/22
admissions_2122 = (gdf_units[['Postcode', 'Admissions 21/22']]
                   .set_index('Postcode')
                   .rename(columns={'Admissions 21/22': 'admissions_2122'})
                   .astype(int))
admissions_2122
admissions_2122
Postcode
RM70AG 981
E11BB 861
SW66SX 1147
SE59RW 824
BR68ND 847
... ...
SY231ER 126
SA148QF 162
SA312AF 207
SA612PZ 208
SA66NL 615

113 rows × 1 columns

# Join with geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(admissions_2122))
# Create map
create_map(gdf=gdf_ivt_catchment.reset_index(),
           col='admissions_2122',
           col_readable='Admissions to each unit (21/22)',
           gdf_units=gdf_units,
           base_map=False,
           show_labels=True,
           cmap='Reds')
_images/74bb80de4b5675e9f9c8aae04472023b5a1627d42bc5769e1ee853682da28865.png
'Range: 126 to 1848'
'Median: 702.0'
_images/89e71d7ff64e91a8d2c78a8d53d2cf44609cc380049af1f007271691fe9fccd0.png

Population-weighted mean income domain score by IVT catchment#

# Weighted mean income domain score by closest IVT unit
df_ivt_ids = weighted_av(
    df_lsoa, group_by='closest_ivt_unit',
    col='income_domain_score', new_col='income_domain_weighted_mean')
df_ivt_ids
mean income_domain_weighted_mean
closest_ivt_unit
B152TH 0.180572 0.182951
B714HJ 0.195096 0.202612
BA13NG 0.083828 0.083469
BA214AT 0.096528 0.095511
BB23HH 0.172242 0.177377
... ... ...
WD180HB 0.079041 0.080071
WF14DG 0.148147 0.148666
WR51DD 0.095798 0.093969
WV100QP 0.168500 0.172614
YO318HE 0.085214 0.084118

113 rows × 2 columns

# Add mean IoD to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    df_ivt_ids['income_domain_weighted_mean']))
# Create map
create_map(gdf=gdf_ivt_catchment.reset_index(),
           col='income_domain_weighted_mean',
           col_readable=('Population weighted mean income deprivation domain ' +
                         'score\n(higher score indicates more deprivation) ' +
                         '(2019)'),
           gdf_units=gdf_units,
           base_map=False,
           show_labels=True,
           cmap='Reds')
_images/ae55c841b25a789b48e9481a8855e5b43734f666cdbba1aba05e237d61541fc4.png
'Range: 0.058 to 0.209'
'Median: 0.122'
_images/12e75c61e2c8fee562ff6fa0f6496b2b2929f94e54e69321bc57113e522f376c.png
4:80: E501 line too long (80 > 79 characters)

Population-weighted mean IMD score by IVT catchment#

# Weighted mean IMD by closest IVT unit
df_ivt_imd = weighted_av(
    df_lsoa, group_by='closest_ivt_unit',
    col='imd_2019_score', new_col='imd_weighted_mean')
df_ivt_imd
mean imd_weighted_mean
closest_ivt_unit
B152TH 31.070594 31.488294
B714HJ 32.588740 33.655519
BA13NG 13.353572 13.335982
BA214AT 16.946994 16.858199
BB23HH 31.070595 31.789134
... ... ...
WD180HB 12.004248 12.153052
WF14DG 26.534618 26.646876
WR51DD 16.780844 16.451837
WV100QP 26.119438 26.753959
YO318HE 14.738411 14.606578

113 rows × 2 columns

# Add mean IMD to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    df_ivt_imd['imd_weighted_mean']))
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='imd_weighted_mean',
    col_readable=('Population weighted mean Index of Multiple Deprivation ' +
                  '(IMD) score\n(higher score indicates more deprived area) ' +
                  '(IMD 2019)'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/bb3645cadf88431d29007f4551451b5df4cc3fdfffa61e0e494e4b647accc5a0.png
'Range: 10.093 to 36.601'
'Median: 20.032'
_images/5627c2daa00773eb02e500e7ef6d245d2e0651bd92070a0dd9f03447837f5a75.png

IVT catchment size#

# Calculate populations
pop = df_lsoa.groupby('closest_ivt_unit')['population_all'].sum()

# Print the range in population size of each catchment area
print(f'The range of people across the {pop.shape[0]} catchments is: ' +
      f'{round(pop.min(),1)} to ' +
      f'{round(pop.max(),1)}')

# Print the median population size
print(f'Median number of people: {round(pop.median(),2)}')

# Create histogram to show distribution in population size of catchment areas
plt.hist(pop)
plt.xlabel('IVT catchment population size')
plt.ylabel('Number of catchment areas')
plt.show()
The range of people across the 113 catchments is: 87086 to 1479138
Median number of people: 483031.0
_images/351f3009afcc2330d9e0bdb44873ab52f70f46dd97bbfe5d386031781b3163fa.png

Ethnicity - proportion ethnic minority by IVT catchment#

Fluck et al. 2023 Nature publication using SSNAP data to compare stroke outcomes in caucasion population versus ethnic minorities, finding people from ethnic minorities have earlier onset and more adverse outcomes.

In the maps, I have present the proportion of people in each IVT catchment area who are an ethnic minorities (language as per government recommendations).

# Show the ethnic groups
[col for col in df_lsoa if col.startswith('ethnic_group')]
['ethnic_group_all_categories_ethnic_group',
 'ethnic_group_white_total',
 'ethnic_group_white_english_welsh_scottish_northern_irish_british',
 'ethnic_group_white_irish',
 'ethnic_group_white_gypsy_or_irish_traveller',
 'ethnic_group_white_other_white',
 'ethnic_group_mixed_multiple_ethnic_group_total',
 'ethnic_group_mixed_multiple_ethnic_group_white_and_black_caribbean',
 'ethnic_group_mixed_multiple_ethnic_group_white_and_black_african',
 'ethnic_group_mixed_multiple_ethnic_group_white_and_asian',
 'ethnic_group_mixed_multiple_ethnic_group_other_mixed',
 'ethnic_group_asian_asian_british_total',
 'ethnic_group_asian_asian_british_indian',
 'ethnic_group_asian_asian_british_pakistani',
 'ethnic_group_asian_asian_british_bangladeshi',
 'ethnic_group_asian_asian_british_chinese',
 'ethnic_group_asian_asian_british_other_asian',
 'ethnic_group_black_african_caribbean_black_british_total',
 'ethnic_group_black_african_caribbean_black_british_african',
 'ethnic_group_black_african_caribbean_black_british_caribbean',
 'ethnic_group_black_african_caribbean_black_british_other_black',
 'ethnic_group_other_ethnic_group_total',
 'ethnic_group_other_ethnic_group_arab',
 'ethnic_group_other_ethnic_group_any_other_ethnic_group']
# Add count of non-white (combines mixed, asian, black + other categories)
df_lsoa['ethnic_group_not_white'] = (
    df_lsoa['ethnic_group_all_categories_ethnic_group'] -
    df_lsoa['ethnic_group_white_total'])
# Find proportion with bad or very bad general health by IVT catchment
ethnicity = find_proportion(
    subgroup_col='ethnic_group_not_white',
    overall_col='ethnic_group_all_categories_ethnic_group',
    proportion_col='ethnic_minority_proportion',
    group_col='closest_ivt_unit',
    df=df_lsoa
)

ethnicity.head()
ethnic_group_not_white ethnic_group_all_categories_ethnic_group ethnic_minority_proportion
closest_ivt_unit
B152TH 296338 906743 0.326816
B714HJ 307488 958616 0.320762
BA13NG 17349 448923 0.038646
BA214AT 5261 261069 0.020152
BB23HH 86542 473321 0.182840
# Add proportion to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    ethnicity['ethnic_minority_proportion']))
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='ethnic_minority_proportion',
    col_readable=('Proportion of people who are an ethnic minority ' +
                  '(2011 census)'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/5b32db56f1717cfa454704a57e7f52b00a9c712d9ec9a1f031adeaf2947e4793.png
'Range: 0.015 to 0.53'
'Median: 0.05'
_images/ad007182bdbd1c0f1d8fec6a05b9baaecf4c18388494a0d529330787886acdf4.png

General health - proportion with bad/very bad health by IVT catchment#

# Find proportion with bad or very bad general health by IVT catchment
bad_health = find_proportion(
    subgroup_col='bad_or_very_bad_health',
    overall_col='all_categories_general_health',
    proportion_col='bad_health_proportion',
    group_col='closest_ivt_unit',
    df=df_lsoa
)

bad_health.head()
bad_or_very_bad_health all_categories_general_health bad_health_proportion
closest_ivt_unit
B152TH 55724 906743 0.061455
B714HJ 65605 958616 0.068437
BA13NG 19467 448923 0.043364
BA214AT 12458 261069 0.047719
BB23HH 32885 473321 0.069477
# Add proportion to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    bad_health['bad_health_proportion']))
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='bad_health_proportion',
    col_readable=('Proportion of people with bad or very bad general health ' +
                  '(2011 census)'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/c6025d0ba2b0ca97ac8c60a2f1594a6ac7f8b4dbc37729736366e57d34d023db.png
'Range: 0.032 to 0.103'
'Median: 0.055'
_images/081f614e2846f6ddc9de5f9e24c09d1049bb5ef85a320a86913cbdebba15248c.png

Long-term health problem or disability - proportion by IVT catchment#

# Add column with count of long-term health problem (limited a little or alot)
df_lsoa['long_term_health_count'] = (
    df_lsoa['all_categories_long_term_health_problem_or_disability'] -
    df_lsoa['day_to_day_activities_not_limited'])
# Find proportion with long-term health problem by IVT catchment
long_term = find_proportion(
    subgroup_col='long_term_health_count',
    overall_col='all_categories_long_term_health_problem_or_disability',
    proportion_col='long_term_health_proportion',
    group_col='closest_ivt_unit',
    df=df_lsoa
)

long_term.head()
long_term_health_count all_categories_long_term_health_problem_or_disability long_term_health_proportion
closest_ivt_unit
B152TH 157057 889227 0.176622
B714HJ 180021 947524 0.189991
BA13NG 70366 438597 0.160434
BA214AT 45407 254179 0.178642
BB23HH 93679 468472 0.199967
# Add proportion to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    long_term['long_term_health_proportion']))
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='long_term_health_proportion',
    col_readable=('Proportion of people with long-term health problem or ' +
                  'disability\n(which limits day-to-day activities a ' +
                  'little or a lot) (2011 census)'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/0dfa7cc8df6f39c42c9046f02c85220a85e14bb6e2ee5f17452a4b8ef1080191.png
'Range: 0.115 to 0.262'
'Median: 0.182'
_images/629b55eac41dbe34f96168a094ddbfc7deeba02950ef7e5fff7496fc9a8b29fb.png

Rural/urban - proportion of people living in rural areas by IVT catchment#

The LSOA rural/urban classification had 8 categories. The local authority rural urban classification is based on the percentage of residents living in rural areas (and then has 6 classes). Hence, for each IVT catchment, I have presented the percentage of residents living in rural areas.

df_lsoa['rural_urban_2011'].value_counts()
rural_urban_2011
Urban city and town                                15724
Urban major conurbation                            11523
Rural town and fringe                               3189
Rural village and dispersed                         2490
Urban minor conurbation                             1208
Rural village and dispersed in a sparse setting      327
Rural town and fringe in a sparse setting            197
Urban city and town in a sparse setting               94
Name: count, dtype: int64
# Extract if is rural or if is urban (each classification starts with that)
df_lsoa['ruc_overall'] = df_lsoa['rural_urban_2011'].str[:5]

# Add column where True if rural
df_lsoa['rural'] = (df_lsoa['ruc_overall'] == 'Rural').map({
    True: 'True', False: 'False'})

# Show columns
df_lsoa[['ruc_overall', 'rural']].value_counts()
ruc_overall  rural
Urban        False    28549
Rural        True      6203
Name: count, dtype: int64
# Find proportion people in rural LSOA for each IVT unit
ruc = find_count_and_proportion(
    df_lsoa=df_lsoa,
    area='closest_ivt_unit',
    col='rural',
    count='population_all',
    prop='proportion_rural'
)

# Preview dataframe
ruc.head()
False True proportion_rural
closest_ivt_unit
B152TH 926409.0 36603.0 0.038009
B714HJ 981575.0 30814.0 0.030437
BA13NG 365427.0 117604.0 0.243471
BA214AT 96824.0 177608.0 0.647184
BB23HH 422949.0 63171.0 0.129949
# Add proportion to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    ruc['proportion_rural']))
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='proportion_rural',
    col_readable=('Proportion of people in rural areas in each IVT ' +
                  'catchment area (based on RUC 2011)'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='viridis')
_images/4f79591daeb94de774d63df7980f086b20960f706d9ae674f250989e33b82c0b.png
'Range: 0.0 to 0.838'
'Median: 0.22'
_images/70fe835b768d662c0a5bf6700dcb89ae284d3b45408e2b2a20d0b8c5d2004a3c.png

Age - proportion aged 65+ by IVT catchment#

The age band reflect the bottom of the band - so age_band_65 is people aged 65, 66, 67, 68 and 69.

# Find count of people age 65 plus
df_lsoa['age_65_plus_count'] = df_lsoa[[
    'age_band_all_65', 'age_band_all_70', 'age_band_all_75',
    'age_band_all_80', 'age_band_all_85', 'age_band_all_90']].sum(axis=1)
# Find proportion age 65 plus by IVT catchment
age_65 = find_proportion(
    subgroup_col='age_65_plus_count',
    overall_col='population_all',
    proportion_col='age_65_plus_proportion',
    group_col='closest_ivt_unit',
    df=df_lsoa
)

age_65.head()
age_65_plus_count population_all age_65_plus_proportion
closest_ivt_unit
B152TH 144555 963012 0.150107
B714HJ 159109 1012389 0.157162
BA13NG 102838 483031 0.212901
BA214AT 72947 274432 0.265811
BB23HH 88140 486120 0.181313
# Add proportion to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    age_65['age_65_plus_proportion']))
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='age_65_plus_proportion',
    col_readable=('Proportion of people aged 65 or older (mid 2020)'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/60c41ea63820d8a21b39d750533391f4f00107d2259661b6ec1e2d70636c7231.png
'Range: 0.082 to 0.287'
'Median: 0.203'
_images/42575e92375bc0ee5f5295e842becb0c5994d1af0031fe8db8fb702230c4b65b.png

Thrombolysis rates for each IVT unit#

Thrombolysis rate here is the rate for the closest IVT unit (not specific to LSOAs - so its proportion of all admissions to that unit, not just of those from that LSOA)

# Extract rate for each unit, with unit as index
ivt_rate = (df_lsoa[['closest_ivt_unit', 'ivt_rate']]
            .drop_duplicates()
            .set_index('closest_ivt_unit')['ivt_rate'])
# Add rate to geopandas dataframe
gdf_ivt_catchment = (gdf_ivt_catchment.join(ivt_rate))
# Create map
create_map(
    gdf=gdf_ivt_catchment,
    col='ivt_rate',
    col_readable='Thrombolysis rate',
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Blues')
_images/ef80b9a59e42be8b61355804956798d84b4da5f62ee1e16b965407f95373f4c1.png
'Range: 1.9 to 27.5'
'Median: 10.6'
_images/9072ae4eec197dd99107ac22c546d2d208c523df80b5c338c1c9ef7fe7c6a9bb.png

Population density#

Calculation of population density based on archive/05_create_proportion_over_65_maps and these stack overflow posts: post 1 and post 2

# Find population of each IVT catchment
pop = df_lsoa.groupby('closest_ivt_unit')['population_all'].sum()
# Add total population size to geopandas
gdf_ivt_catchment = (gdf_ivt_catchment.join(pop))
# Change the projection to a Cartesian system (EPSG:3857) and get area in
# square kilometers (through /10**6)
gdf_ivt_catchment['polygon_area'] = (
    gdf_ivt_catchment.to_crs('epsg:3857').area / 10**6)

# Divide population by area to get population density
gdf_ivt_catchment['population_density'] = (
    gdf_ivt_catchment['population_all'] / gdf_ivt_catchment['polygon_area'])
# Create map
create_map(
    gdf=gdf_ivt_catchment,
    col='population_density',
    col_readable='Population density',
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/7330d8f0e998839262900e86172b195c8cee7e37a06a142e5f7dc3a81120b4c6.png
'Range: 9.482 to 4023.447'
'Median: 152.691'
_images/9594e2a672d158adb4fbedbaf060aa87da2f203afb25b455ed1da7c6647ee2bd.png

Average travel time to IVT unit weighted by number of over 65 year olds in catchment#

# Find mean time to IVT weighted by 65+ population size of each LSOA
ivt_time = weighted_av(
    df=df_lsoa, group_by='closest_ivt_unit',
    col='closest_ivt_unit_time', new_col='weighted_ivt_time',
    weight='age_65_plus_count')

ivt_time
mean weighted_ivt_time
closest_ivt_unit
B152TH 18.971119 19.881289
B714HJ 18.613081 19.355652
BA13NG 24.218246 25.199316
BA214AT 23.251572 24.275645
BB23HH 19.500654 20.295453
... ... ...
WD180HB 19.495918 19.740606
WF14DG 18.663793 18.843356
WR51DD 20.794239 22.107489
WV100QP 18.264433 19.183549
YO318HE 28.439159 30.242922

113 rows × 2 columns

# Add time to geopandas dataframe
# Add total population size to geopandas
gdf_ivt_catchment = (gdf_ivt_catchment.join(ivt_time['weighted_ivt_time']))
# Create map
create_map(
    gdf=gdf_ivt_catchment,
    col='weighted_ivt_time',
    col_readable='Average time to IVT unit (weighted by LSOA 65+ population)',
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/de3ea5771f99f44b33a174667a33b2e87f18360cc26cd2a60044dd4f309b40ee.png
'Range: 12.891 to 35.582'
'Median: 20.089'
_images/20bea7db7311d70939f393a5be33bd619cc2034cbf87016cac86e2db48ac0939.png

Proportion of 65+ year olds living within 30 minutes of their closest IVT unit#

Created as in archive/06_create_ave_travel_time_maps.

# Create column indicating if IVT unit is within 30 minutes
df_lsoa['ivt_within_30'] = (df_lsoa['closest_ivt_unit_time'] < 30).map(
    {True: 'True', False: 'False'})

df_lsoa.ivt_within_30.value_counts()
ivt_within_30
True     30903
False     3849
Name: count, dtype: int64
# Find proportion of over 65 year olds who live within 30 min for each unit
within_30 = find_count_and_proportion(
    df_lsoa=df_lsoa,
    area='closest_ivt_unit',
    col='ivt_within_30',
    count='age_65_plus_count',
    prop='proportion_over_65_within_30'
)

within_30.head()
False True proportion_over_65_within_30
closest_ivt_unit
B152TH 6329.0 138226.0 0.956217
B714HJ 20680.0 138429.0 0.870026
BA13NG 38113.0 64725.0 0.629388
BA214AT 27527.0 45420.0 0.622644
BB23HH 6443.0 81697.0 0.926900
# Add to geopandas
gdf_ivt_catchment = (gdf_ivt_catchment.join(
    within_30['proportion_over_65_within_30']))
# Create map
create_map(
    gdf=gdf_ivt_catchment.reset_index(),
    col='proportion_over_65_within_30',
    col_readable=('Proportion of people aged 65 or older who live within 30 ' +
                  'minutes of their closest IVT unit'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds_r')
_images/fef689f14febb5e16dd26dfe6e0637d69c89f58fbc0d46557f031e806f6fe248.png
'Range: 0.359 to 1.0'
'Median: 0.907'
_images/b338dce91e134acbbcefc458e104725f5ad657e451620b7902d06260a59a0c3a.png

3 Ambulance service statistics#

Ambulance service map#

# Add ambulance service as a column for the mapping function
gdf_amb_catchment['ambulance_service'] = gdf_amb_catchment.index
# Create map
create_map(
    gdf=gdf_amb_catchment,
    col='ambulance_service',
    col_readable='Ambulance service map (based on CCG or county of each LSOA)',
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='viridis')
_images/18a73766c8eeb67f96d9e4871faeee5c4fc7399bd678306ce69218656297b49f.png

Population-weighted mean time to IVT unit by ambulance trust#

This is the population weighted mean time to the closest unit providing thrombolysis (IVT), which may also be referred to as a hyper acute stroke service (HASU). It is population weighted as we have mean time per LSOA, so we weight by the number of individuals in each LSOA.

# Weight mean time to closest IVT unit by ambulance service
df_amb_ivt_time = weighted_av(
    df_lsoa, group_by='ambulance_service',
    col='closest_ivt_unit_time', new_col='ivt_time_weighted_mean')
df_amb_ivt_time.sort_values(by='ivt_time_weighted_mean')
mean ivt_time_weighted_mean
ambulance_service
London 16.379483 16.304175
North East 17.977248 17.842487
Isle of Wight 18.789888 18.605475
North West 18.845056 18.851069
West Midlands 19.284198 19.177783
Yorkshire 19.553911 19.277658
South Central 19.810180 19.808518
Welsh 20.272813 20.130240
East of England 20.745680 20.717194
South East Coast 20.879474 20.826622
East Midlands 21.243307 21.218572
South West 21.566687 21.636569
# Add mean to geopandas dataframe
gdf_amb_catchment = (gdf_amb_catchment.join(
    df_amb_ivt_time['ivt_time_weighted_mean']))
# Create map
create_map(
    gdf=gdf_amb_catchment.reset_index(),
    col='ivt_time_weighted_mean',
    col_readable='Weighted mean time to closest IVT unit by ambulance trust',
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/b34fd825f2bc32ddc419ad0311be29d3bd132e613c6336ecd712fae586542e78.png
'Range: 16.304 to 21.637'
'Median: 19.543'
_images/672e2a49fe4193df4cd2d220ebd1c5459ec1b3007d843778dfc64909dcdd556b.png

Population-weighted mean time to MT unit by ambulance trust#

MT unit might also be referred to as Comprehensive Stroke Centre (CSC).

This considers transfers - i.e. combines (1) time from home to closest IVT unit, and (2) transfer time to its closest MT unit (which, if the same unit, is 0)

# Weight mean time to thrombectomy unit by ambulance service
df_amb_mt_time = weighted_av(
    df_lsoa, group_by='ambulance_service',
    col='total_mt_time', new_col='mt_time_weighted_mean')
df_amb_mt_time.sort_values(by='mt_time_weighted_mean')
mean mt_time_weighted_mean
ambulance_service
London 25.917497 25.450542
West Midlands 37.916748 37.530963
Yorkshire 38.175378 38.060249
North West 39.155098 38.608323
North East 43.506337 43.203993
South Central 47.813957 47.961602
East Midlands 51.961945 52.201445
South East Coast 60.718024 60.631874
South West 63.863551 64.049319
Welsh 64.987114 64.959504
East of England 66.040353 65.884606
Isle of Wight 136.289888 136.105475
# Add mean to geopandas dataframe
gdf_amb_catchment = (gdf_amb_catchment.join(
    df_amb_mt_time['mt_time_weighted_mean']))
# Create map
create_map(
    gdf=gdf_amb_catchment.reset_index(),
    col='mt_time_weighted_mean',
    col_readable='Weighted mean time to thrombectomy unit by ambulance trust',
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/5d2c38856e8d697b468d2d4aff2e2285a9409a7db03f3b86763310934ae980af.png
'Range: 25.451 to 136.105'
'Median: 50.082'
_images/745af19600ce406f263cd2c89ae30abdafb3a4de12a6addcefe6b1eb73913ddf.png

Population-weighted mean transfer time between IVT and MT unit#

# Weight mean transfer time by ambulance trust
df_amb_transfer_time = weighted_av(
    df_lsoa, group_by='ambulance_service',
    col='closest_mt_transfer_time', new_col='mt_transfer_time_weighted_mean')
df_amb_transfer_time.sort_values(by='mt_transfer_time_weighted_mean')
mean mt_transfer_time_weighted_mean
ambulance_service
London 9.538014 9.146366
West Midlands 18.632549 18.353180
Yorkshire 18.621468 18.782592
North West 20.310042 19.757255
North East 25.529089 25.361506
South Central 28.003777 28.153084
East Midlands 30.718638 30.982873
South East Coast 39.838550 39.805252
South West 42.296864 42.412750
Welsh 44.714301 44.829263
East of England 45.294673 45.167412
Isle of Wight 117.500000 117.500000
# Add mean to geopandas dataframe
gdf_amb_catchment = (gdf_amb_catchment.join(
    df_amb_transfer_time['mt_transfer_time_weighted_mean']))
# Create map
create_map(
    gdf=gdf_amb_catchment.reset_index(),
    col='mt_transfer_time_weighted_mean',
    col_readable=('Weighted mean transfer time between IVT and MT unit ' +
                  'by ambulance trust'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds')
_images/9bbae78da1ac18a0e546f1df75cb1693e47f4b136ea94ea6cac92e00d0282f31.png
'Range: 9.146 to 117.5'
'Median: 29.568'
_images/93db6473d7f33b1f9fcc46db41d79b0bc3ef4544730e39ea8d1d03cdb99d8c35.png

Proportion of patients whose nearest hospital is an MT unit#

# Add column with whether closest hospital offers thrombectomy
df_lsoa['closest_hospital_is_mt'] = (
    df_lsoa['closest_mt_transfer_time'] == 0).map({True: 'True',
                                                   False: 'False'})
df_lsoa.head()
LSOA admissions closest_ivt_unit closest_ivt_unit_time closest_mt_unit closest_mt_unit_time closest_mt_transfer closest_mt_transfer_time total_mt_time ivt_rate ... local_authority_district_22 LAD22NM country ethnic_group_not_white long_term_health_count ruc_overall rural age_65_plus_count ivt_within_30 closest_hospital_is_mt
0 Welwyn Hatfield 010F 0.666667 SG14AB 18.7 NW12BU 36.9 CB20QQ 39.1 57.8 6.8 ... Welwyn Hatfield Welwyn Hatfield England 279 194 Urban False 212 True False
1 Welwyn Hatfield 012A 4.000000 SG14AB 19.8 NW12BU 36.9 CB20QQ 39.1 58.9 6.8 ... Welwyn Hatfield Welwyn Hatfield England 388 247 Urban False 203 True False
2 Welwyn Hatfield 002F 2.000000 SG14AB 18.7 NW12BU 38.0 CB20QQ 39.1 57.8 6.8 ... Welwyn Hatfield Welwyn Hatfield England 224 125 Urban False 193 True False
3 Welwyn Hatfield 002E 0.666667 SG14AB 18.7 NW12BU 36.9 CB20QQ 39.1 57.8 6.8 ... Welwyn Hatfield Welwyn Hatfield England 332 129 Urban False 156 True False
4 Welwyn Hatfield 010A 3.333333 SG14AB 18.7 NW12BU 36.9 CB20QQ 39.1 57.8 6.8 ... Welwyn Hatfield Welwyn Hatfield England 200 276 Urban False 256 True False

5 rows × 122 columns

# Find proportion of people whose closest unit offers MT, for each trust
mt_counts = find_count_and_proportion(
    df_lsoa=df_lsoa,
    area='ambulance_service',
    col='closest_hospital_is_mt',
    count='population_all',
    prop='proportion_closest_is_mt'
)

mt_counts.head()
False True proportion_closest_is_mt
ambulance_service
East Midlands 4128292.0 1018729.0 0.197926
East of England 5433023.0 853208.0 0.135726
Isle of Wight 142296.0 0.0 0.000000
London 2861272.0 6141216.0 0.682169
North East 1970544.0 710219.0 0.264932
# Add proportion to geopandas dataframe
gdf_amb_catchment = (gdf_amb_catchment.join(
    mt_counts['proportion_closest_is_mt']))
# Create map
create_map(
    gdf=gdf_amb_catchment.reset_index(),
    col='proportion_closest_is_mt',
    col_readable=('Proportion of people in each ambulance trust whose ' +
                  'stroke unit is one that offers thrombectomy'),
    gdf_units=gdf_units,
    base_map=False,
    show_labels=True,
    cmap='Reds_r')
_images/59f0cc16252df05cf86613b613835fbb7bcd19bbcbb889fc1174d57d287161d1.png
'Range: 0.0 to 0.682'
'Median: 0.241'
_images/016ce8cc2dd797349d31c68c12c6a26f9383c89911ae4363e40914c245973679.png