Interactive Geographic Analysis: Population Density and GDP per Person

Purpose. This notebook is a reproducible teaching example for mapping country-level geographic data. It asks:

How does estimated GDP per person relate to population density across countries, and what spatial patterns are visible on a map?

The notebook follows a complete analysis workflow: define the question, load real data, clean and validate it, explore it with tables and interactive maps, fit a simple descriptive model, and interpret the result without treating correlation as causation.

What we use

• Natural Earth Admin 0 countries, 1:110m for country polygons and published estimates of population and GDP.
• A projected, equal-area coordinate reference system (EPSG:6933) for area and density calculations.
• geopandas for spatial data, folium for interactive maps, and scikit-learn for a transparent log-log regression.

The source is downloaded from GitHub at run time and cached locally under data/, so a reader can rerun the notebook and inspect the raw input.

1. Assumptions and analysis choices

1. A country polygon is treated as the unit of analysis; this is an ecological, country-level comparison, not an individual-level study.
2. Natural Earth's POP_EST and GDP_MD_EST are estimates, not a current official statistical release. GDP is in millions of US dollars and population is a count.
3. Population density is computed as estimated population divided by polygon area in square kilometres.
4. Because both density and GDP per person are strongly right-skewed, the model uses natural logarithms. Countries with non-positive or missing values cannot be log-transformed and are excluded from the model with an explicit count reported below.
5. The map is for exploration. It is not evidence of causal effects, and country boundaries and estimates have uncertainty.

# 2. Reproducible setup
from pathlib import Path
import json
import warnings

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import seaborn as sns
import folium
from branca.colormap import linear
from IPython.display import display, Markdown
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_absolute_error

warnings.filterwarnings("ignore", category=UserWarning)
sns.set_theme(style="whitegrid", context="notebook")

PROJECT_ROOT = Path.cwd()
DATA_DIR = PROJECT_ROOT / "data"
OUTPUT_DIR = PROJECT_ROOT / "outputs"
DATA_DIR.mkdir(exist_ok=True)
OUTPUT_DIR.mkdir(exist_ok=True)

SOURCE_URL = (
    "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/"
    "master/geojson/ne_110m_admin_0_countries.geojson"
)
RAW_PATH = DATA_DIR / "ne_110m_admin_0_countries.geojson"
print(f"Project root: {PROJECT_ROOT.resolve()}")
print(f"Source: {SOURCE_URL}")

Project root: /home/ubuntu/workspace-7030b71e8f122ce281269f4f67d3843c
Source: https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_110m_admin_0_countries.geojson

# 3. Load the real source data, caching it for repeatable reruns
import requests

if not RAW_PATH.exists():
    response = requests.get(SOURCE_URL, timeout=60)
    response.raise_for_status()
    RAW_PATH.write_bytes(response.content)
    print(f"Downloaded {len(response.content):,} bytes to {RAW_PATH}")
else:
    print(f"Using cached file: {RAW_PATH}")

world = gpd.read_file(RAW_PATH)
# Normalize the current Natural Earth GDP field name for the rest of the notebook.
if "GDP_MD_EST" not in world.columns and "GDP_MD" in world.columns:
    world = world.rename(columns={"GDP_MD": "GDP_MD_EST"})
print(f"Loaded {len(world):,} country/territory features")
print(f"Original CRS: {world.crs}")
print("Columns:", ", ".join(world.columns))
display(world[["NAME", "CONTINENT", "POP_EST", "GDP_MD_EST", "geometry"]].head(3))

Using cached file: /home/ubuntu/workspace-7030b71e8f122ce281269f4f67d3843c/data/ne_110m_admin_0_countries.geojson

