{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e8d2e41c",
   "metadata": {},
   "source": [
    "# Interactive Geographic Analysis: Population Density and GDP per Person\n",
    "\n",
    "**Purpose.** This notebook is a reproducible teaching example for mapping country-level geographic data. It asks:\n",
    "\n",
    "> **How does estimated GDP per person relate to population density across countries, and what spatial patterns are visible on a map?**\n",
    "\n",
    "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.\n",
    "\n",
    "## What we use\n",
    "\n",
    "- **Natural Earth Admin 0 countries, 1:110m** for country polygons and published estimates of population and GDP.\n",
    "- A projected, equal-area coordinate reference system (**EPSG:6933**) for area and density calculations.\n",
    "- `geopandas` for spatial data, `folium` for interactive maps, and `scikit-learn` for a transparent log-log regression.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba8d9509",
   "metadata": {},
   "source": [
    "## 1. Assumptions and analysis choices\n",
    "\n",
    "1. A country polygon is treated as the unit of analysis; this is an ecological, country-level comparison, not an individual-level study.\n",
    "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.\n",
    "3. Population density is computed as estimated population divided by polygon area in square kilometres.\n",
    "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.\n",
    "5. The map is for exploration. It is not evidence of causal effects, and country boundaries and estimates have uncertainty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3503e6f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. Reproducible setup\n",
    "from pathlib import Path\n",
    "import json\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import folium\n",
    "from branca.colormap import linear\n",
    "from IPython.display import display, Markdown\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import r2_score, mean_absolute_error\n",
    "\n",
    "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
    "sns.set_theme(style=\"whitegrid\", context=\"notebook\")\n",
    "\n",
    "PROJECT_ROOT = Path.cwd()\n",
    "DATA_DIR = PROJECT_ROOT / \"data\"\n",
    "OUTPUT_DIR = PROJECT_ROOT / \"outputs\"\n",
    "DATA_DIR.mkdir(exist_ok=True)\n",
    "OUTPUT_DIR.mkdir(exist_ok=True)\n",
    "\n",
    "SOURCE_URL = (\n",
    "    \"https://raw.githubusercontent.com/nvkelso/natural-earth-vector/\"\n",
    "    \"master/geojson/ne_110m_admin_0_countries.geojson\"\n",
    ")\n",
    "RAW_PATH = DATA_DIR / \"ne_110m_admin_0_countries.geojson\"\n",
    "print(f\"Project root: {PROJECT_ROOT.resolve()}\")\n",
    "print(f\"Source: {SOURCE_URL}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64646376",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3. Load the real source data, caching it for repeatable reruns\n",
    "import requests\n",
    "\n",
    "if not RAW_PATH.exists():\n",
    "    response = requests.get(SOURCE_URL, timeout=60)\n",
    "    response.raise_for_status()\n",
    "    RAW_PATH.write_bytes(response.content)\n",
    "    print(f\"Downloaded {len(response.content):,} bytes to {RAW_PATH}\")\n",
    "else:\n",
    "    print(f\"Using cached file: {RAW_PATH}\")\n",
    "\n",
    "world = gpd.read_file(RAW_PATH)\n",
    "# Normalize the current Natural Earth GDP field name for the rest of the notebook.\n",
    "if \"GDP_MD_EST\" not in world.columns and \"GDP_MD\" in world.columns:\n",
    "    world = world.rename(columns={\"GDP_MD\": \"GDP_MD_EST\"})\n",
    "print(f\"Loaded {len(world):,} country/territory features\")\n",
    "print(f\"Original CRS: {world.crs}\")\n",
    "print(\"Columns:\", \", \".join(world.columns))\n",
    "display(world[[\"NAME\", \"CONTINENT\", \"POP_EST\", \"GDP_MD_EST\", \"geometry\"]].head(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cdd86b9",
   "metadata": {},
   "source": [
    "### Why inspect columns before transforming?\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4f0eb36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4. Clean and validate the raw data\n",
    "required = {\"NAME\", \"CONTINENT\", \"POP_EST\", \"GDP_MD_EST\", \"geometry\"}\n",
    "missing_required = required.difference(world.columns)\n",
    "assert not missing_required, f\"Missing required columns: {missing_required}\"\n",
    "\n",
    "validation = pd.Series({\n",
    "    \"rows\": len(world),\n",
    "    \"missing_geometries\": int(world.geometry.isna().sum()),\n",
    "    \"invalid_geometries\": int((~world.geometry.is_valid).sum()),\n",
    "    \"duplicate_country_names\": int(world[\"NAME\"].duplicated().sum()),\n",
    "    \"nonpositive_population\": int((pd.to_numeric(world[\"POP_EST\"], errors=\"coerce\") <= 0).sum()),\n",
    "    \"nonpositive_gdp_millions\": int((pd.to_numeric(world[\"GDP_MD_EST\"], errors=\"coerce\") <= 0).sum()),\n",
    "})\n",
    "display(validation.to_frame(\"value\"))\n",
    "\n",
    "# Keep only records with usable country names and geometries; repair only invalid geometry if needed.\n",
    "world = world.loc[world[\"NAME\"].notna() & world.geometry.notna()].copy()\n",
    "if (~world.geometry.is_valid).any():\n",
    "    world[\"geometry\"] = world.geometry.make_valid()\n",
    "\n",
    "# Convert source fields explicitly to numeric so bad strings become visible as NaN.\n",
    "world[\"POP_EST\"] = pd.to_numeric(world[\"POP_EST\"], errors=\"coerce\")\n",
    "world[\"GDP_MD_EST\"] = pd.to_numeric(world[\"GDP_MD_EST\"], errors=\"coerce\")\n",
    "print(f\"Rows after basic cleaning: {len(world):,}\")\n",
    "assert world.geometry.notna().all()\n",
    "assert world.geometry.is_valid.all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9378743",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5. Feature engineering in an equal-area CRS\n",
    "# EPSG:6933 (World Cylindrical Equal Area) avoids using latitude/longitude degrees for area.\n",
    "world_equal_area = world.to_crs(\"EPSG:6933\")\n",
    "world[\"area_km2\"] = world_equal_area.geometry.area / 1_000_000\n",
    "world[\"gdp_per_person_usd\"] = (world[\"GDP_MD_EST\"].astype(float) * 1_000_000) / world[\"POP_EST\"]\n",
    "world[\"population_density_per_km2\"] = world[\"POP_EST\"] / world[\"area_km2\"]\n",
    "\n",
    "# Valid modeling set: positive values are required for log transformations.\n",
    "model_mask = (\n",
    "    world[\"POP_EST\"].gt(0)\n",
    "    & world[\"GDP_MD_EST\"].gt(0)\n",
    "    & world[\"area_km2\"].gt(0)\n",
    "    & world[\"gdp_per_person_usd\"].gt(0)\n",
    "    & world[\"population_density_per_km2\"].gt(0)\n",
    ")\n",
    "world[\"model_ready\"] = model_mask\n",
    "world[\"log_gdp_per_person\"] = np.where(model_mask, np.log(world[\"gdp_per_person_usd\"]), np.nan)\n",
    "world[\"log_density\"] = np.where(model_mask, np.log(world[\"population_density_per_km2\"]), np.nan)\n",
    "\n",
    "print(f\"Model-ready countries/territories: {model_mask.sum():,} of {len(world):,}\")\n",
    "display(world.loc[model_mask, [\"NAME\", \"CONTINENT\", \"area_km2\", \"gdp_per_person_usd\", \"population_density_per_km2\"]].describe().T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ebc3b01",
   "metadata": {},
   "source": [
    "## 2. Exploratory tables\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "965ad901",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 6. Useful tables for data understanding\n",
    "model_df = world.loc[model_mask].copy()\n",
    "summary_cols = [\"NAME\", \"CONTINENT\", \"gdp_per_person_usd\", \"population_density_per_km2\", \"area_km2\"]\n",
    "\n",
    "print(\"Five highest estimated GDP per person\")\n",
    "display(model_df.nlargest(5, \"gdp_per_person_usd\")[summary_cols].round(2))\n",
    "print(\"Five highest estimated population density\")\n",
    "display(model_df.nlargest(5, \"population_density_per_km2\")[summary_cols].round(2))\n",
    "\n",
    "continent_summary = (\n",
    "    model_df.groupby(\"CONTINENT\", dropna=False)\n",
    "    .agg(countries=(\"NAME\", \"count\"), median_gdp_per_person=(\"gdp_per_person_usd\", \"median\"),\n",
    "         median_density=(\"population_density_per_km2\", \"median\"))\n",
    "    .sort_values(\"median_gdp_per_person\", ascending=False)\n",
    ")\n",
    "display(continent_summary.round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9551842a",
   "metadata": {},
   "source": [
    "## 3. Interactive geographic exploration\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81e78b16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 7. Interactive choropleth: GDP per person\n",
    "map_gdp = folium.Map(location=[20, 0], zoom_start=2, tiles=None, control_scale=True)\n",
    "\n",
    "folium.Choropleth(\n",
    "    geo_data=world.to_json(),\n",
    "    data=world.loc[model_mask, [\"NAME\", \"log_gdp_per_person\"]],\n",
    "    columns=[\"NAME\", \"log_gdp_per_person\"],\n",
    "    key_on=\"feature.properties.NAME\",\n",
    "    fill_color=\"YlGnBu\",\n",
    "    nan_fill_color=\"lightgray\",\n",
    "    line_color=\"white\",\n",
    "    line_weight=0.3,\n",
    "    legend_name=\"log(GDP per person, USD)\",\n",
    "    highlight=True,\n",
    ").add_to(map_gdp)\n",
    "\n",
    "folium.GeoJson(\n",
    "    world.to_json(),\n",
    "    name=\"Country details\",\n",
    "    style_function=lambda feature: {\"fillOpacity\": 0, \"weight\": 0.4, \"color\": \"#555\"},\n",
    "    tooltip=folium.GeoJsonTooltip(\n",
    "        fields=[\"NAME\", \"CONTINENT\", \"POP_EST\", \"GDP_MD_EST\"],\n",
    "        aliases=[\"Country\", \"Continent\", \"Population estimate\", \"GDP estimate (USD millions)\"],\n",
    "        localize=True,\n",
    "        sticky=False,\n",
    "    ),\n",
    ").add_to(map_gdp)\n",
    "folium.LayerControl().add_to(map_gdp)\n",
    "map_gdp.save(OUTPUT_DIR / \"interactive_gdp_per_person_map.html\")\n",
    "map_gdp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79537fae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 8. Interactive choropleth: population density\n",
    "map_density = folium.Map(location=[20, 0], zoom_start=2, tiles=None, control_scale=True)\n",
    "\n",
    "folium.Choropleth(\n",
    "    geo_data=world.to_json(),\n",
    "    data=world.loc[model_mask, [\"NAME\", \"log_density\"]],\n",
    "    columns=[\"NAME\", \"log_density\"],\n",
    "    key_on=\"feature.properties.NAME\",\n",
    "    fill_color=\"OrRd\",\n",
    "    nan_fill_color=\"lightgray\",\n",
    "    line_color=\"white\",\n",
    "    line_weight=0.3,\n",
    "    legend_name=\"log(population density per km²)\",\n",
    "    highlight=True,\n",
    ").add_to(map_density)\n",
    "folium.GeoJson(\n",
    "    world.to_json(),\n",
    "    name=\"Country details\",\n",
    "    style_function=lambda feature: {\"fillOpacity\": 0, \"weight\": 0.4, \"color\": \"#555\"},\n",
    "    tooltip=folium.GeoJsonTooltip(\n",
    "        fields=[\"NAME\", \"CONTINENT\", \"population_density_per_km2\"],\n",
    "        aliases=[\"Country\", \"Continent\", \"Population density per km²\"],\n",
    "        localize=True, sticky=False,\n",
    "    ),\n",
    ").add_to(map_density)\n",
    "folium.LayerControl().add_to(map_density)\n",
    "map_density.save(OUTPUT_DIR / \"interactive_population_density_map.html\")\n",
    "map_density"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60017c06",
   "metadata": {},
   "source": [
    "### How to read the maps\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5580c6f1",
   "metadata": {},
   "source": [
    "## 4. Descriptive analysis\n",
    "\n",
    "We fit a simple ordinary least squares model:\n",
    "\n",
    "\\[\\log(\text{GDP per person}) = \beta_0 + \beta_1\\log(\text{population density}) + \\epsilon\\]\n",
    "\n",
    "The log-log slope can be read approximately as an elasticity: a 1% higher population density is associated with about \\(\beta_1\\)% higher estimated GDP per person in this cross-sectional dataset. This is descriptive and does not establish causality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bff8496",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9. Log-log regression and visual diagnostic\n",
    "X = model_df[[\"log_density\"]].to_numpy()\n",
    "y = model_df[\"log_gdp_per_person\"].to_numpy()\n",
    "reg = LinearRegression().fit(X, y)\n",
    "pred = reg.predict(X)\n",
    "\n",
    "model_df[\"predicted_log_gdp_per_person\"] = pred\n",
    "model_df[\"residual\"] = y - pred\n",
    "metrics = pd.Series({\n",
    "    \"observations\": len(model_df),\n",
    "    \"intercept\": reg.intercept_,\n",
    "    \"slope_on_log_density\": reg.coef_[0],\n",
    "    \"R_squared\": r2_score(y, pred),\n",
    "    \"MAE_in_log_units\": mean_absolute_error(y, pred),\n",
    "    \"excluded_from_model\": int((~model_mask).sum()),\n",
    "})\n",
    "display(metrics.to_frame(\"value\").round(4))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 6))\n",
    "sns.scatterplot(data=model_df, x=\"log_density\", y=\"log_gdp_per_person\", hue=\"CONTINENT\",\n",
    "                palette=\"tab10\", alpha=0.8, ax=ax)\n",
    "order = np.argsort(model_df[\"log_density\"].to_numpy())\n",
    "ax.plot(model_df[\"log_density\"].to_numpy()[order], pred[order], color=\"black\", linewidth=2, label=\"OLS fit\")\n",
    "ax.set(title=\"Estimated GDP per person vs. population density\",\n",
    "       xlabel=\"log(population density per km²)\", ylabel=\"log(GDP per person, USD)\")\n",
    "ax.legend(bbox_to_anchor=(1.02, 1), loc=\"upper left\")\n",
    "fig.tight_layout()\n",
    "fig.savefig(OUTPUT_DIR / \"gdp_density_loglog_relationship.png\", dpi=160, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9d82926",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 10. Edge-case checks: influential extremes and residuals\n",
    "# These checks do not remove observations automatically; they help the reader see sensitivity.\n",
    "print(\"Largest absolute residuals\")\n",
    "display(model_df.assign(abs_residual=model_df[\"residual\"].abs())\n",
    "        .nlargest(8, \"abs_residual\")[[\"NAME\", \"CONTINENT\", \"log_density\", \"log_gdp_per_person\", \"residual\"]]\n",
    "        .round(3))\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))\n",
    "sns.histplot(model_df[\"residual\"], kde=True, ax=axes[0], color=\"#4C78A8\")\n",
    "axes[0].axvline(0, color=\"black\", linewidth=1)\n",
    "axes[0].set_title(\"Residual distribution\")\n",
    "axes[0].set_xlabel(\"Observed - predicted log GDP per person\")\n",
    "sns.scatterplot(data=model_df, x=\"log_density\", y=\"residual\", hue=\"CONTINENT\",\n",
    "                palette=\"tab10\", alpha=0.8, ax=axes[1], legend=False)\n",
    "axes[1].axhline(0, color=\"black\", linewidth=1)\n",
    "axes[1].set_title(\"Residuals vs. log density\")\n",
    "axes[1].set_xlabel(\"log(population density per km²)\")\n",
    "fig.tight_layout()\n",
    "fig.savefig(OUTPUT_DIR / \"regression_residual_checks.png\", dpi=160, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb6d9124",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 11. A compact, reproducible result statement\n",
    "slope = reg.coef_[0]\n",
    "r2 = r2_score(y, pred)\n",
    "comparison = \"positive\" if slope > 0 else \"negative\"\n",
    "print(\n",
    "    f\"Among {len(model_df)} model-ready countries/territories, the fitted association is {comparison}: \"\n",
    "    f\"a 1% increase in population density corresponds to about {slope:.2f}% change in estimated GDP per person \"\n",
    "    f\"on average in this log-log specification. The model R² is {r2:.2f}.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ce3fc9d",
   "metadata": {},
   "source": [
    "## 5. Summary, limitations, and next steps\n",
    "\n",
    "### Summary\n",
    "\n",
    "- 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.\n",
    "- Two interactive maps show the geographic distribution of estimated GDP per person and population density.\n",
    "- 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.\n",
    "\n",
    "### Limitations\n",
    "\n",
    "- Natural Earth's attributes are estimates and may be dated relative to the execution date.\n",
    "- Country-level averages conceal within-country inequality and local variation (ecological fallacy).\n",
    "- Polygon area and population density are sensitive to boundary definitions, disputed territories, and the generalized 1:110m geometry.\n",
    "- The regression is observational and omits confounders such as urbanization, institutions, geography, and measurement quality; it should not be interpreted causally.\n",
    "- GDP per person is derived from rounded GDP/population estimates and is not PPP-adjusted.\n",
    "\n",
    "### Next steps\n",
    "\n",
    "1. Replace Natural Earth estimates with a versioned statistical source such as World Bank indicators for a specified year.\n",
    "2. Use higher-resolution administrative boundaries for within-country analysis.\n",
    "3. Add spatial autocorrelation diagnostics and a spatial regression if residuals cluster geographically.\n",
    "4. Test sensitivity to excluding microstates, using alternative density definitions, and changing the year/source.\n",
    "5. Publish the exact raw-data checksum and environment lockfile for stricter reproducibility."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e8ad258",
   "metadata": {},
   "source": [
    "## Reproducibility outputs\n",
    "\n",
    "When executed, the notebook writes the following files under `outputs/`:\n",
    "\n",
    "- `interactive_gdp_per_person_map.html`\n",
    "- `interactive_population_density_map.html`\n",
    "- `gdp_density_loglog_relationship.png`\n",
    "- `regression_residual_checks.png`\n",
    "\n",
    "The raw source is cached under `data/ne_110m_admin_0_countries.geojson`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69f43cc7",
   "metadata": {},
   "source": [
    "## 6. Advanced spatial analysis: global autocorrelation and local clusters\n",
    "\n",
    "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:\n",
    "\n",
    "- **Global Moran’s I** asks whether the whole map is more spatially patterned than a random arrangement.\n",
    "- **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.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21dca894",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 12. Global Moran's I and local cluster detection\n",
    "from libpysal.weights import KNN\n",
    "from esda.moran import Moran, Moran_Local\n",
    "from branca.element import Element\n",
    "\n",
    "cluster_df = world.loc[model_mask].copy()\n",
    "cluster_wgs84 = cluster_df.to_crs(\"EPSG:4326\")\n",
    "neighbor_k = min(4, len(cluster_df) - 1)\n",
    "weights = KNN.from_dataframe(cluster_wgs84, k=neighbor_k)\n",
    "weights.transform = \"R\"\n",
    "\n",
    "np.random.seed(42)\n",
    "y_cluster = cluster_df[\"log_gdp_per_person\"].to_numpy()\n",
    "global_moran = Moran(y_cluster, weights, permutations=999)\n",
    "local_moran = Moran_Local(y_cluster, weights, permutations=999, seed=42)\n",
    "\n",
    "quadrant_labels = {1: \"High–High\", 2: \"Low–High\", 3: \"Low–Low\", 4: \"High–Low\"}\n",
    "cluster_df[\"local_moran_I\"] = local_moran.Is\n",
    "cluster_df[\"local_p_value\"] = local_moran.p_sim\n",
    "cluster_df[\"cluster_type\"] = [\n",
    "    quadrant_labels[q] if p < 0.05 else \"Not significant\"\n",
    "    for q, p in zip(local_moran.q, local_moran.p_sim)\n",
    "]\n",
    "\n",
    "cluster_counts = cluster_df[\"cluster_type\"].value_counts().reindex(\n",
    "    [\"High–High\", \"Low–Low\", \"High–Low\", \"Low–High\", \"Not significant\"], fill_value=0\n",
    ")\n",
    "cluster_summary = pd.DataFrame({\"cluster_type\": cluster_counts.index, \"countries\": cluster_counts.values})\n",
    "cluster_summary.to_csv(OUTPUT_DIR / \"spatial_cluster_summary.csv\", index=False)\n",
    "\n",
    "print(f\"Global Moran's I: {global_moran.I:.3f}\")\n",
    "print(f\"Permutation p-value: {global_moran.p_sim:.3f}\")\n",
    "print(f\"Neighbors per country: {neighbor_k}\")\n",
    "display(cluster_summary)\n",
    "print(\"Significant local clusters\")\n",
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eedcfcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 13. Interactive local cluster map\n",
    "cluster_map = cluster_df.to_crs(\"EPSG:4326\").copy()\n",
    "cluster_colors = {\n",
    "    \"High–High\": \"#ef7254\",\n",
    "    \"Low–Low\": \"#3b6ea8\",\n",
    "    \"High–Low\": \"#f4b942\",\n",
    "    \"Low–High\": \"#8e6bbf\",\n",
    "    \"Not significant\": \"#c8cdd1\",\n",
    "}\n",
    "\n",
    "def cluster_style(feature):\n",
    "    label = feature[\"properties\"].get(\"cluster_type\", \"Not significant\")\n",
    "    color = cluster_colors.get(label, \"#c8cdd1\")\n",
    "    return {\"fillColor\": color, \"color\": \"#ffffff\", \"weight\": 0.5, \"fillOpacity\": 0.82}\n",
    "\n",
    "hotspot_map = folium.Map(location=[20, 0], zoom_start=2, tiles=None, control_scale=True)\n",
    "folium.GeoJson(\n",
    "    cluster_map.to_json(),\n",
    "    name=\"Local Moran clusters\",\n",
    "    style_function=cluster_style,\n",
    "    highlight_function=lambda feature: {\"weight\": 2, \"color\": \"#102a43\", \"fillOpacity\": 0.95},\n",
    "    tooltip=folium.GeoJsonTooltip(\n",
    "        fields=[\"NAME\", \"CONTINENT\", \"cluster_type\", \"local_p_value\", \"local_moran_I\"],\n",
    "        aliases=[\"Country\", \"Continent\", \"Cluster type\", \"Permutation p-value\", \"Local Moran's I\"],\n",
    "        localize=True, sticky=False,\n",
    "    ),\n",
    ").add_to(hotspot_map)\n",
    "folium.LayerControl(collapsed=True).add_to(hotspot_map)\n",
    "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()))\n",
    "hotspot_map.get_root().html.add_child(Element(legend_html))\n",
    "hotspot_map.save(PROJECT_ROOT / \"maps\" / \"interactive_spatial_clusters.html\")\n",
    "hotspot_map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "622caf67",
   "metadata": {},
   "source": [
    "### Interpreting the cluster map\n",
    "\n",
    "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."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}