Source code: Interactive web maps from Python with Folium and GeoPandas .explore()

"""Interactive web maps with Folium - TechnicalGIS tutorial.

Builds two Leaflet maps from Python:
  1. a one-line choropleth with GeoDataFrame.explore()
  2. a hand-crafted Folium map with a choropleth, clustered city markers,
     multiple basemaps, plugins and a layer switcher.
Run with:  python folium_maps.py   (then open output/*.html in a browser)
"""
import urllib.request
from pathlib import Path

import branca.colormap as cm
import folium
import geopandas as gpd
from folium import plugins

HERE = Path(__file__).parent
DATA = HERE / "data"
OUT = HERE / "output"
DATA.mkdir(exist_ok=True)
OUT.mkdir(exist_ok=True)
NE = "https://naciscdn.org/naturalearth/50m/cultural/"
# Keyless basemaps. (CARTO's free tiles now require an API key for most uses.)
ESRI_GRAY = "https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}"
ESRI_GRAY_ATTR = "Tiles © Esri — Esri, DeLorme, NAVTEQ"


def fetch(name: str) -> Path:
    dest = DATA / name
    if not dest.exists():
        req = urllib.request.Request(NE + name, headers={"User-Agent": "technicalgis-tutorial"})
        with urllib.request.urlopen(req) as r:
            dest.write_bytes(r.read())
    return dest


def section(title: str) -> None:
    print(f"\n=== {title} ===")


# --- Step 1: prepare the data ---------------------------------------------------
countries = gpd.read_file(fetch("ne_50m_admin_0_countries.zip"))
countries = countries[["NAME", "CONTINENT", "POP_EST", "GDP_MD", "geometry"]].rename(
    columns={"NAME": "name", "CONTINENT": "continent", "POP_EST": "population", "GDP_MD": "gdp_musd"})
europe = countries[countries["continent"] == "Europe"].copy()
europe["area_km2"] = europe.to_crs("EPSG:3035").area / 1e6       # ETRS89 LAEA Europe, equal area
europe["density"] = (europe["population"] / europe["area_km2"]).round(1)
europe["gdp_per_capita"] = (europe["gdp_musd"] * 1e6 / europe["population"]).round(0)
europe = europe[europe["population"] > 0]

places = gpd.read_file(fetch("ne_50m_populated_places.zip"))[["NAME", "ADM0NAME", "POP_MAX", "geometry"]]
places = places.rename(columns={"NAME": "city", "ADM0NAME": "country", "POP_MAX": "population"})
eu_places = gpd.sjoin(places, europe[["geometry"]], predicate="within").drop(columns="index_right")
section("data")
print(len(europe), "European countries;", len(eu_places), "populated places")
print(europe.sort_values("density", ascending=False)[["name", "population", "area_km2", "density"]]
      .head(5).round(0).to_string(index=False))

# Simplify geometries (~1 km) to keep the HTML small - the browser draws every vertex.
# Clip to a European window first: Russia crosses the antimeridian, which would otherwise
# draw a band across the whole map.
europe = europe.clip((-25, 34, 50, 72))
europe["geometry"] = europe.to_crs("EPSG:3035").simplify(1000).to_crs("EPSG:4326")

# --- Step 2: the one-liner -------------------------------------------------------
m1 = europe.explore(column="density", scheme="Quantiles", k=5, cmap="YlGnBu",
                    tooltip=["name", "population", "density"],
                    legend_kwds={"caption": "People per km²", "colorbar": False, "fmt": "{:.0f}"},
                    tiles=ESRI_GRAY, attr=ESRI_GRAY_ATTR, location=[54, 15], zoom_start=4)
m1.save(OUT / "explore_density.html")

# --- Step 3: a hand-crafted Folium map --------------------------------------------
m = folium.Map(location=[54, 12], zoom_start=4, tiles=None, control_scale=True)
folium.TileLayer("OpenStreetMap", name="OpenStreetMap").add_to(m)
folium.TileLayer(ESRI_GRAY, attr=ESRI_GRAY_ATTR, name="Light grey (Esri)").add_to(m)
folium.TileLayer(
    tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
    attr="Tiles © Esri — Source: Esri, Maxar, Earthstar Geographics, and the GIS User Community",
    name="Esri World Imagery").add_to(m)

colormap = cm.LinearColormap(["#f7fcf0", "#7bccc4", "#0868ac"], vmin=float(europe.gdp_per_capita.quantile(0.05)),
                             vmax=float(europe.gdp_per_capita.quantile(0.95)),
                             caption="GDP per capita (USD, Natural Earth estimates)")

folium.GeoJson(
    europe[["name", "population", "density", "gdp_per_capita", "geometry"]],
    name="GDP per capita",
    style_function=lambda f: {"fillColor": colormap(f["properties"]["gdp_per_capita"]),
                              "color": "#475569", "weight": 0.6, "fillOpacity": 0.7},
    highlight_function=lambda f: {"weight": 2.5, "color": "#e8590c", "fillOpacity": 0.85},
    tooltip=folium.GeoJsonTooltip(fields=["name", "gdp_per_capita", "density"],
                                  aliases=["Country", "GDP per capita (USD)", "People per km²"],
                                  localize=True, sticky=True),
).add_to(m)
colormap.add_to(m)

cluster = plugins.MarkerCluster(name="Cities (clustered)").add_to(m)
for _, r in eu_places.iterrows():
    folium.Marker(
        location=[r.geometry.y, r.geometry.x],
        tooltip=r["city"],
        popup=folium.Popup(f"<b>{r['city']}</b><br>{r['country']}<br>Population: {r['population']:,}", max_width=220),
        icon=folium.Icon(color="darkblue", icon="info-sign"),
    ).add_to(cluster)

big = eu_places[eu_places["population"] >= 2_000_000]
fg = folium.FeatureGroup(name="Cities ≥ 2 million", show=False).add_to(m)
for _, r in big.iterrows():
    folium.CircleMarker([r.geometry.y, r.geometry.x], radius=4 + r["population"] / 1.5e6, color="#e8590c",
                        fill=True, fill_opacity=0.6, weight=1,
                        tooltip=f"{r['city']}: {r['population']:,}").add_to(fg)

plugins.Fullscreen().add_to(m)
plugins.MiniMap(toggle_display=True).add_to(m)
plugins.MeasureControl(primary_length_unit="kilometers").add_to(m)
folium.LayerControl(collapsed=False).add_to(m)
m.save(OUT / "europe_folium.html")

section("saved")
for f in sorted(OUT.glob("*.html")):
    print(f"{f.name:24s} {f.stat().st_size / 1024:,.0f} KB")