Measuring the evolution of spatial dependence and spatial inequality: A tutorial using Python

Motivation

Set up your deepnote

Load libraries

# Load libraries import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns import geopandas as gpd import libpysal from libpysal import weights from libpysal.weights.util import min_threshold_distance import esda from esda.moran import Moran, Moran_Local import splot from splot.esda import moran_scatterplot, plot_moran from splot.libpysal import plot_spatial_weights sns.set_style('darkgrid') #sns.set_style('whitegrid')

Load and merge data

df = pd.read_csv("/work/tutorial2021-evolution-of-spatial-dependence-and-inequality/data.csv")

# Load non-spatial data from uploaded file #data = pd.read_csv("data.csv") # OR Load non-spatial data from online repository data = pd.read_csv("https://raw.githubusercontent.com/quarcs-lab/data-open/master/us_income/data.csv") data

gpd.read_file("/work/tutorial2021-evolution-of-spatial-dependence-and-inequality/map.shp")

# Load spatial data from uploaded files #map = gpd.read_file("map.shp") # OR Load spatial data from online repository map = gpd.read_file("https://github.com/quarcs-lab/data-open/raw/master/us_income/map.zip") map

map.crs

help(pd.merge)

Help on function merge in module pandas.core.reshape.merge:

merge(left: 'DataFrame | Series', right: 'DataFrame | Series', how: 'str' = 'inner', on: 'IndexLabel | None' = None, left_on: 'IndexLabel | None' = None, right_on: 'IndexLabel | None' = None, left_index: 'bool' = False, right_index: 'bool' = False, sort: 'bool' = False, suffixes: 'Suffixes' = ('_x', '_y'), copy: 'bool' = True, indicator: 'bool' = False, validate: 'str | None' = None) -> 'DataFrame'
    Merge DataFrame or named Series objects with a database-style join.
    
    A named Series object is treated as a DataFrame with a single named column.
    
    The join is done on columns or indexes. If joining columns on
    columns, the DataFrame indexes *will be ignored*. Otherwise if joining indexes
    on indexes or indexes on a column or columns, the index will be passed on.
    When performing a cross merge, no column specifications to merge on are
    allowed.
    
    Parameters
    ----------
    left : DataFrame
    right : DataFrame or named Series
        Object to merge with.
    how : {'left', 'right', 'outer', 'inner', 'cross'}, default 'inner'
        Type of merge to be performed.
    
        * left: use only keys from left frame, similar to a SQL left outer join;
          preserve key order.
        * right: use only keys from right frame, similar to a SQL right outer join;
          preserve key order.
        * outer: use union of keys from both frames, similar to a SQL full outer
          join; sort keys lexicographically.
        * inner: use intersection of keys from both frames, similar to a SQL inner
          join; preserve the order of the left keys.
        * cross: creates the cartesian product from both frames, preserves the order
          of the left keys.
    
          .. versionadded:: 1.2.0
    
    on : label or list
        Column or index level names to join on. These must be found in both
        DataFrames. If `on` is None and not merging on indexes then this defaults
        to the intersection of the columns in both DataFrames.
    left_on : label or list, or array-like
        Column or index level names to join on in the left DataFrame. Can also
        be an array or list of arrays of the length of the left DataFrame.
        These arrays are treated as if they are columns.
    right_on : label or list, or array-like
        Column or index level names to join on in the right DataFrame. Can also
        be an array or list of arrays of the length of the right DataFrame.
        These arrays are treated as if they are columns.
    left_index : bool, default False
        Use the index from the left DataFrame as the join key(s). If it is a
        MultiIndex, the number of keys in the other DataFrame (either the index
        or a number of columns) must match the number of levels.
    right_index : bool, default False
        Use the index from the right DataFrame as the join key. Same caveats as
        left_index.
    sort : bool, default False
        Sort the join keys lexicographically in the result DataFrame. If False,
        the order of the join keys depends on the join type (how keyword).
    suffixes : list-like, default is ("_x", "_y")
        A length-2 sequence where each element is optionally a string
        indicating the suffix to add to overlapping column names in
        `left` and `right` respectively. Pass a value of `None` instead
        of a string to indicate that the column name from `left` or
        `right` should be left as-is, with no suffix. At least one of the
        values must not be None.
    copy : bool, default True
        If False, avoid copy if possible.
    indicator : bool or str, default False
        If True, adds a column to the output DataFrame called "_merge" with
        information on the source of each row. The column can be given a different
        name by providing a string argument. The column will have a Categorical
        type with the value of "left_only" for observations whose merge key only
        appears in the left DataFrame, "right_only" for observations
        whose merge key only appears in the right DataFrame, and "both"
        if the observation's merge key is found in both DataFrames.
    
    validate : str, optional
        If specified, checks if merge is of specified type.
    
        * "one_to_one" or "1:1": check if merge keys are unique in both
          left and right datasets.
        * "one_to_many" or "1:m": check if merge keys are unique in left
          dataset.
        * "many_to_one" or "m:1": check if merge keys are unique in right
          dataset.
        * "many_to_many" or "m:m": allowed, but does not result in checks.
    
    Returns
    -------
    DataFrame
        A DataFrame of the two merged objects.
    
    See Also
    --------
    merge_ordered : Merge with optional filling/interpolation.
    merge_asof : Merge on nearest keys.
    DataFrame.join : Similar method using indices.
    
    Notes
    -----
    Support for specifying index levels as the `on`, `left_on`, and
    `right_on` parameters was added in version 0.23.0
    Support for merging named Series objects was added in version 0.24.0
    
    Examples
    --------
    >>> df1 = pd.DataFrame({'lkey': ['foo', 'bar', 'baz', 'foo'],
    ...                     'value': [1, 2, 3, 5]})
    >>> df2 = pd.DataFrame({'rkey': ['foo', 'bar', 'baz', 'foo'],
    ...                     'value': [5, 6, 7, 8]})
    >>> df1
        lkey value
    0   foo      1
    1   bar      2
    2   baz      3
    3   foo      5
    >>> df2
        rkey value
    0   foo      5
    1   bar      6
    2   baz      7
    3   foo      8
    
    Merge df1 and df2 on the lkey and rkey columns. The value columns have
    the default suffixes, _x and _y, appended.
    
    >>> df1.merge(df2, left_on='lkey', right_on='rkey')
      lkey  value_x rkey  value_y
    0  foo        1  foo        5
    1  foo        1  foo        8
    2  foo        5  foo        5
    3  foo        5  foo        8
    4  bar        2  bar        6
    5  baz        3  baz        7
    
    Merge DataFrames df1 and df2 with specified left and right suffixes
    appended to any overlapping columns.
    
    >>> df1.merge(df2, left_on='lkey', right_on='rkey',
    ...           suffixes=('_left', '_right'))
      lkey  value_left rkey  value_right
    0  foo           1  foo            5
    1  foo           1  foo            8
    2  foo           5  foo            5
    3  foo           5  foo            8
    4  bar           2  bar            6
    5  baz           3  baz            7
    
    Merge DataFrames df1 and df2, but raise an exception if the DataFrames have
    any overlapping columns.
    
    >>> df1.merge(df2, left_on='lkey', right_on='rkey', suffixes=(False, False))
    Traceback (most recent call last):
    ...
    ValueError: columns overlap but no suffix specified:
        Index(['value'], dtype='object')
    
    >>> df1 = pd.DataFrame({'a': ['foo', 'bar'], 'b': [1, 2]})
    >>> df2 = pd.DataFrame({'a': ['foo', 'baz'], 'c': [3, 4]})
    >>> df1
          a  b
    0   foo  1
    1   bar  2
    >>> df2
          a  c
    0   foo  3
    1   baz  4
    
    >>> df1.merge(df2, how='inner', on='a')
          a  b  c
    0   foo  1  3
    
    >>> df1.merge(df2, how='left', on='a')
          a  b  c
    0   foo  1  3.0
    1   bar  2  NaN
    
    >>> df1 = pd.DataFrame({'left': ['foo', 'bar']})
    >>> df2 = pd.DataFrame({'right': [7, 8]})
    >>> df1
        left
    0   foo
    1   bar
    >>> df2
        right
    0   7
    1   8
    
    >>> df1.merge(df2, how='cross')
       left  right
    0   foo      7
    1   foo      8
    2   bar      7
    3   bar      8

