<jupyter_start><jupyter_text>External validationThis notebook compares classifiction of Great Britain using spatial signatures with three other classifications - WorldPop urban footprints, MODUM and Copernicus Urban Atlas. This comparison serves as an external validation of signatures.<jupyter_code>import geopandas as gpd
import pandas as pd
import seaborn as sns
import scipy as sp
import numpy as np
import xarray
import rioxarray
import rasterstats
import matplotlib.pyplot as plt
from itertools import product
import tobler
import dask_geopandas
import urbangrammar_graphics as ugg


from download import download
from geocube.api.core import make_geocube
def cramers_v(confusion_matrix):
    chi2 = sp.stats.chi2_contingency(confusion_matrix)[0]
    n = confusion_matrix.sum().sum()
    phi2 = chi2/n
    r,k = confusion_matrix.shape
    phi2corr = max(0, phi2-((k-1)*(r-1))/(n-1))
    rcorr = r-((r-1)**2)/(n-1)
    kcorr = k-((k-1)**2)/(n-1)
    return np.sqrt(phi2corr/min((kcorr-1),(rcorr-1)))<jupyter_output><empty_output><jupyter_text>WorldPop urban classificationDownload and read data.<jupyter_code>soton_path = download("https://eprints.soton.ac.uk/446482/1/GB_100m_mbclust6_3x3maj.tif",
                      '../../urbangrammar_samba/spatial_signatures/validation/B_100m_mbclust6_3x3maj.tif')
soton = rioxarray.open_rasterio("../../urbangrammar_samba/spatial_signatures/validation/B_100m_mbclust6_3x3maj.tif")
signatures = gpd.read_parquet("../../urbangrammar_samba/spatial_signatures/signatures/signatures_combined_levels_simplified.pq")<jupyter_output>/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:283: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)<jupyter_text>Reproject.<jupyter_code>soton.rio.crs
soton_osgb = soton.rio.reproject("EPSG:27700")<jupyter_output><empty_output><jupyter_text>Get signature type as float.<jupyter_code>signatures.signature_type = (signatures.kmeans10gb * 10) + signatures.level2
signatures<jupyter_output><empty_output><jupyter_text>Rasterise signatures.<jupyter_code>%%time
signatures_raster = make_geocube(signatures, measurements=["signature_type"], like=soton_osgb)<jupyter_output>/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:283: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)<jupyter_text>Save raster to file.<jupyter_code>signatures_raster.signature_type.rio.to_raster("signatures_raster.tif")<jupyter_output>/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:283: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)<jupyter_text>Create cross tabulation. There may be a smarter way but this one works.<jupyter_code>spsig_vals = signatures.signature_type.unique()
soton_vals = [1, 2, 3, 4, 5, 6, 15]
crosstab = pd.DataFrame(columns=soton_vals, index=spsig_vals)

for sps, so in product(spsig_vals, soton_vals):
    crosstab.loc[sps, so] =  np.logical_and(signatures_raster.signature_type == sps, soton_osgb[0] == so).data.sum()
