ENVS 563 Computational Essay¶

# Import required libraries

# Mathematical functions and matrix manipulation
import numpy as np

# General data handling
import pandas as pd
import geopandas as gpd

# Creating plots
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
# Adding basemaps to choropleths
import contextily as cx

# Creating spatial weights matrices
from pysal.lib import weights
# Calculating moran's I
import esda

# k means clustering for geodemographic classification
from sklearn import cluster
# davies bouldin score for clustering evaluation
from sklearn.metrics import davies_bouldin_score

# Read dataset from course website as a geodataframe to make use of geometry
url = 'https://darribas.org/gds_course/_downloads/a2bdb4c2a088e602c3bd6490ab1d26fa/toronto_socio-economic.gpkg'
neis = gpd.read_file(url)

# Create a new area column using local projected coordinate reference system
neis['areaGeometry_km2'] = neis.to_crs('EPSG:4536').area / 1000**2
# Drop original area column due to discrepancy with calculated area
neis = neis.drop(columns = ['Shape__Area'])

# Create list of all columns holding age variables to simplify handling
# As pop85+ column holds only relatively small counts and pop65+ column already includes pop85+ counts, pop85+ will not be used
ageCols = ['pop_0-14_yearsold', 'pop_15-24_yearsold', 'pop_25-54_yearsold', 'pop_55-64_yearsold', 'pop_65+_yearsold']

# Create a new total population count column
neis['popTotal'] = neis[ageCols].sum(axis=1)
# Drop original population column due to count discrepancy
neis = neis.drop(columns = ['population2016'])

# Create new aggregated groups as percentages of the total population of each neighbourhood
neis['youngPeopleRate%'] = 100 * (neis['pop_0-14_yearsold'] + neis['pop_15-24_yearsold']) / neis['popTotal']
neis['workingAgeRate%'] = 100 * (neis['pop_25-54_yearsold'] + neis['pop_55-64_yearsold']) / neis['popTotal']
neis['retirementAgeRate%'] = 100 * (neis['pop_65+_yearsold']) / neis['popTotal']

# Create population density variable from calculated total population and geometry area
neis['popDensity_km-2'] = neis['popTotal']/neis['areaGeometry_km2']
# Drop original population density column
neis = neis.drop(columns = ['population_sqkm'])
# Create natural log of population density variable due to right skewed population density distribution
neis['natLog(PopDensity_km-2)'] = np.log(neis['popDensity_km-2'])

# Create list of all columns holding bedroom count variables to simplify handling
bedroomCols = ['bedrooms_0', 'bedrooms_1', 'bedrooms_2', 'bedrooms_3', 'bedrooms_4+'] 
# Create count of properties from sum of all bedroom counts
neis['properties'] = neis[bedroomCols].sum(axis=1)

# Create mean bedroom count for each neighbourhood as weighted average
neis['meanBedroomCount'] = (0*neis['bedrooms_0'] + 1*neis['bedrooms_1'] +\
                            2*neis['bedrooms_2'] + 3*neis['bedrooms_3'] +\
                            4*neis['bedrooms_4+'])/neis['properties']

# Create employment rate using only working age population
neis['employmentRate%'] = 100 * neis['employed'] / (neis['pop_25-54_yearsold'] + neis['pop_55-64_yearsold'])

# Read dataset as a geodataframe to enable spatial manipulation
parks = gpd.read_file('envs563/parks/CITY_GREEN_SPACE_WGS84.shp')
# Remove NaN geometry row
parks = parks.iloc[parks['geometry'].dropna().index]

# Read dataset as a geodataframe to enable spatial manipulation
trees = gpd.read_file('envs563/TopographicMapping–TreedArea/TOPO_TREED_AREA_WGS84.shp')

# Read in data
walkscore = pd.read_csv('envs563/walkscore/walkscore.csv')

