Predict on rollout grids

%matplotlib inline
%reload_ext autoreload
%autoreload 2

import os
import sys

sys.path.append("../../../")

import getpass
import pickle
from pathlib import Path

import contextily as cx
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from povertymapping import nightlights, settings
from povertymapping.dhs import generate_dhs_cluster_level_data
from povertymapping.feature_engineering import (
    categorize_wealth_index,
    generate_features,
)
from povertymapping.iso3 import get_region_name
from povertymapping.rollout_grids import get_region_filtered_bingtile_grids

Model Prediction on Rollout Grids: Myanmar

This notebook is the final step in the rollout and runs the final model to create relative wealth estimations over populated areas within the given country. The model predictions will have a spatial resolution of 2.4km.

The predicted relative wealth value gives us the relative wealth level of an area compared to the rest of the country, which fixes the value range from 0 (lowest wealth) to 1 (highest wealth). In between these extremes, each area’s wealth estimate is scaled to a value between 0 and 1.

The predicted relative wealth value is later binned into 5 wealth categories A-E by dividing the distribution into quintiles (every 20th percentile).

Set up Data Access

The following cell will prompt you to enter your EOG username and password. See this page to learn how to set-up your EOG account.

# Log-in using EOG credentials
username = os.environ.get("EOG_USER", None)
username = username if username is not None else input("Username?")
password = os.environ.get("EOG_PASSWORD", None)
password = password if password is not None else getpass.getpass("Password?")

# set save_token to True so that access token gets stored in ~/.eog_creds/eog_access_token
access_token = nightlights.get_eog_access_token(username, password, save_token=True)

2023-04-16 23:00:54.461 | INFO     | povertymapping.nightlights:get_eog_access_token:43 - Loaded access_token from /home/alron/.eog_creds/eog_access_token.txt

Set country-specific parameters

COUNTRY_CODE = "mm"
COUNTRY_OSM = get_region_name(COUNTRY_CODE, code="alpha-2").lower()
OOKLA_YEAR = 2019
NIGHTLIGHTS_YEAR = 2016

rollout_date = "-".join(os.getcwd().split("/")[-2].split("-")[:3])
rollout_grids_path = Path(f"./{rollout_date}-{COUNTRY_CODE}-rollout-grids.geojson")
rollout_grids_path

Path('2023-02-21-mm-rollout-grids.geojson')

Set Model Parameters

# Model to use for prediction
MODEL_SAVE_PATH = Path(f"./{rollout_date}-{COUNTRY_CODE}-single-country-model.pkl")

Load Country Rollout AOI

The rollout area of interest is split into 2.4km grid tiles (zoom level 14), matching the areas used during model training. The grids are also filtered to only include populated areas based on Meta’s High Resolution Settlement Layer (HRSL) data.

Refer to the previous notebook 2_mm_generate_grids.ipynb for documentation on generating this grid.

aoi = gpd.read_file(rollout_grids_path)
# aoi.explore()  # Uncomment to view data in a map

Generate Features For Rollout AOI

%%time
rollout_aoi = aoi.copy()

# Create features dataframe using generate_features module
features = generate_features(
    rollout_aoi,
    country_osm=COUNTRY_OSM,
    ookla_year=OOKLA_YEAR,
    nightlights_year=NIGHTLIGHTS_YEAR,
    scale=False,
    features_only=True,
)

2023-04-16 23:00:59.368 | INFO     | povertymapping.osm:download_osm_country_data:199 - OSM Data: Cached data available for myanmar at /home/alron/.geowrangler/osm/myanmar? True
2023-04-16 23:00:59.369 | DEBUG    | povertymapping.osm:load_pois:161 - OSM POIs for myanmar being loaded from /home/alron/.geowrangler/osm/myanmar/gis_osm_pois_free_1.shp

Inspect the generated features

features.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
Int64Index: 27830 entries, 0 to 27829
Data columns (total 61 columns):
 #   Column                             Non-Null Count  Dtype  
