Web Mapping & Interactive Visualization with Python

Processing and analysis are only half of a geospatial project — at some point the results have to reach a browser. This domain, part of the wider python-geospatial.com reference, covers the visualization layer of the Python stack: turning a validated GeoDataFrame into an interactive map a stakeholder can pan, zoom, and click. It sits downstream of the analytical work — the libraries covered in Mastering Core Geospatial Python Libraries and the pipelines in Geospatial Data Ingestion & Processing Workflows — and connects them to the rendering technologies that run client-side. Four topics structure the work: rapid prototyping with Interactive Maps with Folium, production vector rendering with MapLibre GL Vector Web Maps, scalable delivery through Vector Tile Pipelines with PMTiles, and getting the result in front of users via Geospatial Dashboards & App Deployment. The recurring discipline across all three is knowing exactly where a coordinate transformation happens: analysis stays in a metric projected CRS, and the hand-off to the browser is the single, deliberate moment you drop to EPSG:4326.

Web mapping delivery pipeline A four-stage flow from a projected GeoDataFrame, through Python encoding, to tiles or GeoJSON, ending in a browser map, with a note to reproject to EPSG 4326 before delivery. From GeoDataFrame to the Browser Spatial data GeoDataFrame / COG Python encode Folium / tiling Tiles or GeoJSON PMTiles / MBTiles Browser map MapLibre / Leaflet Analyse in a projected CRS — reproject to EPSG:4326 only at the delivery boundary Web clients assume lon/lat tiles; metric work stays in UTM upstream
The visualization layer sits downstream of analysis: data is encoded in Python, materialized as tiles or GeoJSON, and rendered client-side.

Ecosystem Architecture & Dependency Management

The Python web-mapping ecosystem splits cleanly into two responsibilities: generating map artifacts in Python and rendering them in JavaScript. Unlike the analytical stack in Mastering Core Geospatial Python Libraries, where GeoPandas, Shapely, and Rasterio all bottom out in the same GDAL/PROJ/GEOS C libraries, the visualization layer is a hand-off across a language boundary. Python's job ends when it has emitted a portable artifact — HTML, GeoJSON, or a tile archive — and a JavaScript renderer takes over in the browser. Getting the split right is what lets a 50-row notebook demo and a 5-million-feature production map share the same upstream analysis code.

Three tools anchor the Python side, each targeting a different scale:

Adjacent libraries fill specialist niches: ipyleaflet for bidirectional widgets inside Jupyter, pydeck/lonboard for GPU-accelerated point clouds, and kepler.gl exports for exploratory large-data views. All of them follow the same rule — Python encodes, JavaScript renders — so the CRS and payload discipline below applies regardless of which renderer you pick.

A more useful way to compare them than "which library" is what artifact each one leaves behind, because the artifact determines how the thing is hosted, cached, versioned and debugged long after the code is written:

Python tool Artifact it emits Hosting it needs Where it stops
Folium one self-contained .html with data inline none — a file payload past a few MB
GeoPandas .to_file(..., "GeoJSON") a .geojson asset static host + CORS one download, no zoom filtering
tippecanoe + pmtiles a single .pmtiles archive static host with range requests rebuild required to change data
rio cogeo / GDAL a Cloud Optimized GeoTIFF static host with range requests client must speak COG or use a tile service
FastAPI / TiTiler a live {z}/{x}/{y} endpoint a running process scales with request volume, not data size

ipyleaflet deserves a note because it is the one genuine exception to the one-way hand-off. It is a Jupyter widget, so the browser map and the Python kernel stay connected over the notebook's comm channel: a user drawing a polygon on the map fires a Python callback, and the result can flow straight into a GeoDataFrame. That is impossible in Folium, whose output has no live kernel behind it. The cost is that an ipyleaflet map only exists while the notebook is running — there is no file to send anyone — which makes the two libraries complements rather than competitors.

The browser-side dependency question is easy to defer and expensive to get wrong. Folium vendors its Leaflet references as CDN <script> tags inside the generated HTML, so a saved map silently depends on those CDNs being reachable at view time, on whatever network the reader is on. For anything that must work offline, behind a restrictive corporate proxy, or for longer than a CDN's URL scheme survives, vendor the JavaScript alongside the artifact and rewrite the references — a build step of a few lines that converts a fragile document into a durable one.

