A geocomputational notebook to monitor regional development in Bolivia

Select indicator

poly_id0.7%

asdf_id0.7%

137 others98.6%

Polygon ID0.7%

ASDF ID0.7%

137 others98.6%

0

poly_id

Polygon ID

1

asdf_id

ASDF ID

2

mun

Municipality

3

mun_id

Municipality ID

4

dep

Department

5

dep_id

Department ID

6

dep_mun

Department-municipality

7

shapeID

Municipality Geoquery Polygon ID

8

imds

Municipal Sustainable Development Index

9

rank_imds

Bolivia Index Ranking

10

population_2020

Population 2020

11

urbano_2012

Urbanization rate, 2012 (% of population)

12

sdg1_1_eepr

Extreme energy poverty rate, 2016 (% of houses)

13

sdg1_1_ubn

Unsatisfied basic needs, 2012 (% of population)

14

sdg1_2_mpi

Multidimensional poverty index, 2012

15

sdg1_4_abs

Access to the 3 basic services, 2012 (% of households)

16

sdg2_2_cmc

Chronic malnutrition in children (< 5 years), 2016 (%)

17

sdg2_2_oww

Obesity in women (15-49 years), 2016 (%)

18

sdg2_4_pual

Average area per Production Unit Agriculture and Livestock, 2013 (ha)

19

sdg2_4_td

Tractor density, 2013 (per 1,000 UPAs)

20

sdg3_1_idca

Institutional childbirth coverage, average 2008-2012 (%)

21

sdg3_2_imr

Infant mortality rate (< 1 year), 2016 (per. 1,000 live births)

22

sdg3_2_mrc

Children mortality rate in (< 5 years), 2016 (per. 1,000 live births)

23

sdg3_3_cdir

Chagas disease infestation rate, 2017 (% of households)

24

sdg3_3_di

Dengue incidence, 2018 (per 10,000 population)

139 rows, showing

per page

Page of 6

variables = dataDefinitions['Variable']

INDICATOR1 = 'imds' INDICATOR2 = 'pop2017' INDICATOR3 = 'ln_t400NTLpc2017' INDICATOR4 = 'co2017'

Descriptive statistics

gdf[[INDICATOR1, INDICATOR2, INDICATOR3, INDICATOR4]].describe().round(2)

count

339

mean

51.05

32858.76

std

6.77

117649.33

min

35.7

661.82

25%

47

6425.71

50%

50.5

11627.51

75%

54.85

22497.26

max

80.2

1604457.5

8 rows, showing

per page

Page of 1

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 bandwidth = 46.0

gwr_results = GWR(coords, y, x, gwr_bw).fit() gwr_results.summary()

===========================================================================
Model type                                                         Gaussian
Number of observations:                                                 339
Number of covariates:                                                     2

Global Regression Results
---------------------------------------------------------------------------
Residual sum of squares:                                          10221.458
Log-likelihood:                                                   -1058.379
AIC:                                                               2120.757
AICc:                                                              2122.829
BIC:                                                               8258.096
R2:                                                                   0.341
Adj. R2:                                                              0.339

Variable                              Est.         SE  t(Est/SE)    p-value
------------------------------- ---------- ---------- ---------- ----------
X0                                   8.808      3.216      2.739      0.006
X1                                   3.091      0.234     13.192      0.000

Geographically Weighted Regression (GWR) Results
---------------------------------------------------------------------------
Spatial kernel:                                           Adaptive bisquare
Bandwidth used:                                                      46.000

Diagnostic information
---------------------------------------------------------------------------
Residual sum of squares:                                           6137.357
Effective number of parameters (trace(S)):                           37.677
Degree of freedom (n - trace(S)):                                   301.323
Sigma estimate:                                                       4.513
Log-likelihood:                                                    -971.917
AIC:                                                               2021.188
AICc:                                                              2031.442
BIC:                                                               2169.165
R2:                                                                   0.604
Adjusted R2:                                                          0.554
Adj. alpha (95%):                                                     0.003
Adj. critical t value (95%):                                          3.028

Summary Statistics For GWR Parameter Estimates
---------------------------------------------------------------------------
Variable                   Mean        STD        Min     Median        Max
-------------------- ---------- ---------- ---------- ---------- ----------
X0                        5.465     26.130    -56.622     13.003     62.759
X1                        3.314      1.811     -0.689      2.819      7.676
===========================================================================

# 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)

Mean R2 = 0.6040381097677936
AIC = 2021.1882671117585
AICc = 2031.441802189135

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