# Merge data and map gdf = map.merge(data, on='id', how='left') gdf

Basic maps (y)

gdf.plot();

gdf.explore()

gdf.plot('ln_y_1994');

gdf.explore('ln_y_1994')

gdf.plot(column='ln_y_1994', cmap='coolwarm', figsize=(15, 10), legend=True);

gdf.explore( column='ln_y_1994', tooltip=['state', 'region', 'ln_y_1994'], scheme='naturalbreaks', k=5, cmap='coolwarm', legend=True, tiles='Stamen Toner' )

gdf.explore( column='ln_y_1994', tooltip=['state', 'region', 'ln_y_1994'], scheme='naturalbreaks', k=5, cmap='coolwarm', legend=True, tiles='CartoDB dark_matter' )

gdf.explore( column='ln_y_1994', tooltip=['state', 'region', 'ln_y_1994'], scheme='naturalbreaks', k=5, cmap='coolwarm', legend=True, tiles='Stamen Terrain', style_kwds =dict(color="gray", weight=0.5), legend_kwds=dict(colorbar=False) )

fig, ax = plt.subplots(figsize=(15,10)) gdf.plot(column="ln_y_1994", scheme='NaturalBreaks', k=5, cmap='coolwarm', legend=True, ax=ax, legend_kwds={'bbox_to_anchor':(0.82, 0.92)}) plt.title('Spatial distribution of log GDP per capita in 1994: Five natural breaks') plt.tight_layout() ax.axis("off") #plt.savefig('myMap.png',dpi=150, bbox_inches='tight') plt.show()

