<jupyter_start><jupyter_text>CO2 emissions cluster analysisThe goal of this analysis is to use [BEIS measures of the CO2 emissions](https://data.gov.uk/dataset/723c243d-2f1a-4d27-8b61-cdb93e5b10ff/uk-local-authority-and-regional-carbon-dioxide-emissions-national-statistics-2005-to-2019) under a local authorities control (the subset dataset) to:1. Determine for each council which other councils are most similar. 2. Group councils into clusters that explain broadly the problem of their emissions and why they are similar. This first attempt uses as dimensions the five grand totals for Industry, Commerical, Domestic, Public Sector, Transport.  The data is conformed to the current 2021 local authorities, and reduced to just principal authorities.Although the data has some higher level authorities, their totals are then recreated using the sum of their constituent parts. All authorities are then divided by the 2020 population to get a per capita figure. Go back and forth on if this is a good idea, but essentally lets us isolate the question of 'big' authorities for a seperate comparison system.One thing the per capita approach causes problems with is the City of London, which shows up as such an outlier I exclude it from clustering approaches.<jupyter_code>data_year = 2019

df = pd.read_excel(Path("data", "2005-19_UK_local_and_regional_CO2_emissions.xlsx"),
                   sheet_name="Subset dataset", header=1)
df = df.loc[df["Year"] == data_year].loc[~df["Code"].isna()]

used_cols = ["Industry Total", "Commercial Total", "Domestic Total", "Public Sector Total", "Transport Total"]

df = df[["Code"] + used_cols]

# transform to principal modern councils
df = (df
      .la.create_code_column(from_type="gss", source_col="Code")
      .drop(columns=["Code"])
      .la.to_current()
      .la.just_principals())

higher_df = df.la.to_higher(aggfunc="sum")

df = pd.concat([df, higher_df])

original_df = df.copy()

df = df.la.get_council_info(["region", "pop_2020"])

#convert to per capita emissions
for c in used_cols:
      df[c] = (df[c]/df["pop_2020"])*100

df = df.drop(columns=["pop_2020"])
df["region"] = df["region"].fillna("Unknown")
df = df.set_index("local-authority-code")
df.head()<jupyter_output><empty_output><jupyter_text>Comparsions to physical distanceComparing the distance between all councils (where the input variables are normalised so there is a mean of 0 and a SD of 1) to the physical distance between councils there is very little agreement on the closest councils. The table below shows the agreement between councils in a top k comparison (where k = 1 is the nearest council, k = 10 is the nearest 10 councils). There is on average only a 11% agreement in the top 25 councils.<jupyter_code>relevant_cols = ["Industry Total", "Commercial Total", "Domestic Total", "Public Sector Total", "Transport Total"]
pol_distance = df.space.self_distance("local-authority-code", relevant_cols, normalize=True)

pol_distance.to_csv(Path("data","outputs","council_clusters_v1_distance.csv"), index=False)

real_distance = pd.read_csv(Path("data", "outputs", "geographic_distance.csv" ))

jd = space.join_distance({"pollution_distance": pol_distance,
                         "real_distance": real_distance})
jd.joint_space.agreement()

plot = jd.joint_space.plot(sample=1000, title="1000 samples of distances between authorities")<jupyter_output><empty_output><jupyter_text>There is Similar dispersal of distances, but no clear pattern (although not many far away, and *also* similar pollution).Remember all the coordinates were normalised before use, so this is *distance* between spaces on normalised dimensions. This also shows that some councils can have quite a high pollution distance, but be very *close*. This validates the general approach of exploring alternatives to distance.   ClusteringUsing k means clustering, this section tries to group local authorities into a set of clusters of different profiles of emissions. Because this has been adjusted per capita, the City of London gets very high numbers, and so fits badly into other clusters. It is dropped from these purposes (or can be said to be a cluster of 1).<jupyter_code># drop london because it ends up on its own anyway and distorts the scores
cdf = df.drop(index=["LND"]).reset_index()
cluster = space.Cluster(cdf, "local-authority-code", cols = relevant_cols, label_cols=["region"], normalize=True)
cluster.find_k(2,15)<jupyter_output><empty_output><jupyter_text>Chosing 6 clusters The relatively slow decent reflects that there is not clear clusters in the data, cloosing 6 as a high point in the silhouette score.<jupyter_code># Choosing 6 clusters, but in general this isn't that much different. Data is not clearly clustered.
cluster.set_k(6)

nice_name_lookup = la.get_la_df().set_index("local-authority-code")["official-name"].to_dict()
stats = cluster.stats(label_lookup=nice_name_lookup)
stats<jupyter_output><empty_output><jupyter_text>Broadly, this splits the data into three 'Mainstream' categories, and three 'Outlier' categories. Where the majority are in the mainstream categories, and the outliers are generally better defined (but much smaller) clustesr. Using the label and raincloud tool, descriptions were created of how these labels relate to the variables being grouped.<jupyter_code># examine variables individually
cluster.raincloud_tool()
cluster.label_tool()
cluster.assign_name(n=1, name="Mainstream B", desc="Above average industry/transport/domestic")
cluster.assign_name(n=2, name="Mainstream A", desc="Below average commercial/industry/transport/domestic")
cluster.assign_name(n=3, name="Outlier A", desc="Above average transport, industry, domestic")
cluster.assign_name(n=4, name="Outlier B", desc="Above average commercial/public sector")
cluster.assign_name(n=5, name="Outlier C", desc="Very above average domestic and transport")
cluster.assign_name(n=6, name="Mainstream C", desc="Above average commercial/public sector, below average transport/domestic")
cluster.stats(label_lookup=nice_name_lookup)<jupyter_output><empty_output><jupyter_text>Area and populationBreaking down the labels by area and population, While Mainstream A is the largest group by population, Mainstream B is more rural and contains the largest amount of the land. Similarly, Outlier A has a very low density.<jupyter_code>df = cluster.df_with_labels()
df = df.la.get_council_info(["area", "pop_2020"])
df = df.la.just_principals()
pt = df.pivot_table(["area", "pop_2020"], index=["label"], aggfunc="sum")
pt["area"] = (pt["area"] / pt["area"].sum()).round(3)
pt["pop_2020"] = (pt["pop_2020"] / pt["pop_2020"].sum()).round(3)
pt = pt.reset_index()
pt.style.format({
    'area': '{:.1%}'.format,
    'pop_2020': '{:.1%}'.format
})<jupyter_output><empty_output><jupyter_text>Relabelling clustersReviewing the map below gave a clearer sense of the urban/rural divide in clusters (and makes the (in retrospect, obvious) case for including a measure based on RUC in future cluster label comparisons. Based on these the stand-in categories were renamed with more human descriptions.<jupyter_code>cluster.assign_name(n=1, name="Rural - Mainstream", desc="Above average industry/transport/domestic. Low density.")
cluster.assign_name(n=2, name="Urban - Mainstream", desc="Below average commercial/industry/transport/domestic. High density.")
cluster.assign_name(n=3, name="Rural - Industrial", desc="Above average transport, industry, domestic. Low density.")
cluster.assign_name(n=4, name="Urban - High Commercial", desc="Very high commercial/public sector.")
cluster.assign_name(n=5, name="High Domestic Counties", desc="Very above average domestic and transport in county councils.")
cluster.assign_name(n=6, name="Urban - Commercial", desc="Above average commercial/public sector, below average domestic/transport. High density.")
cluster.stats(label_lookup=nice_name_lookup)
#export with pre per-capita numbers

ldf = cluster.df_with_labels().set_index("local-authority-code")
ldf = ldf[["label", "label_desc"]]
ndf = original_df.copy().set_index("local-authority-code")
ndf = ndf.join(ldf, how="left").reset_index()
ndf = ndf.la.get_council_info(["official-name", "local-authority-type-name"])
ndf.to_csv(Path("data", "outputs", "council_clusters_v1.csv"), index=False)<jupyter_output><empty_output><jupyter_text>Is there any benefit to increasing the number of variables?The `pollution.ks` notebook explores using *all* the sub-variables. Broadly this produces a set of coherent clusters, but still has the same gradual decline and puts most councils in a single mainstream category. The approach in this notebook is preferable for the moment.   How does this work in practice?In terms of how this information is presented, there are several considerations:* All variables have now been transformed several times, so the raw 'distance' means very little. * Per capita emissions is very abstract anyway. So in the example below, there are no units and everything is presented in relative terms:* Per capita emissions are grouped into deciles (Where 10 is high).* Distance is reframed as a percentage match score, where a perfect match is 0 distance, and a 0% match is the distance of the furthest but 10 councils. This is because per capita makes London so far away that it distorts all other distances. So going 10 from the end is stopping this outlier distortion.The below shows an example for Leeds (an interactive version can be seen using the `lookup_demo.ipynb`).<jupyter_code>dmd = lambda x: display(md(x))

# prepare decile guide
cdf = pd.read_csv(Path("data", "outputs", "council_clusters_v1.csv"))
cdf = cdf.la.get_council_info(["pop_2020"]).set_index(
    "local-authority-code")

for c in cdf.columns:
    if "Total" in c:
        name = " ".join(c.split(" ")[:-1]) +" Emissions Decile"
        per_capita = cdf[c] / cdf["pop_2020"]
        cdf[name] = pd.qcut(per_capita, 10,
                                    labels=range(1, 11), duplicates="drop")

ddf = pd.read_csv(Path("data", "outputs", "council_clusters_v1_distance.csv")) #"council_clusters_v1_distance.csv"


def standardise_distance(df):
    df = df.copy()
    # use tenth from last because the last point might be an extreme outlier (in this case london)
    tenth_from_last_score = df["distance"].sort_values().tail(10).iloc[0]
    df["match"] = 1 - (df["distance"] / tenth_from_last_score)
    df["match"] = df["match"].round(3) * 100
    df = df.sort_values("match", ascending=False)
    df = df.drop(columns="local-authority-code_A")
    return df

ddf = ddf.groupby("local-authority-code_A").apply(standardise_distance)

def display_council(council_slug):
    wanted_columns = ["official-name", "match",
                    "local-authority-type-name", "label"]

    matches = ddf.loc[council_slug].sort_values("match", ascending=False).head(10)
    matches = matches.set_index("local-authority-code_B").join(cdf).reset_index()
    matches = matches[wanted_columns]

    details = cdf.loc[council_slug].to_dict()

    dmd(f'## {details["official-name"]}')
    dmd(f'**Type**: {details["local-authority-type-name"]}')
    dmd(f'**Cluster**: {details["label"]}')
    dmd(f'**Cluster description**: {details["label_desc"]}')

    deciles = cdf.loc[cdf.index == council_slug][[x for x in cdf.columns if "Decile" in x]]
    deciles = deciles.transpose().rename(columns={council_slug:"Decile"})
    dmd("### Decile profiles")
    dmd("Decile 10 is the highest emissions, decile 1 is the lowest.")

    display(deciles)

    dmd("### Nearest neighbours")
    display(matches)

display_council("LDS")<jupyter_output><empty_output><jupyter_text>Map of clustersAnd finally, a map for lower-tier authorities of their respective clusters.<jupyter_code>df = ndf[["local-authority-code", "label"]]
# remove london and label for HDCs
df = df.loc[df["local-authority-code"] != "LND"]
df = df.loc[df["label"] != "High Domestic Counties"]
df["label"] = df["label"].fillna("Other")
subtitle = ""
la_layer = lower_tier_layer().merge(df, how="left")


palette = ["colour_dark_blue",
           "colour_berry",
           "colour_orange",
           "colour_blue",
           "colour_purple",
           "colour_yellow"]

#labels = list(df["label"].unique())
#labels.sort(key=lambda x: x if x != "Other" else "Z", reverse=True)
labels = ["Urban - Mainstream",
          "Urban - Commercial",
          "Urban - High Commercial",
          "Rural - Mainstream",
          "Rural - Industrial"]

colors = ["colour_dark_blue", "colour_blue", "colour_orange", "colour_berry", "colour_yellow"]

color_scale = altair_theme.color_scale(domain=labels, named_palette=colors )

# create chart layer
chart = gdf_chart(la_layer).mark_geoshape(
    stroke='black',
    strokeWidth=0
).encode(color=alt.Color("properties.label:O",
                         scale=color_scale,
                         legend=alt.Legend(title=subtitle, orient="top", titleFontSize=18, labelFontSize=15, titleLimit=400))
         ).properties(
    title="Lower tier authorities clustered by CO2 emissions profiles (v1)"
)

# put between layers
chart = lower_tier_sandwich(chart)

# final properties
chart = chart.properties(width=800, height=600)
caption = "Source: mySociety analysis of BEIS dataset; House of Commons LA cartogram (Open Parliament Licence)"
chart.display(logo=True, scale_factor=1, caption=caption)
save_chart(chart, Path("data", "outputs", "cluster_map.png"),
           scale_factor=2, logo=True, caption=caption)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11584869346361684, 'fraction_numerical': 0.014138722333259308, 'mean_word_length': 5.205787781350482, 'pattern_counts': {'":': 2, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '284512', 'license': 'MIT', 'n_tokens_mistral': 4246, 'n_tokens_neox': 3951, 'n_words': 1293, 'path': 'notebooks/old_approach/emissions.v1.ipynb', 'repo_name': 'mysociety/la_distance', 'size': '13509'}