# A number of the neighbourhoods have slightly different names in the walkscore dataset
# Create a dictionary to map walkscore df names to neis gdf names
walkscoreMap = {'Danforth Village - East York':'Danforth East York',
              'Cabbagetown-South St.James Town':'Cabbagetown-South St. James Town',
              'Dovercourt-Wallace Emerson-Juncti':'Dovercourt-Wallace Emerson-Junction',
              'Oakwood-Vaughan':'Oakwood Village',
              'Crescent Town':'Taylor-Massey',
              'Markland Woods':'Markland Wood',
              'Danforth Village - Toronto':'Danforth',
              'Mimico':'Mimico (includes Humber Bay Shores)',
              'North St.James Town':'North St. James Town',
              'Weston-Pellam Park':'Weston-Pelham Park',
              'Corsa Italia-Davenport':'Corso Italia-Davenport',
              'Caledonia-Fairbanks':'Caledonia-Fairbank'
             }

# Replace walkscore names with neis names using mapping
walkscore = walkscore.replace(walkscoreMap)

# Merge walkscore ratings to main neighbourhood dataset using neighbourhood names
neis = pd.merge(neis, walkscore[['Name', 'Walk Score', 'Transit Score', 'Bike Score']],
                left_on='neighbourhood_name', 
                right_on='Name', 
                copy=False)

# Drop walkscore neighbourhood name column due to repetition
neis = neis.drop(columns=['Name'])

# Aggregate walk score and bike score variables into a single active transport score to reduce variable count
neis['activeTransportScore'] = (neis['Walk Score'] + neis['Bike Score']) / 2

# Create spatial weights matrix from neighbourhood geometry based on queen contiguity
wQueen = weights.Queen.from_dataframe(neis)
# Row standardise spatial weight matrix
wQueen.transform = 'R'

def globalSpatialAutocorrelationPlotter(gdf, cols, swm):
    """
    Creates a figure containing 1 row of subplots
    Each subplot displays the moran plot for a variable in cols
    
    Parameters
    ----------
    gdf : geopandas geodataframe
        gdf containing columns of interest for gsa plots
        
    cols : list[string]
        list containing column names from gdf to plot
    
    swm : pysal spatial weights matrix
        spatial weights matrix based on gdf geometry
    """
    
    # Initialise empty geodataframe to store calculated variables
    gdfMoranPlot = gpd.GeoDataFrame()
    
    # Initialise empty dictionary to hold calculated moran's I for each variable
    miDict = {}
    
    # Create lists of new column names: standardised columns and spatial lag of standardised columns
    stdCols = [col + '_Standardised' for col in cols]
    slStdCols = [col + '_SpatialLag' for col in stdCols]       
    
    # Iterate over cols list
    for i, col in enumerate(cols):
        # Standardise col and add to moran plot gdf
        gdfMoranPlot[stdCols[i]] = (gdf[col]-gdf[col].mean())/gdf[col].std()
        # Calculate spatial lag of standardised col and add to moran plot gdf
        gdfMoranPlot[slStdCols[i]] = weights.lag_spatial(swm, gdfMoranPlot[stdCols[i]])
        # Calculate moran's I for standardised col and add to moran's I dictionary
        miDict[i] = str(esda.Moran(gdfMoranPlot[stdCols[i]], swm).I.round(3))
    
    # Figure to display moran plots
    # Share both axes across subplots to aid comparison
    fig, axs = plt.subplots(1, len(cols), figsize=(18, 6), sharex=True, sharey=True) 
    
    # Iterate over subplots array
    for i, ax in enumerate(axs):
        # Add moran plot from calculated columns stored in moran plot gdf
        sns.regplot(ax=ax, x=gdfMoranPlot[stdCols[i]], y=gdfMoranPlot[slStdCols[i]], ci=None)        
        # Add moran's I statistic in lower right quadrant
        ax.text(2, -1.5, 'Moran\'s I:\n' + miDict[i], fontsize=14)
        # Add lines at 0 on both axes to indicate mean
        ax.axvline(0, c='k', alpha=0.5)
        ax.axhline(0, c='k', alpha=0.5)
        # Add subplot title
        ax.set_title(cols[i])
        
    # Add figure title
    fig.suptitle('Global spatial autocorrelation of selected variables', fontsize=18)    
    return

