A geocomputational notebook to monitor regional development in Bolivia
Select indicator
variables = dataDefinitions['Variable']
INDICATOR1 = 'imds'
INDICATOR2 = 'pop2017'
INDICATOR3 = 'ln_t400NTLpc2017'
INDICATOR4 = 'co2017'
INDICATOR1
imds
INDICATOR2
pop2017
INDICATOR3
ln_t400NTLpc2017
INDICATOR4
co2017
Descriptive statistics
gdf[[INDICATOR1, INDICATOR2, INDICATOR3, INDICATOR4]].describe().round(2)
Regional distribution
px.strip(gdf, x = INDICATOR1, y = 'dep', color = 'dep', hover_name= 'mun', hover_data= ['rank_imds'],
labels=dict(imds = 'Sustinable development index',
rank_imds = 'Development ranking',
mun = 'Municipality',
dep = 'Department')
)
px.box(gdf, x = INDICATOR1, y = 'dep', color = 'dep', hover_name= 'mun', hover_data= ['rank_imds'],
labels=dict(imds = 'Sustinable development index',
rank_imds = 'Development ranking',
mun = 'Municipality',
dep = 'Department')
)
px.histogram(gdf, x = INDICATOR1, color = 'dep', hover_name= 'mun', marginal='rug', hover_data= ['rank_imds'],
labels=dict(imds = 'Sustinable development index',
rank_imds = 'Development ranking',
mun = 'Municipality',
dep = 'Department')
)
The role of population
px.treemap(gdf, color = INDICATOR1, values = "pop2020", path = ["dep", "mun"], hover_name = "mun",
hover_data= ['rank_imds'],
labels=dict(imds = 'Development',
rank_imds = 'Development ranking',
pop2020 = 'Population in 2020',
mun = 'Municipality',
dep = 'Department')
)
px.sunburst(gdf, color = INDICATOR1, values = "pop2020", path = ["dep", "mun"], hover_name = "mun",
hover_data= ['rank_imds'],
labels=dict(imds = 'Development',
rank_imds = 'Development ranking',
pop2020 = 'Population in 2020',
mun = 'Municipality',
dep = 'Department')
)
Nighttime lights and development
px.scatter(gdf,
x = INDICATOR3,
y = INDICATOR1,
color = 'dep',
symbol = 'dep',
hover_name= 'mun',
trendline= 'ols',
trendline_scope= 'overall',
hover_data= ['rank_imds'],
labels=dict(imds = 'Sustainable development index',
rank_imds = 'Development ranking',
ln_t400NTLpc2017 = 'Trend (log) nighttime lights per capita',
mun = 'Municipality',
dep = 'Department')
)
px.scatter(gdf,
x = INDICATOR3,
y = INDICATOR1,
color = 'dep',
symbol = 'dep',
hover_name= 'mun',
trendline= 'ols',
marginal_x="box",
marginal_y="box",
hover_data= ['rank_imds'],
labels=dict(imds = 'Sustainable development index',
rank_imds = 'Development ranking',
ln_t400NTLpc2017 = 'Trend (log) nighttime lights per capita',
mun = 'Municipality',
dep = 'Department')
)
Geographic distribution
fig, ax = plt.subplots(figsize=(12,8))
gdf.plot(column= INDICATOR1,
scheme='FisherJenks',
k = 3,
cmap='coolwarm',
edgecolor='k',
linewidth=0.5,
alpha=0.8,
legend=True,
ax=ax,
legend_kwds={'bbox_to_anchor':(1.00, 0.92)})
cx.add_basemap(ax, crs=gdf.crs.to_string(), source = cx.providers.Stamen.TonerHybrid, attribution=False)
plt.title('Spatial distribution of regional development: Three natural breaks')
plt.tight_layout()
ax.axis("off")
#plt.savefig('myMap.png',dpi=150, bbox_inches='tight')
plt.show()
gdf.explore(
column= INDICATOR1,
tooltip=['dep', 'mun', 'imds', 'rank_imds'],
k = 3,
scheme='FisherJenks',
cmap='coolwarm',
legend=True,
tiles= 'Stamen Terrain',
style_kwds =dict(color="gray", weight=0.4),
legend_kwds=dict(colorbar=False)
)
Spatial dependence
#W = weights.Queen.from_dataframe(gdf)
W = weights.KNN.from_dataframe(gdf, k=4)
W.transform = 'r'
plot_spatial_weights(W, gdf);
myLIST = ['asdf_id', 'mun', 'dep', 'dep_mun', INDICATOR1]
df_MORAN = gdf[myLIST]
df_MORAN["WxINDICATOR1"] = weights.lag_spatial(W, df_MORAN.iloc[: , -1])
globalMoran = Moran(gdf[INDICATOR1], W)
MoranI = globalMoran.I
MoranI = "{:.2f}".format(MoranI)
localMoran = Moran_Local(gdf[INDICATOR1], W, permutations = 999, seed=12345)
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12,5))
moran_scatterplot(localMoran, p=0.10, aspect_equal=False, zstandard = False, ax=axes[0])
lisa_cluster(localMoran, gdf, p=0.10, legend_kwds={'bbox_to_anchor':(0.02, 0.90)}, ax=axes[1])
#cx.add_basemap(axes[1], crs=gdf.crs.to_string(), source=cx.providers.Stamen.TonerLite, attribution=False)
#cx.add_basemap(axes[1], crs=gdf.crs.to_string(), source=cx.providers.Stamen.TonerLabels, attribution=False)
axes[0].set_xlabel('Indicator')
axes[0].set_ylabel('Spatial lag of indicator')
axes[0].set_title(f"(a) Moran scatterplot (Moran's I = {MoranI})")
axes[1].set_title("(b) Spatial clusters and outliers (p<0.10)")
plt.show()
Spatial inequality
sns.histplot(x=gdf[INDICATOR1], kde=True);
sns.histplot(x=gdf[INDICATOR3], kde=True);
inequality.gini.Gini(gdf[INDICATOR1].values).g
inequality.gini.Gini(gdf[INDICATOR3].values).g
Gini_Spatial(gdf[INDICATOR1], W).wcg_share
Gini_Spatial(gdf[INDICATOR1], W).p_sim
Gini_Spatial(gdf[INDICATOR3], W).wcg_share
Gini_Spatial(gdf[INDICATOR3], W).p_sim
Spatial heterogeneity
y = gdf[INDICATOR1].values.reshape((-1,1))
x = gdf[INDICATOR3].values.reshape((-1,1))
u = gdf['COORD_X']
v = gdf['COORD_Y']
coords = list(zip(u,v))
gwr_selector = Sel_BW(coords, y, x, spherical = True)
gwr_bw = gwr_selector.search(criterion='AICc')
print('GWR bandwidth =', gwr_bw)
gwr_results = GWR(coords, y, x, gwr_bw).fit()
gwr_results.summary()
# As reference, here is the (mean) R2, AIC, and AICc
print('Mean R2 =', gwr_results.R2)
print('AIC =', gwr_results.aic)
print('AICc =', gwr_results.aicc)
# Add R2 to GeoDataframe
gdf['gwr_R2'] = gwr_results.localR2
# Local R2
gdf.explore(
column= 'gwr_R2',
tooltip=['dep', 'mun', 'gwr_R2', 'rank_imds', INDICATOR1, INDICATOR3],
k = 5,
scheme='FisherJenks',
cmap='coolwarm',
legend=True,
tiles= 'CartoDB dark_matter',
style_kwds =dict(color="gray", weight=0.4),
legend_kwds=dict(colorbar=False)
)
# Add coefficients to data frame
gdf['gwr_intercept'] = gwr_results.params[:,0]
gdf['gwr_slope'] = gwr_results.params[:,1]
# Filter t-values: standard alpha = 0.05
gwr_filtered_t = gwr_results.filter_tvals(alpha = 0.05)
# Filter t-values: corrected alpha due to multiple testing
gwr_filtered_tc = gwr_results.filter_tvals()
# Slope heterogeneity
gdf.explore(
column= 'gwr_slope',
tooltip=['dep', 'mun', 'gwr_slope', 'rank_imds', INDICATOR1, INDICATOR3],
k = 5,
scheme='FisherJenks',
cmap='coolwarm',
legend=True,
tiles= 'CartoDB dark_matter',
style_kwds =dict(color="gray", weight=0.4),
legend_kwds=dict(colorbar=False)
)
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18,6))
gdf.plot(column='gwr_slope', cmap = 'coolwarm', linewidth=0.01, scheme = 'FisherJenks', k=5, legend=True, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[0])
gdf.plot(column='gwr_slope', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[1])
gdf[gwr_filtered_t[:,1] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[1])
gdf.plot(column='gwr_slope', cmap = 'coolwarm', linewidth=0.05, scheme = 'FisherJenks', k=5, legend=False, legend_kwds={'bbox_to_anchor':(1.10, 0.96)}, ax=axes[2])
gdf[gwr_filtered_tc[:,1] == 0].plot(color='white', linewidth=0.05, edgecolor='black', ax=axes[2])
plt.tight_layout()
axes[0].axis("off")
axes[1].axis("off")
axes[2].axis("off")
axes[0].set_title('(a) GWR: Nighttime lights (BW: ' + str(gwr_bw) +'), all coeffs', fontsize=12)
axes[1].set_title('(b) GWR: Nighttime lights (BW: ' + str(gwr_bw) +'), significant coeffs', fontsize=12)
axes[2].set_title('(c) GWR: Nighttime lights (BW: ' + str(gwr_bw) +'), significant coeffs and corr. p-values', fontsize=12)
plt.show()