Loaded 177 country/territory features
Original CRS: EPSG:4326
Columns: featurecla, scalerank, LABELRANK, SOVEREIGNT, SOV_A3, ADM0_DIF, LEVEL, TYPE, TLC, ADMIN, ADM0_A3, GEOU_DIF, GEOUNIT, GU_A3, SU_DIF, SUBUNIT, SU_A3, BRK_DIFF, NAME, NAME_LONG, BRK_A3, BRK_NAME, BRK_GROUP, ABBREV, POSTAL, FORMAL_EN, FORMAL_FR, NAME_CIAWF, NOTE_ADM0, NOTE_BRK, NAME_SORT, NAME_ALT, MAPCOLOR7, MAPCOLOR8, MAPCOLOR9, MAPCOLOR13, POP_EST, POP_RANK, POP_YEAR, GDP_MD_EST, GDP_YEAR, ECONOMY, INCOME_GRP, FIPS_10, ISO_A2, ISO_A2_EH, ISO_A3, ISO_A3_EH, ISO_N3, ISO_N3_EH, UN_A3, WB_A2, WB_A3, WOE_ID, WOE_ID_EH, WOE_NOTE, ADM0_ISO, ADM0_DIFF, ADM0_TLC, ADM0_A3_US, ADM0_A3_FR, ADM0_A3_RU, ADM0_A3_ES, ADM0_A3_CN, ADM0_A3_TW, ADM0_A3_IN, ADM0_A3_NP, ADM0_A3_PK, ADM0_A3_DE, ADM0_A3_GB, ADM0_A3_BR, ADM0_A3_IL, ADM0_A3_PS, ADM0_A3_SA, ADM0_A3_EG, ADM0_A3_MA, ADM0_A3_PT, ADM0_A3_AR, ADM0_A3_JP, ADM0_A3_KO, ADM0_A3_VN, ADM0_A3_TR, ADM0_A3_ID, ADM0_A3_PL, ADM0_A3_GR, ADM0_A3_IT, ADM0_A3_NL, ADM0_A3_SE, ADM0_A3_BD, ADM0_A3_UA, ADM0_A3_UN, ADM0_A3_WB, CONTINENT, REGION_UN, SUBREGION, REGION_WB, NAME_LEN, LONG_LEN, ABBREV_LEN, TINY, HOMEPART, MIN_ZOOM, MIN_LABEL, MAX_LABEL, LABEL_X, LABEL_Y, NE_ID, WIKIDATAID, NAME_AR, NAME_BN, NAME_DE, NAME_EN, NAME_ES, NAME_FA, NAME_FR, NAME_EL, NAME_HE, NAME_HI, NAME_HU, NAME_ID, NAME_IT, NAME_JA, NAME_KO, NAME_NL, NAME_PL, NAME_PT, NAME_RU, NAME_SV, NAME_TR, NAME_UK, NAME_UR, NAME_VI, NAME_ZH, NAME_ZHT, FCLASS_ISO, TLC_DIFF, FCLASS_TLC, FCLASS_US, FCLASS_FR, FCLASS_RU, FCLASS_ES, FCLASS_CN, FCLASS_TW, FCLASS_IN, FCLASS_NP, FCLASS_PK, FCLASS_DE, FCLASS_GB, FCLASS_BR, FCLASS_IL, FCLASS_PS, FCLASS_SA, FCLASS_EG, FCLASS_MA, FCLASS_PT, FCLASS_AR, FCLASS_JP, FCLASS_KO, FCLASS_VN, FCLASS_TR, FCLASS_ID, FCLASS_PL, FCLASS_GR, FCLASS_IT, FCLASS_NL, FCLASS_SE, FCLASS_BD, FCLASS_UA, geometry

	NAME	CONTINENT	POP_EST	GDP_MD_EST	geometry