def kMeansGeodemographicClassification(gdf, cols, kRange, kMap):
    """
    Creates a geodemographic classification using k means clustering algorithm
    Displays results of classification based on kMap clusters
    Plots 2 additional subplots showing inertia and davies-bouldin scores for values of k in kRange
    
    Parameters
    ----------
    gdf : geopandas geodataframe
        gdf containing columns of interest for k means clustering and geometry
        
    cols : list[string]
        list containing names of columns to be used for k means clustering
    
    kRange : range
        range of k values for which to evaluate inertia and davies-bouldin scores
        
    kMap : integer
        value of k used to create map
    
    Returns
    -------
    gdfClstr['kcls'] : pandas series
        series holding label for each row in gdf resulting from clustering for kMap clusters
    """
    
    # Initialise empty lists to store inertia and davies bouldin scores for each k
    inertias = []    
    dbs = []    
    # Create list of standardised column names
    stdCols = [col + 'Std' for col in cols]    
    
    # Create geodataframe to store standardised columns for clustering
    # Add original gdf geometry for plotting
    gdfClstr = gpd.GeoDataFrame({'geometry':gdf['geometry']})
    
    # Add standardised of columns of interest to clustering geodataframe
    for i, col in enumerate(cols):
        gdfClstr[stdCols[i]] = (gdf[col]-gdf[col].mean())/gdf[col].std()
    
    # Iterate over each value of k in range
    for k in kRange:                
        # Create k means clustering object for current k, random state for reproducibility
        kmeans = cluster.KMeans(n_clusters=k, random_state=321)
        # Identify cluster members using standardised columns of interest
        kcls = kmeans.fit(gdfClstr[stdCols])                
        # Add cluster labels to clustering geodataframe
        gdfClstr['kcls'] = kcls.labels_

        # Add inertia score to inertia score list directly
        inertias.append(kcls.inertia_)
        # Calculate davies bouldin score and add to list
        dbs.append(davies_bouldin_score(gdfClstr[stdCols], gdfClstr['kcls']))
        
        # If current k matches kMap, store current labels
        if k == kMap:
            kMapLabels = kcls.labels_.copy()

    # Overwrite clustering labels from final iteration with labels associated with kMap value of k, for plotting map
    gdfClstr['kcls'] = kMapLabels
    
    # Figure to display k means clustering results
    fig = plt.figure(figsize=(18,10))
    # Use gridspec to manage subplot positions
    gs = gridspec.GridSpec(3, 4)
            
    # Map geodemographic classification based on kMap value of k
    # Position upper centre
    ax1 = fig.add_subplot(gs[:2, :])
    
    gdfClstr.plot(ax=ax1,
                  column='kcls',
                  legend=True,
                  categorical=True,
                  cmap = 'Set1',
                  alpha=0.9,
                  legend_kwds={'loc':'lower right'}
                 )
    # Add title and remove axes
    ax1.set_title(f"Classification based on {kMap} clusters")
    ax1.set_axis_off()
    # Add basemap for context
    cx.add_basemap(ax1, 
                   crs=gdfClstr.crs,
                   source = cx.providers.Stamen.TonerBackground)
    
    # Plot inertia values for each k in kRange
    # Position lower left
    ax2 = fig.add_subplot(gs[2, :2])
    ax2.plot(kRange, inertias, '-o', color='b')
    # Add subplot title and x axis label
    ax2.set_title('Inertia')
    ax2.set_xlabel('k')
    # Add lines at kMap value
    ax2.axvline(kMap, c='k', alpha=0.5)
    ax2.axhline(inertias[kRange.index(kMap)], c='k', alpha=0.5)
    
    # Plot davies bouldin scores for each k in kRange
    # Position lower right
    ax3 = fig.add_subplot(gs[2, 2:])
    ax3.plot(kRange, dbs, '-o', color='b')
    # Add subplot title and x axis label
    ax3.set_title('Davies-bouldin score')
    ax3.set_xlabel('k')
    # Add lines at kMap value
    ax3.axvline(kMap, c='k', alpha=0.5)
    ax3.axhline(dbs[kRange.index(kMap)], c='k', alpha=0.5)
    
    # Add figure title
    fig.suptitle('K means clustering based geodemographic classification of Toronto neighbourhoods', fontsize=18)        
    
    # Return cluster labels from kMap clusters
    return gdfClstr['kcls']