Install the rendering-adjacent toolchain in an isolated environment. Pin versions, because Folium's Leaflet bundle and the tile-generation libraries evolve independently:

# environment.yml (conda-forge channel)
# name: webmap
# dependencies:
#   - python=3.11
#   - geopandas=0.14.*
#   - folium=0.16.*
#   - pmtiles=3.2.*
#   - gdal=3.8.*          # provides ogr2ogr for tile generation
#   - pip
#   - pip: ["pyogrio>=0.7"]

import importlib

WEBMAP_PACKAGES = ["geopandas", "folium", "pmtiles", "shapely", "pyproj"]


def verify_webmap_stack() -> None:
    missing = [pkg for pkg in WEBMAP_PACKAGES if not importlib.util.find_spec(pkg)]
    if missing:
        raise RuntimeError(f"Missing web-mapping dependencies: {', '.join(missing)}")
    print("Web mapping stack ready.")


verify_webmap_stack()

The tippecanoe tile generator (used for large vector tile builds) is a C++ binary, not a Python package — install it via your system package manager or a container layer and call it as a subprocess. Treat it like GDAL: a native dependency the Python code orchestrates rather than imports.

Core Concepts & Data Model

Every web map is a stack of layers drawn over a base map, positioned by a viewport (center plus zoom). The data you supply is either raster (pre-rendered image tiles) or vector (GeoJSON or vector tiles the client styles on the fly). The single most important rule: web clients expect geographic coordinates. Leaflet and MapLibre consume EPSG:4326 longitude/latitude and internally project to EPSG:3857 (Web Mercator) for tile placement. Your analysis should happen in a projected CRS — see Coordinate Systems with PyProj — and only convert to 4326 at the moment of export.

import folium
import geopandas as gpd

# floodplain_boundary was analysed in a metric UTM CRS upstream
floodplain_boundary = gpd.read_file("floodplain_utm.gpkg")

# Reproject to WGS84 ONLY for web delivery
floodplain_wgs84 = floodplain_boundary.to_crs(epsg=4326)

centroid = floodplain_wgs84.geometry.union_all().centroid
fmap = folium.Map(location=[centroid.y, centroid.x], zoom_start=11, tiles="CartoDB positron")
folium.GeoJson(floodplain_wgs84, name="Floodplain").add_to(fmap)
fmap.save("floodplain_map.html")

Note the [centroid.y, centroid.x] ordering — Leaflet and Folium take [lat, lon], the opposite of the (x, y) ordering used throughout the analytical libraries. This axis flip is the most common first bug in web mapping code, and it mirrors the always-xy pitfalls documented in Coordinate Systems with PyProj: the same longitude/latitude pair is written (x, y) in the geometry engine and [y, x] in the map widget.

The second conceptual fork is how the data reaches the renderer. There are two data models, and choosing between them is the central architectural decision of any web map:

Inline GeoJSON versus vector tiles: browser cost against dataset size A line chart of client cost against dataset size. Inline GeoJSON cost rises linearly with feature count, while vector-tile cost stays flat and low. The two curves cross at roughly five megabytes, or tens of thousands of features: below the crossover inline GeoJSON is simplest, above it vector tiles win. Two data models, two cost curves inline-GeoJSON zone vector-tile zone client cost per frame dataset size / feature count crossover ≈ 5 MB / tens of thousands of features Inline GeoJSON — cost grows with features Vector tiles — constant per-frame cost
Inline GeoJSON pays a per-frame cost that scales with feature count; vector tiles fetch only the current viewport, holding that cost flat. The curves cross near 5 MB — the practical threshold for switching to tiles.

Underneath both models sits the tile pyramid, and a little of its arithmetic explains most of the surprises. The Web Mercator grid is a square: at zoom 0 the whole world is one 256-pixel tile, and each zoom level splits every tile into four, so zoom z holds 4^z tiles. Zoom 14 — a typical maximum for a city-scale vector layer — is 268 million tiles worldwide, which is why nobody pre-renders the whole world and why tippecanoe only cuts tiles where features actually exist. Ground resolution follows directly: about 156,543 metres per pixel at zoom 0, halving each level, and multiplied by the cosine of the latitude. That last factor is the reason a 60-metre buffer looks generous in Oslo and invisible in Nairobi at the same zoom.

