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")
df
# 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()