SIEW (2021) On the spatial concentration of economic activity in East Asia and Pacific: New evidence from satellite night-time lights
This is a computational notebook to replicate the results of the Global Moran scatterplot, Local Moran map and Directional LISAs for the above-mentioned paper, to execute, use this link: https://deepnote.com/project/Global-and-Local-Moran-Dynamic-LISAs-UcR4wDs6Qt-a78r_YNpx9g/%2Fsatellite.ipynb
Load required libraries
import geopandas as gpd
import contextily as cx
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from pysal.viz import mapclassify
import numpy as np
import esda
from pysal.lib import weights
from splot.esda import moran_scatterplot, lisa_cluster, plot_local_autocorrelation
from esda.moran import Moran_Local
from giddy.directional import Rose
from pysal.lib import weights
import libpysal
import splot
from splot.giddy import dynamic_lisa_composite
Load satellite night-time lights data and skim data
satellite = gpd.read_file("satellite.gpkg")
satellite_df = gpd.read_file("satellite.csv")
satellite_df
Plot the 24 countries in East Asia and Pacific
satellite.plot(column="Country",
categorical=True,
legend=True,
figsize=(10, 10),
cmap="Set3")
plt.axis('off')
# save
plt.savefig("Countries.png", dpi = 72)
Compute spatial weights matrices
knn7= weights.KNN.from_dataframe(satellite, k=7)
knn8= weights.KNN.from_dataframe(satellite, k=8)
knn9= weights.KNN.from_dataframe(satellite, k=9)
#queen = weights.Queen.from_dataframe(satellite)
#rook = weights.Rook.from_dataframe(satellite)
#distB = weights.DistanceBand.from_dataframe(satellite, 1182479)
#distC = weights.DistanceBand.from_dataframe(satellite, 1182479, binary=False)
knn7.transform = 'r'
knn8.transform = 'r'
knn9.transform = 'r'
#queen.transform = 'r'
#rook.transform = 'r'
#distB.transform = 'r'
#distC.transform = 'r' #kernel ran out of memory, computed on GeoDa instead
/opt/conda/lib/python3.8/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected:
There are 29 disconnected components.
warnings.warn(message)
/opt/conda/lib/python3.8/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected:
There are 25 disconnected components.
warnings.warn(message)
/opt/conda/lib/python3.8/site-packages/libpysal/weights/weights.py:172: UserWarning: The weights matrix is not fully connected:
There are 21 disconnected components.
warnings.warn(message)
1) VIIRS satellite night-time lights radiance intensity
Compute Global Moran's I and p-values
mi_12_7 = esda.Moran(satellite['Sat_2012'], knn7)
mi_12_7_p = mi_12_7.p_sim
mi_12_7 = mi_12_7.I
mi_12_8 = esda.Moran(satellite['Sat_2012'], knn8)
mi_12_8_p = mi_12_8.p_sim
mi_12_8 = mi_12_8.I
mi_12_9 = esda.Moran(satellite['Sat_2012'], knn9)
mi_12_9_p = mi_12_9.p_sim
mi_12_9 = mi_12_9.I
mi_18_7 = esda.Moran(satellite['Sat_2018'], knn7)
mi_18_7_p = mi_18_7.p_sim
mi_18_7 = mi_18_7.I
mi_18_8 = esda.Moran(satellite['Sat_2018'], knn8)
mi_18_8_p = mi_18_8.p_sim
mi_18_8 = mi_18_8.I
mi_18_9 = esda.Moran(satellite['Sat_2018'], knn9)
mi_18_9_p = mi_18_9.p_sim
mi_18_9 = mi_18_9.I
data = {'Year': ['2012', '2018'],
'7 k-NN': [mi_12_7, mi_18_7], '7 p-value': [mi_12_7_p, mi_18_7_p],
'8 k-NN': [mi_12_8, mi_18_8], '8 p-value': [mi_12_8_p, mi_18_8_p],
'9 k-NN': [mi_12_9, mi_18_9], '9 p-value': [mi_12_9_p, mi_18_9_p]}
df = pd.DataFrame(data)
df
Compute Local Moran's I and p-value
lisa_12 = esda.Moran_Local(satellite['Sat_2012'], knn8, seed = 123456)
lisa_18 = esda.Moran_Local(satellite['Sat_2018'], knn8, seed = 123456)
satellite['sig_12'] = lisa_12.p_sim < 0.01 # 1% as threshold for statistical significance
satellite['sig_18'] = lisa_18.p_sim < 0.01
Compute quadrants
satellite['quad_12'] = lisa_12.q
satellite['quad_18'] = lisa_18.q
Plot Global Moran scatterplot
fig, ax = moran_scatterplot(lisa_12, p=0.01)
plt.show()
# save
plt.savefig("Global Moran 2012", dpi = 72)
fig, ax = moran_scatterplot(lisa_18, p=0.01)
plt.show()
# save
plt.savefig("Global Moran 2018", dpi = 72)
Calculating sum of regions in each quadrants
hh = (satellite['quad_12']==1) & (satellite['sig_12']==True)
hh_12 = hh.sum()
ll = (satellite['quad_12']==3) & (satellite['sig_12']==True)
ll_12 = ll.sum()
lh = (satellite['quad_12']==2) & (satellite['sig_12']==True)
lh_12 = lh.sum()
hl = (satellite['quad_12']==4) & (satellite['sig_12']==True)
hl_12 = hl.sum()
hh = (satellite['quad_18']==1) & (satellite['sig_18']==True)
hh_18 = hh.sum()
ll = (satellite['quad_18']==3) & (satellite['sig_18']==True)
ll_18 = ll.sum()
lh = (satellite['quad_18']==2) & (satellite['sig_18']==True)
lh_18 = lh.sum()
hl = (satellite['quad_18']==4) & (satellite['sig_18']==True)
hl_18 = hl.sum()
data = {'Quadrant': ['High-High', 'Low-Low', 'Low-High', 'High-Low'],
'2012': [hh_12, ll_12, lh_12, hl_12],
'2018': [hh_18, ll_18, lh_18, hl_18],
'Differences': [hh_18-hh_12, ll_18-ll_12, lh_18-lh_12, hl_18-hl_12]}
df = pd.DataFrame(data)
# save
df.to_csv('Quadrants.csv')
df
Standardise variables
satellite['sat_12r'] = (satellite['Sat_2012'] / satellite['Sat_2012'].mean())
satellite['sat_18r'] = (satellite['Sat_2018'] / satellite['Sat_2018'].mean())
Retrieve data for two points in time and calculate rose diagram
y1 = satellite['sat_12r'].values
y2 = satellite['sat_18r'].values
Y = np.array([y1, y2]).T
rose = Rose(Y, knn8, k=5)
Plot Directional LISAs
Conditioned on the starting relative economic density
rose.plot(Y[:,0])
# save
plt.savefig("Directional LISAs - own.png", dpi = 72)
Conditioned on the spatial lag of starting relative economic density
rose.plot(attribute=rose.lag[:,0])
# save
plt.savefig("Directional LISAs - spatial lag.png", dpi = 72)
2) Population density
Standardise Variables
satellite['Pop_10r'] = (satellite['Pop_2010'] / satellite['Pop_2010'].mean())
satellite['Pop_20r'] = (satellite['Pop_2020'] / satellite['Pop_2020'].mean())
Retrieve data for two points in time and calculate rose diagram
y1 = satellite['Pop_10r'].values
y2 = satellite['Pop_20r'].values
Y = np.array([y1, y2]).T
rose = Rose(Y, knn8, k=5)
Plot Directional LISAs
Conditioned on the starting relative economic density
rose.plot(Y[:,0])
# save
plt.savefig("P-Directional LISAs - own.png", dpi = 72)
Conditioned on the spatial lag of starting relative economic density
rose.plot(attribute=rose.lag[:,0])
# save
plt.savefig("P-Directional LISAs - spatial lag.png", dpi = 72)