import math

EQUATOR_RES = 156543.03392   # metres per pixel at zoom 0, 256 px tiles


def ground_resolution(zoom: int, latitude: float) -> float:
    """Metres per screen pixel at a given zoom and latitude."""
    return EQUATOR_RES * math.cos(math.radians(latitude)) / (2 ** zoom)


def tiles_covering(bounds_wgs84, zoom: int) -> int:
    """How many tiles a WGS84 bounding box needs at one zoom level."""
    minx, miny, maxx, maxy = bounds_wgs84
    n = 2 ** zoom

    def to_tile(lon, lat):
        x = int((lon + 180.0) / 360.0 * n)
        lat_rad = math.radians(max(min(lat, 85.0511), -85.0511))
        y = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
        return x, y

    x0, y1 = to_tile(minx, miny)      # note: tile y increases southward
    x1, y0 = to_tile(maxx, maxy)
    return (abs(x1 - x0) + 1) * (abs(y1 - y0) + 1)


berlin_bounds = (13.09, 52.34, 13.76, 52.68)
print(f"{ground_resolution(14, 52.5):.2f} m/px")          # 5.81 m/px
print(tiles_covering(berlin_bounds, 14), "tiles at z14")  # 336 tiles at z14

Two conventions on top of that grid cause real bugs. The first is tile addressing: the XYZ scheme used by OpenStreetMap, Leaflet and MapLibre counts rows from the north, while the older TMS scheme counts from the south, so y_tms = 2^z − 1 − y_xyz. Mixing them produces a map that is correct in longitude and mirrored in latitude — features in the right column, the wrong hemisphere. MBTiles archives store TMS internally while serving XYZ, which is exactly where the confusion usually enters. The second is the latitude clamp: Web Mercator cannot represent the poles, so the grid is truncated at ±85.0511° to keep the world square. Data above that line — Arctic sea ice, northern shipping routes, Antarctic research stations — is not merely distorted, it is outside the tile grid, and a polar dataset needs a different projection and a renderer that supports one.

Everything else — base map choice, tooltips, layer toggles, legends — layers on top of this single decision. Get the data model right for your scale first; styling is cheap to change afterwards.

Key Operations & Vectorized Workflows

Four operations cover the majority of practical web-mapping work:

1. Render a vector overlay. Convert a GeoDataFrame to GeoJSON and add it as a styled layer. For interactive choropleths, see Folium Choropleth from a GeoDataFrame.

2. Cluster dense point sets. Thousands of raw markers freeze the browser; cluster them server-side or with Leaflet's clustering plugin, covered in Clustering Map Markers with Folium.

3. Stream GeoJSON to a vector renderer. MapLibre GL JS styles raw GeoJSON sources — see Serving GeoJSON to MapLibre GL JS.

4. Build tiles for scale. Past ~5 MB of GeoJSON, switch to vector tiles. Generating PMTiles from GeoParquet shows the cloud-native path.

5. Attach interactivity. Bind tooltips and popups so features are self-describing on hover or click — this is what makes a map a communication tool rather than a picture.

import geopandas as gpd

sensors = gpd.read_file("air_quality_sensors.gpkg").to_crs(epsg=4326)

# Trim precision before serialization: 6 decimals ≈ 0.11 m, plenty for web
sensors["geometry"] = sensors.geometry.set_precision(1e-6)

# Drop heavy attribute columns the client never displays
web_cols = ["sensor_id", "pm25", "geometry"]
sensors[web_cols].to_file("sensors_web.geojson", driver="GeoJSON")

Precision trimming and column pruning routinely cut payload size by half — the cheapest performance win in web mapping. The two levers are independent: set_precision quantizes coordinates to a grid (6 decimal degrees is ~0.11 m, finer than any web screen resolves), and dropping unused columns removes attribute bytes the client never reads.

For the interactivity step, keep field selection explicit so you never ship raw internal columns to the browser:

import folium

fmap = folium.Map(location=[52.52, 13.40], zoom_start=11, tiles="CartoDB positron")