---  ------                             --------------  -----  
 0   poi_count                          27830 non-null  float64
 1   atm_count                          27830 non-null  float64
 2   atm_nearest                        27830 non-null  float64
 3   bank_count                         27830 non-null  float64
 4   bank_nearest                       27830 non-null  float64
 5   bus_station_count                  27830 non-null  float64
 6   bus_station_nearest                27830 non-null  float64
 7   cafe_count                         27830 non-null  float64
 8   cafe_nearest                       27830 non-null  float64
 9   charging_station_count             27830 non-null  float64
 10  charging_station_nearest           27830 non-null  float64
 11  courthouse_count                   27830 non-null  float64
 12  courthouse_nearest                 27830 non-null  float64
 13  dentist_count                      27830 non-null  float64
 14  dentist_nearest                    27830 non-null  float64
 15  fast_food_count                    27830 non-null  float64
 16  fast_food_nearest                  27830 non-null  float64
 17  fire_station_count                 27830 non-null  float64
 18  fire_station_nearest               27830 non-null  float64
 19  food_court_count                   27830 non-null  float64
 20  food_court_nearest                 27830 non-null  float64
 21  fuel_count                         27830 non-null  float64
 22  fuel_nearest                       27830 non-null  float64
 23  hospital_count                     27830 non-null  float64
 24  hospital_nearest                   27830 non-null  float64
 25  library_count                      27830 non-null  float64
 26  library_nearest                    27830 non-null  float64
 27  marketplace_count                  27830 non-null  float64
 28  marketplace_nearest                27830 non-null  float64
 29  pharmacy_count                     27830 non-null  float64
 30  pharmacy_nearest                   27830 non-null  float64
 31  police_count                       27830 non-null  float64
 32  police_nearest                     27830 non-null  float64
 33  post_box_count                     27830 non-null  float64
 34  post_box_nearest                   27830 non-null  float64
 35  post_office_count                  27830 non-null  float64
 36  post_office_nearest                27830 non-null  float64
 37  restaurant_count                   27830 non-null  float64
 38  restaurant_nearest                 27830 non-null  float64
 39  social_facility_count              27830 non-null  float64
 40  social_facility_nearest            27830 non-null  float64
 41  supermarket_count                  27830 non-null  float64
 42  supermarket_nearest                27830 non-null  float64
 43  townhall_count                     27830 non-null  float64
 44  townhall_nearest                   27830 non-null  float64
 45  road_count                         27830 non-null  float64
 46  fixed_2019_mean_avg_d_kbps_mean    27830 non-null  float64
 47  fixed_2019_mean_avg_u_kbps_mean    27830 non-null  float64
 48  fixed_2019_mean_avg_lat_ms_mean    27830 non-null  float64
 49  fixed_2019_mean_num_tests_mean     27830 non-null  float64
 50  fixed_2019_mean_num_devices_mean   27830 non-null  float64
 51  mobile_2019_mean_avg_d_kbps_mean   27830 non-null  float64
 52  mobile_2019_mean_avg_u_kbps_mean   27830 non-null  float64
 53  mobile_2019_mean_avg_lat_ms_mean   27830 non-null  float64
 54  mobile_2019_mean_num_tests_mean    27830 non-null  float64
 55  mobile_2019_mean_num_devices_mean  27830 non-null  float64
 56  avg_rad_min                        27830 non-null  float64
 57  avg_rad_max                        27830 non-null  float64
 58  avg_rad_mean                       27830 non-null  float64
 59  avg_rad_std                        27830 non-null  float64
 60  avg_rad_median                     27830 non-null  float64
dtypes: float64(61)
memory usage: 14.2 MB

Run Model on AOI

Load Model

with open(MODEL_SAVE_PATH, "rb") as f:
    model = pickle.load(f)

Make Predictions

rollout_aoi["Predicted Relative Wealth Index"] = model.predict(features.values)

Binning predictions into wealth categories

Afterwards, we label the predicted relative wealth by binning them into 5 categories: A, B, C, D, and E where A is the highest and E is the lowest.

We can create these wealth categories by splitting the output Predicted Relative Wealth Index distribution into 5 equally sized quintiles, i.e. every 20th percentile.

This categorization may be modified to suit the context of the target country.

rollout_aoi["Predicted Wealth Category (quintile)"] = categorize_wealth_index(
    rollout_aoi["Predicted Relative Wealth Index"], split_quantile=False
).astype(str)

Save Output

%%time
rollout_aoi.to_file(
    f"{rollout_date}-{COUNTRY_CODE}-rollout-output.geojson",
    driver="GeoJSON",
    index=False,
)

CPU times: user 31.6 s, sys: 656 ms, total: 32.3 s
Wall time: 33 s

# Join back raw features and save
rollout_output_with_features = rollout_aoi.join(features)
rollout_output_with_features.to_file(
    f"{rollout_date}-{COUNTRY_CODE}-rollout-output-with-features.geojson",
    driver="GeoJSON",
    index=False,
)

Visualizations

Inspect predicted wealth index and output dataframe

rollout_aoi[["Predicted Relative Wealth Index"]].hist()

array([[<AxesSubplot: title={'center': 'Predicted Relative Wealth Index'}>]],
      dtype=object)