0	Fiji	Oceania	889953.0	5496	MULTIPOLYGON (((180 -16.06713, 180 -16.55522, ...
1	Tanzania	Africa	58005463.0	63177	POLYGON ((33.90371 -0.95, 34.07262 -1.05982, 3...
2	W. Sahara	Africa	603253.0	907	POLYGON ((-8.66559 27.65643, -8.66512 27.58948...

Why inspect columns before transforming?

A robust spatial workflow does not assume that a downloaded file has the expected schema. The next checks make the assumptions visible: required columns, duplicate labels, missing geometries, geometry validity, and value ranges.

# 4. Clean and validate the raw data
required = {"NAME", "CONTINENT", "POP_EST", "GDP_MD_EST", "geometry"}
missing_required = required.difference(world.columns)
assert not missing_required, f"Missing required columns: {missing_required}"

validation = pd.Series({
    "rows": len(world),
    "missing_geometries": int(world.geometry.isna().sum()),
    "invalid_geometries": int((~world.geometry.is_valid).sum()),
    "duplicate_country_names": int(world["NAME"].duplicated().sum()),
    "nonpositive_population": int((pd.to_numeric(world["POP_EST"], errors="coerce") <= 0).sum()),
    "nonpositive_gdp_millions": int((pd.to_numeric(world["GDP_MD_EST"], errors="coerce") <= 0).sum()),
})
display(validation.to_frame("value"))

# Keep only records with usable country names and geometries; repair only invalid geometry if needed.
world = world.loc[world["NAME"].notna() & world.geometry.notna()].copy()
if (~world.geometry.is_valid).any():
    world["geometry"] = world.geometry.make_valid()

# Convert source fields explicitly to numeric so bad strings become visible as NaN.
world["POP_EST"] = pd.to_numeric(world["POP_EST"], errors="coerce")
world["GDP_MD_EST"] = pd.to_numeric(world["GDP_MD_EST"], errors="coerce")
print(f"Rows after basic cleaning: {len(world):,}")
assert world.geometry.notna().all()
assert world.geometry.is_valid.all()

	value
rows	177
missing_geometries	0
invalid_geometries	2
duplicate_country_names	0
nonpositive_population	0
nonpositive_gdp_millions	0

Rows after basic cleaning: 177

# 5. Feature engineering in an equal-area CRS
# EPSG:6933 (World Cylindrical Equal Area) avoids using latitude/longitude degrees for area.
world_equal_area = world.to_crs("EPSG:6933")
world["area_km2"] = world_equal_area.geometry.area / 1_000_000
world["gdp_per_person_usd"] = (world["GDP_MD_EST"].astype(float) * 1_000_000) / world["POP_EST"]
world["population_density_per_km2"] = world["POP_EST"] / world["area_km2"]

# Valid modeling set: positive values are required for log transformations.
model_mask = (
    world["POP_EST"].gt(0)
    & world["GDP_MD_EST"].gt(0)
    & world["area_km2"].gt(0)
    & world["gdp_per_person_usd"].gt(0)
    & world["population_density_per_km2"].gt(0)
)
world["model_ready"] = model_mask
world["log_gdp_per_person"] = np.where(model_mask, np.log(world["gdp_per_person_usd"]), np.nan)
world["log_density"] = np.where(model_mask, np.log(world["population_density_per_km2"]), np.nan)

print(f"Model-ready countries/territories: {model_mask.sum():,} of {len(world):,}")
display(world.loc[model_mask, ["NAME", "CONTINENT", "area_km2", "gdp_per_person_usd", "population_density_per_km2"]].describe().T)

Model-ready countries/territories: 177 of 177

	count	mean	std	min	25%	50%	75%	max
area_km2	177.0	832560.843222	2.163569e+06	2415.601364	46171.676608	185256.073256	621832.314600	1.702059e+07
gdp_per_person_usd	177.0	16193.031208	2.567632e+04	261.218430	1816.545232	5789.643726	17828.417039	2.000000e+05
population_density_per_km2	177.0	117.132808	1.617607e+02	0.000364	24.805541	69.290306	124.935390	1.218960e+03

2. Exploratory tables

Before mapping, inspect the extremes. Extreme values are not automatically errors: a small, densely populated territory can be substantively meaningful, while a very large country can have low density. They do, however, influence a regression and should be visible to the reader.

# 6. Useful tables for data understanding
model_df = world.loc[model_mask].copy()
summary_cols = ["NAME", "CONTINENT", "gdp_per_person_usd", "population_density_per_km2", "area_km2"]

print("Five highest estimated GDP per person")
display(model_df.nlargest(5, "gdp_per_person_usd")[summary_cols].round(2))
print("Five highest estimated population density")
display(model_df.nlargest(5, "population_density_per_km2")[summary_cols].round(2))

continent_summary = (
    model_df.groupby("CONTINENT", dropna=False)
    .agg(countries=("NAME", "count"), median_gdp_per_person=("gdp_per_person_usd", "median"),
         median_density=("population_density_per_km2", "median"))
    .sort_values("median_gdp_per_person", ascending=False)
)
display(continent_summary.round(2))

Five highest estimated GDP per person

	NAME	CONTINENT	gdp_per_person_usd	population_density_per_km2	area_km2
159	Antarctica	Antarctica	200000.00	0.00	12337771.52
128	Luxembourg	Europe	114703.11	256.62	2415.60
23	Fr. S. Antarctic Lands	Seven seas (open ocean)	114285.71	0.01	11589.61
20	Falkland Is.	South America	82989.99	0.21	16373.42
127	Switzerland	Europe	81993.68	185.72	46171.68

Five highest estimated population density

	NAME	CONTINENT	gdp_per_person_usd	population_density_per_km2	area_km2
99	Bangladesh	Asia	1855.74	1218.96	133758.41
79	Palestine	Asia	3473.84	930.54	5035.06
140	Taiwan	Asia	47818.31	685.74	34369.12
77	Lebanon	Asia	7583.60	679.58	10088.12
169	Rwanda	Africa	819.99	540.42	23365.25

	countries	median_gdp_per_person	median_density
CONTINENT
Antarctica	1	200000.00	0.00
Seven seas (open ocean)	1	114285.71	0.01
Europe	39	23252.05	90.93
North America	18	9385.74	89.71
South America	13	6977.69	19.57
Oceania	7	6175.61	18.89
Asia	47	4697.64	104.46
Africa	51	1219.37	50.27

3. Interactive geographic exploration

The first map colors countries by log GDP per person. A logarithmic color scale helps prevent a few high-income observations from flattening most of the map. Hover over a country to inspect the underlying values. The map is saved as standalone HTML as well as displayed inline.

# 7. Interactive choropleth: GDP per person
map_gdp = folium.Map(location=[20, 0], zoom_start=2, tiles=None, control_scale=True)

folium.Choropleth(
    geo_data=world.to_json(),
    data=world.loc[model_mask, ["NAME", "log_gdp_per_person"]],
    columns=["NAME", "log_gdp_per_person"],
    key_on="feature.properties.NAME",
    fill_color="YlGnBu",
    nan_fill_color="lightgray",
    line_color="white",
    line_weight=0.3,
    legend_name="log(GDP per person, USD)",
    highlight=True,
).add_to(map_gdp)

folium.GeoJson(
    world.to_json(),
    name="Country details",
    style_function=lambda feature: {"fillOpacity": 0, "weight": 0.4, "color": "#555"},
    tooltip=folium.GeoJsonTooltip(
        fields=["NAME", "CONTINENT", "POP_EST", "GDP_MD_EST"],
        aliases=["Country", "Continent", "Population estimate", "GDP estimate (USD millions)"],
        localize=True,
        sticky=False,
    ),
).add_to(map_gdp)
folium.LayerControl().add_to(map_gdp)
map_gdp.save(OUTPUT_DIR / "interactive_gdp_per_person_map.html")
map_gdp

Make this Notebook Trusted to load map: File -> Trust Notebook

# 8. Interactive choropleth: population density
map_density = folium.Map(location=[20, 0], zoom_start=2, tiles=None, control_scale=True)

folium.Choropleth(
    geo_data=world.to_json(),
    data=world.loc[model_mask, ["NAME", "log_density"]],
    columns=["NAME", "log_density"],
    key_on="feature.properties.NAME",
    fill_color="OrRd",
    nan_fill_color="lightgray",
    line_color="white",
    line_weight=0.3,
    legend_name="log(population density per km²)",
    highlight=True,
).add_to(map_density)
folium.GeoJson(
    world.to_json(),
    name="Country details",
    style_function=lambda feature: {"fillOpacity": 0, "weight": 0.4, "color": "#555"},
    tooltip=folium.GeoJsonTooltip(
        fields=["NAME", "CONTINENT", "population_density_per_km2"],
        aliases=["Country", "Continent", "Population density per km²"],
        localize=True, sticky=False,
    ),
).add_to(map_density)
folium.LayerControl().add_to(map_density)
map_density.save(OUTPUT_DIR / "interactive_population_density_map.html")
map_density

Make this Notebook Trusted to load map: File -> Trust Notebook

How to read the maps

Compare broad spatial patterns rather than interpreting individual colors too literally. The 110m Natural Earth layer is generalized for small-scale mapping, so small islands and narrow borders are simplified. The grey areas are observations that failed the modeling-value checks, not zero values.

4. Descriptive analysis

We fit a simple ordinary least squares model:

[\log( ext{GDP per person}) = eta_0 + eta_1\log( ext{population density}) + \epsilon]

The log-log slope can be read approximately as an elasticity: a 1% higher population density is associated with about (eta_1)% higher estimated GDP per person in this cross-sectional dataset. This is descriptive and does not establish causality.

# 9. Log-log regression and visual diagnostic
X = model_df[["log_density"]].to_numpy()
y = model_df["log_gdp_per_person"].to_numpy()
reg = LinearRegression().fit(X, y)
pred = reg.predict(X)

model_df["predicted_log_gdp_per_person"] = pred
model_df["residual"] = y - pred
metrics = pd.Series({
    "observations": len(model_df),
    "intercept": reg.intercept_,
    "slope_on_log_density": reg.coef_[0],
    "R_squared": r2_score(y, pred),
    "MAE_in_log_units": mean_absolute_error(y, pred),
    "excluded_from_model": int((~model_mask).sum()),
})
display(metrics.to_frame("value").round(4))

fig, ax = plt.subplots(figsize=(9, 6))
sns.scatterplot(data=model_df, x="log_density", y="log_gdp_per_person", hue="CONTINENT",
                palette="tab10", alpha=0.8, ax=ax)
order = np.argsort(model_df["log_density"].to_numpy())
ax.plot(model_df["log_density"].to_numpy()[order], pred[order], color="black", linewidth=2, label="OLS fit")
ax.set(title="Estimated GDP per person vs. population density",
       xlabel="log(population density per km²)", ylabel="log(GDP per person, USD)")
ax.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "gdp_density_loglog_relationship.png", dpi=160, bbox_inches="tight")
plt.show()

	value
observations	177.0000
intercept	9.3061
slope_on_log_density	-0.1586
R_squared	0.0370
MAE_in_log_units	1.2027
excluded_from_model	0.0000

# 10. Edge-case checks: influential extremes and residuals
# These checks do not remove observations automatically; they help the reader see sensitivity.
print("Largest absolute residuals")
display(model_df.assign(abs_residual=model_df["residual"].abs())
        .nlargest(8, "abs_residual")[["NAME", "CONTINENT", "log_density", "log_gdp_per_person", "residual"]]
        .round(3))

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
sns.histplot(model_df["residual"], kde=True, ax=axes[0], color="#4C78A8")
axes[0].axvline(0, color="black", linewidth=1)
axes[0].set_title("Residual distribution")
axes[0].set_xlabel("Observed - predicted log GDP per person")
sns.scatterplot(data=model_df, x="log_density", y="residual", hue="CONTINENT",
                palette="tab10", alpha=0.8, ax=axes[1], legend=False)
axes[1].axhline(0, color="black", linewidth=1)
axes[1].set_title("Residuals vs. log density")
axes[1].set_xlabel("log(population density per km²)")
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "regression_residual_checks.png", dpi=160, bbox_inches="tight")
plt.show()

Largest absolute residuals

	NAME	CONTINENT	log_density	log_gdp_per_person	residual
128	Luxembourg	Europe	5.548	11.650	3.224
154	Eritrea	Africa	3.932	5.828	-2.855
127	Switzerland	Europe	5.224	11.314	2.837
66	Central African Rep.	Africa	2.032	6.148	-2.836
75	Burundi	Africa	6.086	5.565	-2.775
12	Somalia	Africa	3.047	6.138	-2.685
133	Ireland	Europe	4.437	11.273	2.671
85	Qatar	Asia	5.521	11.036	2.606

# 11. A compact, reproducible result statement
slope = reg.coef_[0]
r2 = r2_score(y, pred)
comparison = "positive" if slope > 0 else "negative"
print(
    f"Among {len(model_df)} model-ready countries/territories, the fitted association is {comparison}: "
    f"a 1% increase in population density corresponds to about {slope:.2f}% change in estimated GDP per person "
    f"on average in this log-log specification. The model R² is {r2:.2f}."
)

Among 177 model-ready countries/territories, the fitted association is negative: a 1% increase in population density corresponds to about -0.16% change in estimated GDP per person on average in this log-log specification. The model R² is 0.04.

5. Summary, limitations, and next steps

Summary

• The workflow loads and caches a real public polygon dataset, validates its schema and geometries, calculates area in an equal-area CRS, and makes missing/non-positive values explicit.
• Two interactive maps show the geographic distribution of estimated GDP per person and population density.
• A log-log regression provides a compact descriptive summary of the cross-country relationship; the notebook reports the number of modeled and excluded observations and includes residual checks.

Limitations

• Natural Earth's attributes are estimates and may be dated relative to the execution date.
• Country-level averages conceal within-country inequality and local variation (ecological fallacy).
• Polygon area and population density are sensitive to boundary definitions, disputed territories, and the generalized 1:110m geometry.
• The regression is observational and omits confounders such as urbanization, institutions, geography, and measurement quality; it should not be interpreted causally.
• GDP per person is derived from rounded GDP/population estimates and is not PPP-adjusted.

Next steps

1. Replace Natural Earth estimates with a versioned statistical source such as World Bank indicators for a specified year.
2. Use higher-resolution administrative boundaries for within-country analysis.
3. Add spatial autocorrelation diagnostics and a spatial regression if residuals cluster geographically.
4. Test sensitivity to excluding microstates, using alternative density definitions, and changing the year/source.
5. Publish the exact raw-data checksum and environment lockfile for stricter reproducibility.

Reproducibility outputs

When executed, the notebook writes the following files under outputs/:

• interactive_gdp_per_person_map.html
• interactive_population_density_map.html
• gdp_density_loglog_relationship.png
• regression_residual_checks.png

The raw source is cached under data/ne_110m_admin_0_countries.geojson.

6. Advanced spatial analysis: global autocorrelation and local clusters

A choropleth can make nearby countries look similar, but visual similarity is not the same as spatial dependence. To test whether high or low GDP-per-person values cluster geographically, we use two complementary statistics:

• Global Moran’s I asks whether the whole map is more spatially patterned than a random arrangement.
• Local Moran’s I identifies country-level spatial associations. A significant High–High result is a hot-spot-like cluster of a high-value country surrounded by high-value neighbors; Low–Low is the corresponding cold spot. High–Low and Low–High are spatial outliers.

Because country polygons include islands and disconnected territories, we use a reproducible 4-nearest-neighbor graph based on geographic centroids. This guarantees every analyzed observation has neighbors, but it is a modeling choice rather than a natural law. Significance is assessed with 999 random permutations and a fixed seed of 42.

# 12. Global Moran's I and local cluster detection
from libpysal.weights import KNN
from esda.moran import Moran, Moran_Local
from branca.element import Element

cluster_df = world.loc[model_mask].copy()
cluster_wgs84 = cluster_df.to_crs("EPSG:4326")
neighbor_k = min(4, len(cluster_df) - 1)
weights = KNN.from_dataframe(cluster_wgs84, k=neighbor_k)
weights.transform = "R"

np.random.seed(42)
y_cluster = cluster_df["log_gdp_per_person"].to_numpy()
global_moran = Moran(y_cluster, weights, permutations=999)
local_moran = Moran_Local(y_cluster, weights, permutations=999, seed=42)

quadrant_labels = {1: "High–High", 2: "Low–High", 3: "Low–Low", 4: "High–Low"}
cluster_df["local_moran_I"] = local_moran.Is
cluster_df["local_p_value"] = local_moran.p_sim
cluster_df["cluster_type"] = [
    quadrant_labels[q] if p < 0.05 else "Not significant"
    for q, p in zip(local_moran.q, local_moran.p_sim)
]

cluster_counts = cluster_df["cluster_type"].value_counts().reindex(
    ["High–High", "Low–Low", "High–Low", "Low–High", "Not significant"], fill_value=0
)
cluster_summary = pd.DataFrame({"cluster_type": cluster_counts.index, "countries": cluster_counts.values})
cluster_summary.to_csv(OUTPUT_DIR / "spatial_cluster_summary.csv", index=False)

print(f"Global Moran's I: {global_moran.I:.3f}")
print(f"Permutation p-value: {global_moran.p_sim:.3f}")
print(f"Neighbors per country: {neighbor_k}")
display(cluster_summary)
print("Significant local clusters")
display(cluster_df.loc[cluster_df["cluster_type"] != "Not significant", ["NAME", "cluster_type", "local_p_value", "local_moran_I"]].sort_values(["cluster_type", "NAME"]).head(30).round(4))

Global Moran's I: 0.609
Permutation p-value: 0.001
Neighbors per country: 4

/usr/local/lib/python3.12/dist-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

	cluster_type	countries
0	High–High	19
1	Low–Low	44
2	High–Low	3
3	Low–High	2
4	Not significant	109

Significant local clusters

	NAME	cluster_type	local_p_value	local_moran_I
129	Belgium	High–High	0.001	2.2730
3	Canada	High–High	0.041	1.1746
153	Czechia	High–High	0.026	0.8929
142	Denmark	High–High	0.001	2.4112
151	Finland	High–High	0.024	1.3860
43	France	High–High	0.006	1.5120
121	Germany	High–High	0.001	2.3334
22	Greenland	High–High	0.001	2.0762
144	Iceland	High–High	0.001	2.4366
133	Ireland	High–High	0.009	2.1250
141	Italy	High–High	0.008	1.3883
128	Luxembourg	High–High	0.001	2.9702
130	Netherlands	High–High	0.001	2.3002
21	Norway	High–High	0.002	2.2801
131	Portugal	High–High	0.036	0.8494
85	Qatar	High–High	0.024	1.5801
110	Sweden	High–High	0.002	2.0325
127	Switzerland	High–High	0.001	2.5989
143	United Kingdom	High–High	0.001	2.1619
23	Fr. S. Antarctic Lands	High–Low	0.002	-2.3171
164	Libya	High–Low	0.023	-0.1685
106	Turkmenistan	High–Low	0.021	-0.0943
107	Iran	Low–High	0.019	-0.0537
95	North Korea	Low–High	0.022	-0.9342
82	Algeria	Low–Low	0.025	0.2432
74	Angola	Low–Low	0.039	0.4186
99	Bangladesh	Low–Low	0.047	0.6194
54	Benin	Low–Low	0.008	1.0895
100	Bhutan	Low–Low	0.027	0.3472
65	Burkina Faso	Low–Low	0.004	1.5032

# 13. Interactive local cluster map
cluster_map = cluster_df.to_crs("EPSG:4326").copy()
cluster_colors = {
    "High–High": "#ef7254",
    "Low–Low": "#3b6ea8",
    "High–Low": "#f4b942",
    "Low–High": "#8e6bbf",
    "Not significant": "#c8cdd1",
}

def cluster_style(feature):
    label = feature["properties"].get("cluster_type", "Not significant")
    color = cluster_colors.get(label, "#c8cdd1")
    return {"fillColor": color, "color": "#ffffff", "weight": 0.5, "fillOpacity": 0.82}

hotspot_map = folium.Map(location=[20, 0], zoom_start=2, tiles=None, control_scale=True)
folium.GeoJson(
    cluster_map.to_json(),
    name="Local Moran clusters",
    style_function=cluster_style,
    highlight_function=lambda feature: {"weight": 2, "color": "#102a43", "fillOpacity": 0.95},
    tooltip=folium.GeoJsonTooltip(
        fields=["NAME", "CONTINENT", "cluster_type", "local_p_value", "local_moran_I"],
        aliases=["Country", "Continent", "Cluster type", "Permutation p-value", "Local Moran's I"],
        localize=True, sticky=False,
    ),
).add_to(hotspot_map)
folium.LayerControl(collapsed=True).add_to(hotspot_map)
legend_html = """<div style='position: fixed; bottom: 28px; left: 28px; z-index: 9999; background: white; padding: 12px 14px; border: 1px solid #9aa5ad; font: 12px Arial; color: #102a43'><b>Local Moran cluster</b><br>{items}</div>""".format(items="".join(f"<div><span style='display:inline-block;width:11px;height:11px;background:{color};margin-right:7px'></span>{label}</div>" for label, color in cluster_colors.items()))
hotspot_map.get_root().html.add_child(Element(legend_html))
hotspot_map.save(PROJECT_ROOT / "maps" / "interactive_spatial_clusters.html")
hotspot_map

Make this Notebook Trusted to load map: File -> Trust Notebook

Interpreting the cluster map

Clusters are local spatial associations under the chosen 4-nearest-neighbor graph and 0.05 permutation threshold. They are useful for asking better follow-up questions, not for declaring that geography causes income. A country can be a high-value spatial outlier even when its neighbors are relatively low, and a grey country may simply lack enough local evidence under this specification.