gdf.ln_y_1994.describe()

sns.boxplot(data=gdf, x="ln_y_1994");

fig, ax = plt.subplots(figsize=(15,10)) gdf.plot(column="ln_y_1994", scheme='BoxPlot', cmap='coolwarm', legend=True, ax=ax, legend_kwds={'bbox_to_anchor':(0.82, 0.92)}) plt.title('Spatial distribution of log GDP per capita in 1994: Which regions are below/above the median?') plt.tight_layout() ax.axis("off") #plt.savefig('myMap.png',dpi=150, bbox_inches='tight') plt.show()

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15,10)) gdf.plot(column="ln_y_1929", scheme='BoxPlot', cmap='coolwarm', legend=True, ax=axes[0], legend_kwds={'bbox_to_anchor':(0.99, 0.40)}) gdf.plot(column="ln_y_1994", scheme='BoxPlot', cmap='coolwarm', legend=True, ax=axes[1], legend_kwds={'bbox_to_anchor':(0.99, 0.40)}) plt.tight_layout() axes[0].axis("off") axes[1].axis("off") axes[0].set_title('Log GDP per capita in 1929') axes[1].set_title('Log GDP per capita in 1994') plt.show()

Spatial weights (W)

W = weights.Queen.from_dataframe(gdf) #W = weights.Rook.from_dataframe(gdf) #W = weights.KNN.from_dataframe(gdf, k=4) #W = weights.KNN.from_dataframe(gdf, k=6) #W = weights.KNN.from_dataframe(gdf, k=8) # Compute inverse distance squared based on maximum nearest neighbor distance between the n observations #maxMin = weights.min_threshold_dist_from_shapefile("map.shp") #W = weights.DistanceBand.from_dataframe(gdf, threshold = maxMin, binary=False, alpha=-2) # Row-standardize W W.transform = 'r'