folium.GeoJson(
    "sensors_web.geojson",
    name="Air quality",
    tooltip=folium.GeoJsonTooltip(fields=["sensor_id", "pm25"], aliases=["Sensor", "PM2.5"]),
    marker=folium.CircleMarker(radius=5, fill=True),
).add_to(fmap)

folium.LayerControl().add_to(fmap)  # requires name= on every layer above
fmap.save("air_quality.html")

The LayerControl toggle only lists layers that were given a name=; unnamed overlays are invisible to it, a gotcha covered again below.

6. Measure the payload before you ship it. Every architectural decision on this page reduces to one number — the bytes a browser must download to draw the first viewport — and it takes four lines to measure rather than guess. Put the measurement in the pipeline, not in a notebook you ran once, so that a data refresh which quietly triples the feature count fails the build instead of the reader's browser.

import geopandas as gpd

def payload_report(gdf: gpd.GeoDataFrame, label: str) -> int:
    """Bytes of GeoJSON, and where that lands against the delivery thresholds."""
    payload = gdf.to_json().encode("utf-8")
    mb = len(payload) / 1e6
    verdict = (
        "inline GeoJSON is fine" if mb < 5
        else "borderline — trim, or move to tiles" if mb < 15
        else "tiles required"
    )
    print(f"{label:<20} {len(gdf):>8,} features  {mb:>6.2f} MB  → {verdict}")
    return len(payload)

payload_report(sensors, "air quality sensors")
# air quality sensors     4,182 features    1.37 MB  → inline GeoJSON is fine

Two habits make that number honest. Measure after the trimming steps, not before, or you will move to a tile pipeline you did not need. And measure the serialised form rather than the in-memory GeoDataFrame, because the relationship between the two is not fixed: a layer of long, dense linestrings serialises to many times its pandas footprint, while a layer of points with wide attribute tables barely grows.

There is a seventh operation that is not technical and gets skipped anyway: satisfy the basemap's licence. OpenStreetMap-derived tiles require visible credit to OpenStreetMap contributors under the ODbL, commercial providers require an API key and forbid tile scraping, and most free public endpoints are explicitly not for production traffic. Folium enforces part of this by refusing to build a custom TileLayer without an attr string; nothing enforces the rest. Decide the tile source at design time — self-hosted, commercial, or none at all — because retrofitting a basemap change late means re-checking every colour choice against a different background.

CRS / Projection Considerations

The web-mapping CRS contract is narrow but unforgiving. Tiles are addressed in the Web Mercator grid (EPSG:3857), but the data you hand to a client should be tagged EPSG:4326; the renderer performs the Mercator projection. Three failure modes recur:

import geopandas as gpd

parcels = gpd.read_file("parcels.gpkg")  # EPSG:32633 (UTM 33N)

# Compute density in the projected CRS where area is in metres²
parcels["density"] = parcels["population"] / (parcels.geometry.area / 1e6)

# THEN reproject the finished layer for the browser
parcels_web = parcels.to_crs(epsg=4326)
assert parcels_web.crs.to_epsg() == 4326

Beyond those three, two rules from the GeoJSON specification itself (RFC 7946) catch pipelines that were otherwise correct. The first is the antimeridian: a geometry crossing 180° longitude, written as a single ring whose coordinates jump from +179 to −179, is drawn by every web renderer as a band stretching the wrong way around the entire world. The specification's instruction is to split such geometries at the antimeridian into a MultiPolygon or MultiLineString, and the practical advice is to do it deliberately in Python rather than hope — anything Pacific-facing, Russian, Fijian or New Zealand-adjacent will eventually hit this. The second is winding order: RFC 7946 requires exterior rings counter-clockwise and holes clockwise. Leaflet ignores winding, MapLibre mostly does, but tile cutters and some server-side renderers do not, and a wrongly wound ring turns into an inverted fill that shades everything except the polygon.

import geopandas as gpd
import shapely

territories = gpd.read_file("territorial_waters.gpkg").to_crs(epsg=4326)

# Flag rings that span more than half the globe — the antimeridian signature
minx, maxx = territories.bounds["minx"], territories.bounds["maxx"]
suspect = territories[(maxx - minx) > 180]
print(f"{len(suspect)} geometries appear to cross the antimeridian")   # 3 geometries