crosstab
crosstab.columns = [str(c) for c in crosstab.columns]
crosstab.to_parquet("../../urbangrammar_samba/spatial_signatures/validation/ctab_worldpop.pq")
crosstab = pd.read_parquet("../../urbangrammar_samba/spatial_signatures/validation/ctab_worldpop.pq")<jupyter_output><empty_output><jupyter_text>Plot percentages.<jupyter_code>types = {
    0: "Countryside agriculture",
    10: "Accessible suburbia",
    30: "Open sprawl",
    40: "Wild countryside",
    50: "Warehouse/Park land",
    60: "Gridded residential quarters",
    70: "Urban buffer",
    80: "Disconnected suburbia",
    20: "Dense residential neighbourhoods",
    21: "Connected residential neighbourhoods",
    22: "Dense urban neighbourhoods",
    90: "Local urbanity",
    91: "Concentrated urbanity",
    92: "Regional urbanity",
    94: "Metropolitan urbanity",
    95: "Hyper concentrated urbanity",
}
crosstab = crosstab.drop([x for x in crosstab.index if x not in types.keys()])
crosstab.index = [types[x] for x in crosstab.index]
crosstab
order = [
    "Wild countryside",
    "Countryside agriculture",
    "Urban buffer",
    "Open sprawl",
    "Disconnected suburbia",
    "Accessible suburbia",
    "Warehouse/Park land",
    "Gridded residential quarters",
    "Connected residential neighbourhoods",
    "Dense residential neighbourhoods",
    "Dense urban neighbourhoods",
    "Local urbanity",
    "Regional urbanity",
    "Metropolitan urbanity",
    "Concentrated urbanity",
    "Hyper concentrated urbanity",
]
crosstab = crosstab.loc[order]
crosstab.columns.name = "WorldPop class"
crosstab.index.name = "Spatial Signature type"
fig, ax = plt.subplots(figsize=(8, 8))
sns.heatmap(crosstab.divide(crosstab.sum(axis=1), axis=0).astype('float') * 100, cmap=sns.light_palette(ugg.HEX[2], n_colors=256), annot=True, fmt='.0f', vmax=100)
plt.savefig("crosstab_worldpop.pdf", bbox_inches="tight")<jupyter_output><empty_output><jupyter_text>Compute chi-squared<jupyter_code>chi = sp.stats.chi2_contingency(crosstab)
print(f"chi2: {chi[0]}, p: {chi[1]}, dof: {chi[2]}, N: {crosstab.sum().sum()}")<jupyter_output>chi2: 13341832.31867016, p: 0.0, dof: 114, N: 22993921<jupyter_text>Compute Cramer's V.<jupyter_code>cramers_v(crosstab)<jupyter_output><empty_output><jupyter_text>Results indicate moderate associtation.  MODUMLoad MODUM data (manually downloaded from CDRC).<jupyter_code>modum = gpd.read_file("../../urbangrammar_samba/spatial_signatures/validation/modumew2016.zip")
modum<jupyter_output><empty_output><jupyter_text>Load ET celss with signature type label as Dask GeoDataFrame.<jupyter_code>cells = dask_geopandas.read_parquet("../../urbangrammar_samba/spatial_signatures/signatures/signatures_combined_tessellation/")
cells<jupyter_output><empty_output><jupyter_text>Interpolate MODUM classes to ET cells.<jupyter_code>def _join_modum(ch):
    tb = ch.total_bounds
    modum_ch = modum.cx[tb[0]:tb[2], tb[1]:tb[3]]
    return tobler.area_weighted.area_join(modum_ch, ch, ["CLUSTER_CO"])

meta = tobler.area_weighted.area_join(modum, cells._meta, ["CLUSTER_CO"])

modum_joined = cells.map_partitions(_join_modum, meta=meta)<jupyter_output>/opt/conda/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3357: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.

Left CRS: None
Right CRS: EPSG:27700

  if (await self.run_code(code, result,  async_=asy)):<jupyter_text>Select only subset of data needed for crosstab.<jupyter_code>modum_data = modum_joined[["signature_type", "CLUSTER_CO"]]<jupyter_output><empty_output><jupyter_text>Compute Dask-based steps.<jupyter_code>%%time
modum_data = modum_data.compute()
modum_data<jupyter_output>/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:283: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)<jupyter_text>Generate cross tabulation<jupyter_code>crosstab = pd.crosstab(modum_data.signature_type, modum_data.CLUSTER_CO)
crosstab
crosstab.to_parquet("../../urbangrammar_samba/spatial_signatures/validation/ctab_modum.pq")
crosstab = pd.read_parquet("../../urbangrammar_samba/spatial_signatures/validation/ctab_modum.pq")
crosstab.index
types = {
    0: "Countryside agriculture",
    1: "Accessible suburbia",
    3: "Open sprawl",
    4: "Wild countryside",
    5: "Warehouse/Park land",
    6: "Gridded residential quarters",
    7: "Urban buffer",
    8: "Disconnected suburbia",
    20: "Dense residential neighbourhoods",
    21: "Connected residential neighbourhoods",
    22: "Dense urban neighbourhoods",
    90: "Local urbanity",
    91: "Concentrated urbanity",
    92: "Regional urbanity",
    94: "Metropolitan urbanity",
    95: "Hyper concentrated urbanity",
}
crosstab = crosstab.drop([x for x in crosstab.index if int(x) not in types.keys()])
crosstab.index = [types[int(x)] for x in crosstab.index]
crosstab
crosstab = crosstab.loc[order]
crosstab
names = [
    "Suburban Landscapes",
    "Railway Buzz",
    "The Old Town",
    "Victorian Terraces",
    "Waterside Settings",
    "Countryside Sceneries",
    "High Street and Promenades",
    "Central Business District",
]
crosstab.columns = names
crosstab.columns.name = 'MODUM class'
crosstab.index.name = "Spatial Signature type"<jupyter_output><empty_output><jupyter_text>Plot percentages.<jupyter_code>fig, ax = plt.subplots(figsize=(8, 8))
sns.heatmap(crosstab.divide(crosstab.sum(axis=1), axis=0).astype('float') * 100, cmap=sns.light_palette(ugg.HEX[2], n_colors=256), annot=True, fmt='.0f', vmax=100)
plt.savefig("crosstab_modum.pdf", bbox_inches="tight")<jupyter_output><empty_output><jupyter_text>Compute chi-squared<jupyter_code>chi = sp.stats.chi2_contingency(crosstab)
print(f"chi2: {chi[0]}, p: {chi[1]}, dof: {chi[2]}, N: {crosstab.sum().sum()}")<jupyter_output>chi2: 13938867.969895832, p: 0.0, dof: 152, N: 13067584<jupyter_text>Compute Cramer's V.<jupyter_code>cramers_v(crosstab)<jupyter_output>/opt/conda/lib/python3.8/site-packages/ipykernel/ipkernel.py:283: DeprecationWarning: `should_run_async` will not call `transform_cell` automatically in the future. Please pass the result to `transformed_cell` argument and any exception that happen during thetransform in `preprocessing_exc_tuple` in IPython 7.17 and above.
  and should_run_async(code)<jupyter_text>Results indicate moderate associtation.  Urban AtlasLink is no longer valid and needs to be generated at https://land.copernicus.eu/local/urban-atlas/urban-atlas-2018/. The zip contains the full UK classification.<jupyter_code>urban_atlas_path = download("https://land.copernicus.eu/land-files/f0dd642f9db5ca4de46b3cc955de135503446c43.zip",
                      '../../urbangrammar_samba/spatial_signatures/validation/urban_atlas/', kind="zip")
