%matplotlib inline
%reload_ext autoreload
%autoreload 2Laos Evaluation
import sys
sys.path.append("../../../")
import os
from povertymapping import settings
from povertymapping.rollout_grids import get_region_filtered_bingtile_grids, compute_raster_stats
from povertymapping.hdx import get_hdx_file
from povertymapping.iso3 import get_iso3_code
from povertymapping.geoboundaries import get_geoboundaries
import geowrangler.vector_zonal_stats as vzs
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import spearmanr/home/alron/unicef-ai4d-poverty-mapping/env/lib/python3.9/site-packages/geopandas/_compat.py:111: UserWarning: The Shapely GEOS version (3.11.1-CAPI-1.17.1) is incompatible with the GEOS version PyGEOS was compiled with (3.11.2-CAPI-1.17.2). Conversions between both will be slow.
warnings.warn(Set country-specific parameters
REGION = 'laos'
ADMIN_LVL = 'ADM2'
ZOOM_LEVEL = 14
GROUP_ZOOM_LEVEL = 8
country_code = get_iso3_code(REGION, code='alpha-2').lower()Read bounds
bounds_file = get_geoboundaries(region=REGION, adm='ADM1')
bounds_gdf = gpd.read_file(bounds_file)
bounds_gdf = bounds_gdf.sort_values('shapeName')bounds_gdf.sort_values('shapeName')#.explore()| shapeISO | shapeName | shapeID | shapeGroup | shapeType | geometry | |
|---|---|---|---|---|---|---|
| 1 | LA-AT | Attapeu | LAO-ADM1-3_0_0-B2 | LAO | ADM1 | MULTIPOLYGON (((107.59147 14.87336, 107.59141 ... |
| 2 | LA-BK | Bokeo | LAO-ADM1-3_0_0-B3 | LAO | ADM1 | MULTIPOLYGON (((101.25529 20.18622, 101.25468 ... |
| 3 | LA-BL | Bolikhamsai | LAO-ADM1-3_0_0-B4 | LAO | ADM1 | MULTIPOLYGON (((105.25881 18.23862, 105.25859 ... |
| 4 | LA-CH | Champasak | LAO-ADM1-3_0_0-B5 | LAO | ADM1 | MULTIPOLYGON (((106.82797 15.01953, 106.82743 ... |
| 5 | LA-HO | Houaphan | LAO-ADM1-3_0_0-B6 | LAO | ADM1 | MULTIPOLYGON (((104.99274 20.09541, 104.99377 ... |
| 6 | LA-KH | Khammouane | LAO-ADM1-3_0_0-B7 | LAO | ADM1 | MULTIPOLYGON (((106.42514 17.02408, 106.42598 ... |
| 7 | LA-LM | Luang Namtha | LAO-ADM1-3_0_0-B8 | LAO | ADM1 | MULTIPOLYGON (((101.77623 21.12528, 101.77610 ... |
| 8 | LA-LP | Luang Prabang | LAO-ADM1-3_0_0-B9 | LAO | ADM1 | MULTIPOLYGON (((103.40047 20.78283, 103.39479 ... |
| 9 | LA-OU | Oudomxay | LAO-ADM1-3_0_0-B10 | LAO | ADM1 | MULTIPOLYGON (((102.35779 20.85826, 102.35774 ... |
| 10 | LA-PH | Phongsaly | LAO-ADM1-3_0_0-B11 | LAO | ADM1 | MULTIPOLYGON (((102.99287 21.59014, 102.99283 ... |
| 11 | LA-SL | Salavan | LAO-ADM1-3_0_0-B12 | LAO | ADM1 | MULTIPOLYGON (((107.15702 16.25670, 107.15709 ... |
| 12 | LA-SV | Savannakhet | LAO-ADM1-3_0_0-B13 | LAO | ADM1 | MULTIPOLYGON (((106.79847 16.36384, 106.79878 ... |
| 14 | LA-VI | Vientiane | LAO-ADM1-3_0_0-B15 | LAO | ADM1 | MULTIPOLYGON (((102.89857 18.34156, 102.89845 ... |
| 13 | LA-VT | Vientiane Capital | LAO-ADM1-3_0_0-B14 | LAO | ADM1 | MULTIPOLYGON (((103.10281 18.14077, 103.10143 ... |
| 0 | LA-XA | Xaignabouli | LAO-ADM1-3_0_0-B1 | LAO | ADM1 | MULTIPOLYGON (((101.74539 19.88053, 101.74843 ... |
| 15 | LA-XN | Xaisomboun | LAO-ADM1-3_0_0-B16 | LAO | ADM1 | MULTIPOLYGON (((103.89908 18.95856, 103.89927 ... |
| 16 | LA-XE | Xekong | LAO-ADM1-3_0_0-B17 | LAO | ADM1 | MULTIPOLYGON (((107.63500 15.31944, 107.63498 ... |
| 17 | LA-XI | Xiangkhouang | LAO-ADM1-3_0_0-B18 | LAO | ADM1 | MULTIPOLYGON (((104.28077 19.11118, 104.28054 ... |
Read reference data
data_df = pd.read_csv("./la-iwi-2017.csv")
data_df = bounds_gdf.merge(data_df, left_on='shapeName', right_on='province_name', how='left')
data_df['rank'] = data_df['iwi_2017'].rank(ascending=False)
data_df| shapeISO | shapeName | shapeID | shapeGroup | shapeType | geometry | province_name | iwi_2017 | rank | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | LA-AT | Attapeu | LAO-ADM1-3_0_0-B2 | LAO | ADM1 | MULTIPOLYGON (((107.59147 14.87336, 107.59141 ... | Attapeu | 66.0 | 11.0 |
| 1 | LA-BK | Bokeo | LAO-ADM1-3_0_0-B3 | LAO | ADM1 | MULTIPOLYGON (((101.25529 20.18622, 101.25468 ... | Bokeo | 72.0 | 7.0 |
| 2 | LA-BL | Bolikhamsai | LAO-ADM1-3_0_0-B4 | LAO | ADM1 | MULTIPOLYGON (((105.25881 18.23862, 105.25859 ... | Bolikhamsai | 82.4 | 2.0 |
| 3 | LA-CH | Champasak | LAO-ADM1-3_0_0-B5 | LAO | ADM1 | MULTIPOLYGON (((106.82797 15.01953, 106.82743 ... | Champasak | 74.3 | 5.0 |
| 4 | LA-HO | Houaphan | LAO-ADM1-3_0_0-B6 | LAO | ADM1 | MULTIPOLYGON (((104.99274 20.09541, 104.99377 ... | Houaphan | 59.8 | 15.0 |
| 5 | LA-KH | Khammouane | LAO-ADM1-3_0_0-B7 | LAO | ADM1 | MULTIPOLYGON (((106.42514 17.02408, 106.42598 ... | Khammouane | 70.7 | 8.0 |
| 6 | LA-LM | Luang Namtha | LAO-ADM1-3_0_0-B8 | LAO | ADM1 | MULTIPOLYGON (((101.77623 21.12528, 101.77610 ... | Luang Namtha | 68.1 | 9.0 |
| 7 | LA-LP | Luang Prabang | LAO-ADM1-3_0_0-B9 | LAO | ADM1 | MULTIPOLYGON (((103.40047 20.78283, 103.39479 ... | Luang Prabang | 63.1 | 12.0 |
| 8 | LA-OU | Oudomxay | LAO-ADM1-3_0_0-B10 | LAO | ADM1 | MULTIPOLYGON (((102.35779 20.85826, 102.35774 ... | Oudomxay | 62.1 | 13.0 |
| 9 | LA-PH | Phongsaly | LAO-ADM1-3_0_0-B11 | LAO | ADM1 | MULTIPOLYGON (((102.99287 21.59014, 102.99283 ... | Phongsaly | 56.6 | 16.0 |
| 10 | LA-SL | Salavan | LAO-ADM1-3_0_0-B12 | LAO | ADM1 | MULTIPOLYGON (((107.15702 16.25670, 107.15709 ... | Salavan | 56.0 | 17.0 |
| 11 | LA-SV | Savannakhet | LAO-ADM1-3_0_0-B13 | LAO | ADM1 | MULTIPOLYGON (((106.79847 16.36384, 106.79878 ... | Savannakhet | 67.4 | 10.0 |
| 12 | LA-VI | Vientiane | LAO-ADM1-3_0_0-B15 | LAO | ADM1 | MULTIPOLYGON (((102.89857 18.34156, 102.89845 ... | Vientiane | 76.5 | 4.0 |
| 13 | LA-VT | Vientiane Capital | LAO-ADM1-3_0_0-B14 | LAO | ADM1 | MULTIPOLYGON (((103.10281 18.14077, 103.10143 ... | Vientiane Capital | 90.1 | 1.0 |
| 14 | LA-XA | Xaignabouli | LAO-ADM1-3_0_0-B1 | LAO | ADM1 | MULTIPOLYGON (((101.74539 19.88053, 101.74843 ... | Xaignabouli | 77.2 | 3.0 |
| 15 | LA-XN | Xaisomboun | LAO-ADM1-3_0_0-B16 | LAO | ADM1 | MULTIPOLYGON (((103.89908 18.95856, 103.89927 ... | NaN | NaN | NaN |
| 16 | LA-XE | Xekong | LAO-ADM1-3_0_0-B17 | LAO | ADM1 | MULTIPOLYGON (((107.63500 15.31944, 107.63498 ... | Xekong | 61.0 | 14.0 |
| 17 | LA-XI | Xiangkhouang | LAO-ADM1-3_0_0-B18 | LAO | ADM1 | MULTIPOLYGON (((104.28077 19.11118, 104.28054 ... | Xiangkhouang | 73.1 | 6.0 |
Read and process model output
out_gdf = gpd.read_file(settings.ROOT_DIR/"notebooks/2023-02-21-cross-country-rollouts/la/2023-02-21-la-rollout-output.geojson")
out_gdf.head()#explore(column='Predicted Relative Wealth Index')| quadkey | shapeName | shapeISO | shapeID | shapeGroup | shapeType | pop_count | Predicted Relative Wealth Index | Predicted Wealth Category (quintile) | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 13220113003212 | Tonpheung | None | LAO-ADM2-3_0_0-B114 | LAO | ADM2 | 1139.814456 | 0.546717 | A | POLYGON ((100.06348 20.34463, 100.06348 20.365... |
| 1 | 13220113003211 | Tonpheung | None | LAO-ADM2-3_0_0-B114 | LAO | ADM2 | 18.897820 | 0.497037 | A | POLYGON ((100.08545 20.36523, 100.08545 20.385... |
| 2 | 13220113003213 | Tonpheung | None | LAO-ADM2-3_0_0-B114 | LAO | ADM2 | 467.272466 | 0.622080 | A | POLYGON ((100.08545 20.34463, 100.08545 20.365... |
| 3 | 13220113003231 | Tonpheung | None | LAO-ADM2-3_0_0-B114 | LAO | ADM2 | 457.210483 | 0.620122 | A | POLYGON ((100.08545 20.32402, 100.08545 20.344... |
| 4 | 13220113003233 | Tonpheung | None | LAO-ADM2-3_0_0-B114 | LAO | ADM2 | 484.636913 | 0.551952 | A | POLYGON ((100.08545 20.30342, 100.08545 20.324... |
# get stats for provinces (ADM1)
out_gdf = vzs.create_zonal_stats(
bounds_gdf,
out_gdf,
overlap_method="intersects",
aggregations=[{"column":'Predicted Relative Wealth Index', "func": "mean", "output":"model_wealth_index_mean"}]
)
out_gdf['rank'] = out_gdf["model_wealth_index_mean"].rank(ascending=False)
out_gdf| shapeISO | shapeName | shapeID | shapeGroup | shapeType | geometry | model_wealth_index_mean | rank | |
|---|---|---|---|---|---|---|---|---|
| 1 | LA-AT | Attapeu | LAO-ADM1-3_0_0-B2 | LAO | ADM1 | MULTIPOLYGON (((107.59147 14.87336, 107.59141 ... | 0.255099 | 17.0 |
| 2 | LA-BK | Bokeo | LAO-ADM1-3_0_0-B3 | LAO | ADM1 | MULTIPOLYGON (((101.25529 20.18622, 101.25468 ... | 0.281684 | 6.0 |
| 3 | LA-BL | Bolikhamsai | LAO-ADM1-3_0_0-B4 | LAO | ADM1 | MULTIPOLYGON (((105.25881 18.23862, 105.25859 ... | 0.296864 | 3.0 |
| 4 | LA-CH | Champasak | LAO-ADM1-3_0_0-B5 | LAO | ADM1 | MULTIPOLYGON (((106.82797 15.01953, 106.82743 ... | 0.283330 | 5.0 |
| 5 | LA-HO | Houaphan | LAO-ADM1-3_0_0-B6 | LAO | ADM1 | MULTIPOLYGON (((104.99274 20.09541, 104.99377 ... | 0.258527 | 14.0 |
| 6 | LA-KH | Khammouane | LAO-ADM1-3_0_0-B7 | LAO | ADM1 | MULTIPOLYGON (((106.42514 17.02408, 106.42598 ... | 0.290346 | 4.0 |
| 7 | LA-LM | Luang Namtha | LAO-ADM1-3_0_0-B8 | LAO | ADM1 | MULTIPOLYGON (((101.77623 21.12528, 101.77610 ... | 0.275370 | 7.0 |
| 8 | LA-LP | Luang Prabang | LAO-ADM1-3_0_0-B9 | LAO | ADM1 | MULTIPOLYGON (((103.40047 20.78283, 103.39479 ... | 0.261385 | 12.0 |
| 9 | LA-OU | Oudomxay | LAO-ADM1-3_0_0-B10 | LAO | ADM1 | MULTIPOLYGON (((102.35779 20.85826, 102.35774 ... | 0.261008 | 13.0 |
| 10 | LA-PH | Phongsaly | LAO-ADM1-3_0_0-B11 | LAO | ADM1 | MULTIPOLYGON (((102.99287 21.59014, 102.99283 ... | 0.256803 | 16.0 |
| 11 | LA-SL | Salavan | LAO-ADM1-3_0_0-B12 | LAO | ADM1 | MULTIPOLYGON (((107.15702 16.25670, 107.15709 ... | 0.258425 | 15.0 |
| 12 | LA-SV | Savannakhet | LAO-ADM1-3_0_0-B13 | LAO | ADM1 | MULTIPOLYGON (((106.79847 16.36384, 106.79878 ... | 0.272405 | 8.0 |
| 14 | LA-VI | Vientiane | LAO-ADM1-3_0_0-B15 | LAO | ADM1 | MULTIPOLYGON (((102.89857 18.34156, 102.89845 ... | 0.305262 | 2.0 |
| 13 | LA-VT | Vientiane Capital | LAO-ADM1-3_0_0-B14 | LAO | ADM1 | MULTIPOLYGON (((103.10281 18.14077, 103.10143 ... | 0.415501 | 1.0 |
| 0 | LA-XA | Xaignabouli | LAO-ADM1-3_0_0-B1 | LAO | ADM1 | MULTIPOLYGON (((101.74539 19.88053, 101.74843 ... | 0.271010 | 9.0 |
| 15 | LA-XN | Xaisomboun | LAO-ADM1-3_0_0-B16 | LAO | ADM1 | MULTIPOLYGON (((103.89908 18.95856, 103.89927 ... | 0.264107 | 11.0 |
| 16 | LA-XE | Xekong | LAO-ADM1-3_0_0-B17 | LAO | ADM1 | MULTIPOLYGON (((107.63500 15.31944, 107.63498 ... | 0.249704 | 18.0 |
| 17 | LA-XI | Xiangkhouang | LAO-ADM1-3_0_0-B18 | LAO | ADM1 | MULTIPOLYGON (((104.28077 19.11118, 104.28054 ... | 0.266895 | 10.0 |
Compare model output vs reference
compare_gdf = out_gdf[['shapeISO','shapeName','model_wealth_index_mean','rank', 'geometry']].merge(\
data_df[['shapeISO','shapeName','iwi_2017','rank']], on='shapeISO', how='outer', suffixes=['_model', '_reference'])#.dropna()
compare_gdf| shapeISO | shapeName_model | model_wealth_index_mean | rank_model | geometry | shapeName_reference | iwi_2017 | rank_reference | |
|---|---|---|---|---|---|---|---|---|
| 0 | LA-AT | Attapeu | 0.255099 | 17.0 | MULTIPOLYGON (((107.59147 14.87336, 107.59141 ... | Attapeu | 66.0 | 11.0 |
| 1 | LA-BK | Bokeo | 0.281684 | 6.0 | MULTIPOLYGON (((101.25529 20.18622, 101.25468 ... | Bokeo | 72.0 | 7.0 |
| 2 | LA-BL | Bolikhamsai | 0.296864 | 3.0 | MULTIPOLYGON (((105.25881 18.23862, 105.25859 ... | Bolikhamsai | 82.4 | 2.0 |
| 3 | LA-CH | Champasak | 0.283330 | 5.0 | MULTIPOLYGON (((106.82797 15.01953, 106.82743 ... | Champasak | 74.3 | 5.0 |
| 4 | LA-HO | Houaphan | 0.258527 | 14.0 | MULTIPOLYGON (((104.99274 20.09541, 104.99377 ... | Houaphan | 59.8 | 15.0 |
| 5 | LA-KH | Khammouane | 0.290346 | 4.0 | MULTIPOLYGON (((106.42514 17.02408, 106.42598 ... | Khammouane | 70.7 | 8.0 |
| 6 | LA-LM | Luang Namtha | 0.275370 | 7.0 | MULTIPOLYGON (((101.77623 21.12528, 101.77610 ... | Luang Namtha | 68.1 | 9.0 |
| 7 | LA-LP | Luang Prabang | 0.261385 | 12.0 | MULTIPOLYGON (((103.40047 20.78283, 103.39479 ... | Luang Prabang | 63.1 | 12.0 |
| 8 | LA-OU | Oudomxay | 0.261008 | 13.0 | MULTIPOLYGON (((102.35779 20.85826, 102.35774 ... | Oudomxay | 62.1 | 13.0 |
| 9 | LA-PH | Phongsaly | 0.256803 | 16.0 | MULTIPOLYGON (((102.99287 21.59014, 102.99283 ... | Phongsaly | 56.6 | 16.0 |
| 10 | LA-SL | Salavan | 0.258425 | 15.0 | MULTIPOLYGON (((107.15702 16.25670, 107.15709 ... | Salavan | 56.0 | 17.0 |
| 11 | LA-SV | Savannakhet | 0.272405 | 8.0 | MULTIPOLYGON (((106.79847 16.36384, 106.79878 ... | Savannakhet | 67.4 | 10.0 |
| 12 | LA-VI | Vientiane | 0.305262 | 2.0 | MULTIPOLYGON (((102.89857 18.34156, 102.89845 ... | Vientiane | 76.5 | 4.0 |
| 13 | LA-VT | Vientiane Capital | 0.415501 | 1.0 | MULTIPOLYGON (((103.10281 18.14077, 103.10143 ... | Vientiane Capital | 90.1 | 1.0 |
| 14 | LA-XA | Xaignabouli | 0.271010 | 9.0 | MULTIPOLYGON (((101.74539 19.88053, 101.74843 ... | Xaignabouli | 77.2 | 3.0 |
| 15 | LA-XN | Xaisomboun | 0.264107 | 11.0 | MULTIPOLYGON (((103.89908 18.95856, 103.89927 ... | Xaisomboun | NaN | NaN |
| 16 | LA-XE | Xekong | 0.249704 | 18.0 | MULTIPOLYGON (((107.63500 15.31944, 107.63498 ... | Xekong | 61.0 | 14.0 |
| 17 | LA-XI | Xiangkhouang | 0.266895 | 10.0 | MULTIPOLYGON (((104.28077 19.11118, 104.28054 ... | Xiangkhouang | 73.1 | 6.0 |
data = compare_gdf.dropna()
coef, p = spearmanr(data['iwi_2017'], data['model_wealth_index_mean'])
print(f'The spearman rank correlation is: {coef:.4f}')The spearman rank correlation is: 0.8382Heatmap
def despine(ax):
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
ax.spines["left"].set_visible(False)
ax.spines["bottom"].set_visible(False)
ax.set_xticks([])
ax.set_yticks([])fig, axs = plt.subplots(1,2, figsize=(9,5))
compare_gdf.plot(
column='rank_reference',
ax=axs[0],
cmap="RdYlGn_r",
legend=True,
legend_kwds={"shrink": 0.3},
missing_kwds={
"color": "0.4",
"label": "Missing values",
}
)
despine(axs[0])
compare_gdf.plot(
column='rank_model',
ax=axs[1],
cmap="RdYlGn_r",
legend=True,
legend_kwds={"shrink": 0.3},
missing_kwds={
"color": "0.4",
"label": "Missing values",
}
)
despine(axs[1])
axs[0].set_title('Reference')
axs[1].set_title('Predicted')
# leg1 = axs[0].get_legend()
# leg1.set_title("Rank")
# leg2 = axs[1].get_legend()
# leg2.set_title("Rank")#
#fig.suptitle('Wealth index rank')
plt.savefig(fname="la_rank_comparison.png", dpi=300)
fig.tight_layout()