rollout_aoi.head()

	quadkey	shapeName	shapeISO	shapeID	shapeGroup	shapeType	quadkey_level10	pop_count	geometry	Predicted Relative Wealth Index	...	mobile_2019_mean_avg_d_kbps_mean	mobile_2019_mean_avg_u_kbps_mean	mobile_2019_mean_avg_lat_ms_mean	mobile_2019_mean_num_tests_mean	mobile_2019_mean_num_devices_mean	avg_rad_min	avg_rad_max	avg_rad_mean	avg_rad_std	avg_rad_median
0	13220212231203	Yangon (South)	None	MMR-ADM2-3_0_0-B72	MMR	ADM2	1322021223	1646.340576	POLYGON ((93.36182 14.13658, 93.36182 14.15788...	0.291516	...	0.0	0.0	0.0	0.0	0.0	-0.049524	0.030141	-0.017631	0.019043	-0.020764
1	13220212231221	Yangon (South)	None	MMR-ADM2-3_0_0-B72	MMR	ADM2	1322021223	749.037048	POLYGON ((93.36182 14.11527, 93.36182 14.13658...	0.313033	...	0.0	0.0	0.0	0.0	0.0	-0.023038	0.057469	0.005178	0.020869	0.000433
2	13220030031312	Kyaukpyu	None	MMR-ADM2-3_0_0-B21	MMR	ADM2	1322003003	1276.558105	POLYGON ((93.47168 18.85431, 93.47168 18.87510...	0.313332	...	0.0	0.0	0.0	0.0	0.0	-0.061539	-0.002709	-0.027197	0.012585	-0.025633
3	13220030031313	Kyaukpyu	None	MMR-ADM2-3_0_0-B21	MMR	ADM2	1322003003	59.303814	POLYGON ((93.49365 18.85431, 93.49365 18.87510...	0.302049	...	0.0	0.0	0.0	0.0	0.0	-0.054370	0.008868	-0.024229	0.013000	-0.022040
4	13220030031333	Kyaukpyu	None	MMR-ADM2-3_0_0-B21	MMR	ADM2	1322003003	674.946777	POLYGON ((93.49365 18.81272, 93.49365 18.83352...	0.313332	...	0.0	0.0	0.0	0.0	0.0	-0.053165	-0.007529	-0.028266	0.011237	-0.029542

5 rows × 72 columns

Create Static Maps

Plot Predicted Relative Wealth Index

plt.cla()
plt.clf()
rollout_aoi_plot = rollout_aoi.to_crs("EPSG:3857")
ax = rollout_aoi_plot.plot(
    "Predicted Relative Wealth Index",
    figsize=(20, 8),
    cmap="viridis",
    legend=True,
    legend_kwds={"shrink": 0.8},
)
cx.add_basemap(ax, source=cx.providers.OpenStreetMap.Mapnik)
ax.set_axis_off()
plt.title("Predicted Relative Wealth Index")
plt.tight_layout()
plt.savefig(f"{rollout_date}-{COUNTRY_CODE}-predicted-wealth-index.png")
plt.show()

<Figure size 640x480 with 0 Axes>

Plot Predicted Relative Wealth Index Category

plt.cla()
plt.clf()
rollout_aoi_plot = rollout_aoi.to_crs("EPSG:3857")
ax = rollout_aoi_plot.plot(
    "Predicted Wealth Category (quintile)",
    figsize=(20, 8),
    cmap="viridis_r",
    legend=True,
)
cx.add_basemap(ax, source=cx.providers.OpenStreetMap.Mapnik)
ax.set_axis_off()
plt.title("Predicted Wealth Category")
plt.tight_layout()
plt.savefig(f"{rollout_date}-{COUNTRY_CODE}-predicted-wealth-bin.png")
plt.show()

<Figure size 640x480 with 0 Axes>

Create an Interactive Map

cols_of_interest = [
    "quadkey",
    "shapeName",
    "shapeGroup",
    "pop_count",
    "avg_rad_mean",
    "mobile_2019_mean_avg_d_kbps_mean",
    "fixed_2019_mean_avg_d_kbps_mean",
    "poi_count",
    "road_count",
    "Predicted Relative Wealth Index",
    "Predicted Wealth Category (quintile)",
]

# Warning: This can be a bit laggy due to the large amount of tiles being visualized

# Uncomment the ff if you want to viz the raw wealth predictions
# rollout_aoi.explore(column='Predicted Relative Wealth Index', tooltip=cols_of_interest, cmap="viridis")

# Uncomment the ff if you want to view the quintiles
# rollout_aoi.explore(column='Predicted Wealth Category (quintile)', tooltip=cols_of_interest, cmap="viridis_r")

Alternatively, you may also try to visualize this interactively in Kepler by uploading the rollout output geojson file.