urban_atlas_path = '../../urbangrammar_samba/spatial_signatures/validation/urban_atlas/'
import zipfile
from glob import glob

# zip file handler  
a = zipfile.ZipFile(glob(urban_atlas_path + '*')[0])
a.namelist()
fuas = []

for file in glob(urban_atlas_path + '*'):
    a = zipfile.ZipFile(file)
    gdf = gpd.read_file(file + "!" + a.namelist()[0])
    fuas.append(gdf)
    print(file)
urban_atlas = pd.concat(fuas)
urban_atlas.crs
urban_atlas = urban_atlas.to_crs(27700)
urban_atlas = urban_atlas[urban_atlas.class_2018 != "Other roads and associated land"]
urban_atlas
%%time
ddf_urban_atlas = dask_geopandas.from_geopandas(urban_atlas, npartitions=16)
ddf_urban_atlas.geometry = ddf_urban_atlas.geometry.simplify(5)
simplified_urban_atlas = ddf_urban_atlas.compute()
%%time

def _join_ua(ch):
    tb = ch.total_bounds
    urban_atlas_ch = simplified_urban_atlas.cx[tb[0]:tb[2], tb[1]:tb[3]]
    return tobler.area_weighted.area_join(urban_atlas_ch, ch, ["class_2018"])

meta = tobler.area_weighted.area_join(simplified_urban_atlas, cells._meta, ["class_2018"])

ua_joined = cells.map_partitions(_join_ua, meta=meta)
ua_data = ua_joined[["signature_type", "class_2018"]]
%%time
ua_data = ua_data.compute()
ua_covered = ua_data.dropna()
ua_covered
crosstab = pd.crosstab(ua_covered.signature_type, ua_covered.class_2018)
crosstab
crosstab.to_parquet("../../urbangrammar_samba/spatial_signatures/validation/ctab_ua.pq")
crosstab = pd.read_parquet("../../urbangrammar_samba/spatial_signatures/validation/ctab_ua.pq")
crosstab = crosstab.drop([x for x in crosstab.index if int(x) not in types.keys()])
crosstab.index = [types[int(x)] for x in crosstab.index]
crosstab
crosstab.columns.name = "Urban Atlas class"
crosstab.index.name = "Spatial Signature type"
fig, ax = plt.subplots(figsize=(12, 8))
sns.heatmap(crosstab.divide(crosstab.sum(axis=1), axis=0).astype('float') * 100, cmap=sns.light_palette(ugg.HEX[2], n_colors=256), annot=True, fmt='.0f', vmax=100)
plt.savefig("crosstab_ua.pdf", bbox_inches="tight")<jupyter_output><empty_output><jupyter_text>Compute chi-squared<jupyter_code>chi = sp.stats.chi2_contingency(crosstab)
print(f"chi2: {chi[0]}, p: {chi[1]}, dof: {chi[2]}, N: {crosstab.sum().sum()}")
cramers_v(crosstab)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12448774452949818, 'fraction_numerical': 0.027990412124023815, 'mean_word_length': 6.126170798898071, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 77, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '124427', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5092, 'n_tokens_neox': 4661, 'n_words': 900, 'path': 'clusters/validation.ipynb', 'repo_name': 'urbangrammarai/spatial_signatures', 'size': '12933'}