plot_spatial_weights(W, gdf);

/opt/conda/lib/python3.9/site-packages/splot/_viz_libpysal_mpl.py:115: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  centroids_shp = gdf.centroid.values
/opt/conda/lib/python3.9/site-packages/splot/_viz_libpysal_mpl.py:154: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  gdf.centroid.plot(ax=ax, **node_kws)

W.neighbors

W.min_neighbors

W.mean_neighbors

W.max_neighbors

W.histogram #Note: (numbers of neighbors, frequency)

Spatial lag (Wy)

# Create spatial lag variables gdf["Wln_y_1929"] = weights.lag_spatial(W, gdf["ln_y_1929"])

gdf[["ln_y_1929", "Wln_y_1929"]].round(2)

f,ax = plt.subplots(1,2, figsize=(16, 8)) gdf.plot(column='ln_y_1929', ax=ax[0], edgecolor='k', scheme='NaturalBreaks', k=5, cmap='coolwarm', legend=True, legend_kwds={'bbox_to_anchor':(1.05, 0.40)}) ax[0].set_title("Log GDP pc in 1929") ax[0].axis("off") gdf.plot(column='Wln_y_1929', ax=ax[1], edgecolor='k', scheme='NaturalBreaks', k=5, cmap='coolwarm', legend=True, legend_kwds={'bbox_to_anchor':(1.05, 0.40)}) ax[1].set_title("Spatial lag of log GDP pc in 1929") ax[1].axis("off");

Spatial dependence (Wy, y)

g = sns.jointplot(data=gdf, x="ln_y_1929", y="Wln_y_1929", kind="reg") g.set_axis_labels('Log GDPpc in 1929', 'Log GDPpc in neighbouring states in 1929',);

moran1929 = Moran(gdf["ln_y_1929"], W) moran1994 = Moran(gdf["ln_y_1994"], W)

print(moran1929.I, moran1994.I)

0.6529030905375451 0.3570856202115702

print(moran1929.p_sim, moran1994.p_sim)

0.001 0.001

moran_scatterplot(moran1929);

moran_scatterplot(moran1994);

Evolution of spatial dependence

# Create multidimentional array dfa = gdf.loc[:,'ln_y_1929':'ln_y_1994'].values dfa

dfa.shape

# Create array of years (Be careful with the index. Add +1) years = np.arange(1929,1995) years

mits = [Moran(cs, W) for cs in dfa.T] res = np.array([(m.I, m.EI, m.p_sim, m.z_sim) for m in mits])

plt.plot(years, res[:,0]) plt.title("Evolution of spatial dependence") plt.ylabel("Moran's I of log GDP pc");

Evolution of spatial inequality

for row in dfa: plt.plot(years, row) plt.ylabel('Log GDP per capita');

rel_dfa = dfa / dfa.mean(axis=0) for row in rel_dfa: plt.plot(years, row) plt.ylabel('Relative log GDP per capita');

sigma = dfa.std(axis=0) plt.plot(years, sigma) plt.title("Evolution of spatial inequality") plt.ylabel('Stand. Dev.of log GDP pc');

Summary graphs

plt.plot(years, res[:,0], label="Sp. Dependence") plt.plot(years, sigma, label="Sp. Inequality") plt.legend(loc="upper right") plt.savefig('evolution.png',dpi=150, bbox_inches='tight');

# create figure and axis objects with subplots() fig,ax = plt.subplots() # make a plot ax.plot(years, res[:,0], color="royalblue") # set y-axis label ax.set_ylabel("Spatial Dependence", color="royalblue") # twin object for two different y-axis on the sample plot ax2=ax.twinx() # make a plot with different y-axis using second axis object ax2.plot(years, sigma, color="orange") ax2.set_ylabel("Spatial Inequality", color="darkorange") sns.set_style("dark") plt.show()