# Enforce RFC 7946 winding before handing anything to a tile cutter
territories["geometry"] = shapely.force_ccw(territories.geometry.to_numpy())

One question this section keeps raising is which projected CRS to analyse in, and the answer depends on the measurement. Areas and densities want an equal-area projection — a national equal-area grid, or an Albers/Lambert Azimuthal projection centred on the study area. Distances and buffers want a local UTM zone, chosen automatically where the data spans one zone (see Choosing a UTM Zone Automatically in Python) and replaced by a national grid where it does not. Web Mercator is never the answer for either, and its use as an analysis CRS is common enough that it deserves restating: at 60° latitude it inflates areas by a factor of four.

A final note for anyone integrating an existing service: not every tile source is Web Mercator. WMS endpoints are frequently published in a national grid or in EPSG:4326, and there is a separate 4326-based tile scheme (two tiles at zoom 0 rather than one) used by some agencies. Leaflet can consume a WMS layer directly, reprojecting nothing and letting the server do the work; MapLibre cannot mix projections in one map. When a legacy WMS is a hard requirement, that constraint alone may decide the renderer.

Production Patterns & Performance

The dividing line in production web mapping is payload size. Below a few megabytes, inline GeoJSON is simplest. Above it, you need tiles so the client only downloads the viewport at the current zoom. The cloud-native answer is PMTiles — a single-file tile archive served over HTTP range requests with no tile server, conceptually identical to the Cloud-Native Geospatial Formats used upstream.

The build side of a PMTiles pipeline is a short subprocess orchestration — GeoParquet in, an HTTP-servable archive out:

import subprocess

# tippecanoe reads GeoJSON/FlatGeobuf; export the web-ready layer first
# (parcels already reprojected to EPSG:4326 and precision-trimmed)
subprocess.run(
    [
        "tippecanoe",
        "-o", "parcels.pmtiles",
        "--layer=parcels",
        "--maximum-zoom=14",
        "--minimum-zoom=6",
        "--drop-densest-as-needed",
        "--force",
        "parcels_web.geojson",
    ],
    check=True,
)
# Host parcels.pmtiles on any static bucket; MapLibre reads it via the pmtiles:// protocol

check=True turns a non-zero exit code into a raised CalledProcessError, so a broken tile build fails the pipeline loudly instead of silently shipping an empty map.

Performance on a web map is three separate numbers that people collapse into one. Transfer is bytes over the wire, and compression fixes most of it. Parse and build is the main-thread time spent turning those bytes into geometry and GPU buffers or DOM nodes — this is the one that produces a frozen tab, it is not helped by a faster connection, and it scales with feature and vertex count rather than file size. Frame cost is what happens on every pan and zoom, and it is driven by how many features are visible at once, which is why a dataset that is fine at zoom 14 stutters at zoom 8 where everything is in view. Vector tiles attack all three; simplification attacks the second and third; compression attacks only the first. Diagnosing which number is actually your problem before choosing a remedy saves a lot of wasted pipeline work.

Serving the assets correctly is mostly a matter of getting three headers right, and each has a specific failure signature:

The last production concern is the one that never appears in a benchmark: who can read the map. An interactive map is an image with no alternative text unless you provide one, so pair every published map with the numbers behind it — a table, a download link, or a summary sentence — so the information survives for a reader using a screen reader or a keyboard. Choose a colour ramp that stays ordered in greyscale and for the common forms of colour-vision deficiency, and never encode a category by colour alone where the distinction matters; add a shape, a pattern, or a label. Check contrast against the actual basemap you shipped, not against white, because a palette tuned on a light background disappears over satellite imagery.

Common Mistakes

Frequently Asked Questions

Should I use Folium or MapLibre GL JS? Use Folium for notebooks, internal dashboards, and anything you want as a one-file HTML artifact. Use MapLibre GL JS when you need smooth vector rendering, data-driven styling, or datasets large enough to require vector tiles. Folium renders raster-style Leaflet maps; MapLibre is a WebGL vector renderer. The pragmatic rule is dataset size: inline GeoJSON in Folium below a few megabytes, GeoJSON in MapLibre when you want WebGL styling but the data still fits in memory, and PMTiles vector tiles once the data outgrows a single payload.

Choosing a web-mapping stack by dataset size Three tiers, left to right by increasing dataset size. Small datasets under about five megabytes use Folium with inline GeoJSON in a single HTML file. Medium datasets that fit in browser memory use MapLibre GL JS with a GeoJSON source. Large datasets of millions of features use tippecanoe to build PMTiles vector tiles served from a static bucket. Pick the stack by dataset size SMALL < ~5 MB · tens of thousands of features Folium inline GeoJSON in one self-contained HTML file MEDIUM still fits in browser memory MapLibre GL JS GeoJSON source, WebGL data-driven styling LARGE millions of features PMTiles vector tiles tippecanoe build, served from a static bucket increasing dataset size / feature count
The same upstream analysis feeds three delivery stacks; dataset size alone usually decides which one. Styling and interactivity layer on afterwards.

At what size should I stop using GeoJSON and switch to tiles? A practical threshold is ~5 MB of uncompressed GeoJSON or a few tens of thousands of features. Beyond that, the browser stalls parsing and rendering, and vector tiles (PMTiles/MBTiles) become worth the build step.

Do I need a tile server? Not anymore. PMTiles serves an entire tileset from a single static file using HTTP range requests, so a plain object store or CDN replaces a running tile server for most read-only workloads.

Why does my map show the data in the wrong place? Almost always a CRS or axis-order issue: either the data wasn't reprojected to EPSG:4326, or coordinates were passed as [lon, lat] where Leaflet expects [lat, lon].

Can I put a raster analysis result on a web map? Yes — export it as a Cloud Optimized GeoTIFF and either pre-render image tiles or let a client like MapLibre/Leaflet fetch it. Reproject to EPSG:3857 or supply it as a 4326 source depending on the renderer.

How much of this work belongs in Python and how much in JavaScript? Draw the line at decisions that depend on the data versus decisions that depend on the viewport. Anything derived from the data — classification breaks, aggregation, simplification tolerance, which attributes travel — belongs in Python, where it is testable, reproducible and version-controlled. Anything that depends on what the user is looking at — which features are on screen, how thick a line reads at this zoom, whether a label fits — belongs in the renderer, because Python cannot know it. Pipelines get painful when that line moves: baking zoom-dependent styling into Python means re-rendering for every visual tweak, while pushing classification into JavaScript means the map and the report disagree about what "top quintile" means.

Do I have to rebuild tiles every time the data changes? For PMTiles and any pre-cut archive, yes — the archive is an immutable build artifact, and the workflow is the same as a compiled asset: rebuild, publish under a new versioned name, swap the reference. That is usually fine, because the build is minutes and the data changes daily at most. When the data changes by the second, or when each user sees a different filtered subset, pre-cutting is the wrong model entirely and you want a live tile endpoint generating tiles per request — the path covered in Geospatial Dashboards & App Deployment.

Why does my map look fine locally and break when deployed? Four causes account for most of it, in rough order of frequency. Missing CORS headers on the data file, so the page loads and the layer silently does not. A wrong content type from the static host, which some browsers tolerate and others reject. A host or proxy that does not honour byte-range requests, which turns PMTiles and COG reads into full-file downloads. And relative paths that resolved against your working directory locally but not against the deployed URL structure. All four are invisible in the page itself and obvious in the browser's network panel, which is the first place to look rather than the last.

Can I put a web map in a printed report or a slide deck? Render it to an image rather than trying to print the live page — interactive maps print unpredictably, often with tiles half-loaded. A headless browser screenshot of the saved HTML is the reliable route, and it also gives you a repeatable figure in CI. Bear in mind that the design constraints differ: a static figure needs its own scale bar, north arrow and legend, all of which an interactive map gets away with omitting because the reader can zoom.

How do I keep a large point dataset responsive? Don't ship raw points. Either cluster them — server-side or with Leaflet's marker-cluster plugin, covered in Clustering Map Markers with Folium — or aggregate into bins (H3 hexes, grid cells) upstream so the browser renders a few thousand summary shapes instead of millions of markers. For genuinely large point sets, a WebGL renderer (MapLibre, deck.gl) beats the DOM-based Leaflet marker layer by orders of magnitude.

Related