Rendering WMS GetMap Images with Python and Mapnik

TL;DR: Map the incoming GetMap parameters onto a mapnik.Map: set map.srs from the requested CRS/SRS, add a Layer with a datasource and a Style of symbolizers, convert the BBOX to longitude-first order, call map.zoom_to_box(mapnik.Box2d(*bbox)), render into a mapnik.Image sized by WIDTH/HEIGHT, and return image.tostring("png") with Content-Type: image/png. The two things that most often go wrong are BBOX axis order and transparency encoding.

GetMap is the operation that turns a WMS from a metadata service into a map server. Once you have parsed and validated the request parameters — the discovery half of that flow is covered in How to Parse OGC WMS GetCapabilities XML in Python — the remaining work is translating those parameters into a rendering-engine call and streaming the resulting raster back to the client. Mapnik, the C++ rendering library behind OpenStreetMap’s cartography and many production WMS stacks, exposes a Python binding that makes this translation direct.

The Core Challenge

A GetMap request is a flat bag of HTTP query parameters — LAYERS, STYLES, BBOX, WIDTH, HEIGHT, CRS (or SRS on WMS 1.1.1), FORMAT, and TRANSPARENT. Mapnik, by contrast, has a stateful object model: a Map that owns a projection, an ordered list of Layer objects each bound to a datasource and a set of named styles, and Style objects that hold Rules and Symbolizers. Rendering is the act of mapping one onto the other correctly.

Three of those parameters carry hidden traps. BBOX arrives in an axis order that depends on the WMS version and the CRS — the single most common cause of blank GetMap responses. CRS/SRS must become a projection string that Mapnik understands and that matches the axis assumptions of the BBOX. And FORMAT together with TRANSPARENT decides the pixel encoding, where a wrong choice quietly destroys the alpha channel that overlay clients depend on.

The pipeline is short and linear, but each stage reads from a specific parameter, and skipping a translation step — most often the axis swap feeding zoom_to_box — produces an HTTP 200 response with a wrong or empty image rather than an error.

Production-Ready Code

The function below accepts the already-parsed and validated parameters and returns raw image bytes. It targets the Mapnik 3.x Python bindings (python3-mapnik), which ship with system Mapnik rather than being pure-Python — install libmapnik-dev and python3-mapnik from your distribution, or a platform wheel where available.

"""
wms_getmap_mapnik.py — render a WMS GetMap response with the Mapnik bindings.

Install (Debian/Ubuntu):
    apt-get install libmapnik-dev mapnik-utils python3-mapnik
"""
from __future__ import annotations

from dataclasses import dataclass

import mapnik

# proj4 definitions keyed by the CRS/SRS code a WMS client may send.
PROJ4: dict[str, str] = {
    "EPSG:4326": "+proj=longlat +datum=WGS84 +no_defs",
    "EPSG:3857": (
        "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 "
        "+k=1 +units=m +nadgrids=@null +wktext +no_defs"
    ),
}

# CRS codes whose OGC axis order is latitude-first under WMS 1.3.0.
LAT_FIRST_130: set[str] = {"EPSG:4326"}


@dataclass
class GetMapParams:
    layers: list[str]
    bbox: tuple[float, float, float, float]  # in the request's own axis order
    width: int
    height: int
    crs: str                                 # from CRS (1.3.0) or SRS (1.1.1)
    fmt: str = "image/png"
    version: str = "1.3.0"
    transparent: bool = False


def _bbox_to_xy(p: GetMapParams) -> tuple[float, float, float, float]:
    """Return the BBOX as (minx, miny, maxx, maxy) in longitude/easting-first order.

    mapnik.Box2d is always x-first, but WMS 1.3.0 delivers EPSG:4326 as
    minLat,minLon,maxLat,maxLon. Swap the pairs so the map frames the
    intended region instead of an empty extent.
    """
    a, b, c, d = p.bbox
    if p.version.startswith("1.3") and p.crs in LAT_FIRST_130:
        return b, a, d, c  # lat,lon -> lon,lat
    return a, b, c, d


def build_style() -> mapnik.Style:
    """A minimal polygon-fill + outline style. Replace with SLD-derived rules."""
    style = mapnik.Style()
    rule = mapnik.Rule()

    fill = mapnik.PolygonSymbolizer()
    fill.fill = mapnik.Color("#8db2a3")
    fill.fill_opacity = 0.85
    rule.symbols.append(fill)

    outline = mapnik.LineSymbolizer()
    outline.stroke = mapnik.Color("#33454e")
    outline.stroke_width = 0.6
    rule.symbols.append(outline)

    style.rules.append(rule)
    return style


def render_getmap(p: GetMapParams, shapefile_path: str) -> tuple[str, bytes]:
    """Translate parsed GetMap parameters into a Mapnik render.

    Returns (content_type, image_bytes) ready to write to the HTTP response.
    """
    # 1. The Map owns the output size and the projection of the BBOX + image.
    m = mapnik.Map(p.width, p.height)
    m.srs = PROJ4[p.crs]

    # 2. Transparency: a zero-alpha RGBA background lets overlays composite.
    if p.transparent and "png" in p.fmt.lower():
        m.background = mapnik.Color("rgba(0,0,0,0)")
    else:
        m.background = mapnik.Color("white")

    # 3. Datasource + layer + style. The layer SRS may differ from the map
    #    SRS; Mapnik reprojects features on the fly when they disagree.
    m.append_style("default", build_style())

    layer = mapnik.Layer(",".join(p.layers))
    layer.datasource = mapnik.Shapefile(file=shapefile_path)
    layer.srs = PROJ4["EPSG:4326"]        # source data stored in lon/lat
    layer.styles.append("default")
    m.layers.append(layer)

    # 4. Frame the requested extent (map-SRS x,y order).
    minx, miny, maxx, maxy = _bbox_to_xy(p)
    m.zoom_to_box(mapnik.Box2d(minx, miny, maxx, maxy))

    # 5. Render. scale_factor scales line widths, markers and text for DPI.
    image = mapnik.Image(p.width, p.height)
    scale_factor = 1.0                    # 2.0 for retina / ~181 DPI tiles
    mapnik.render(m, image, scale_factor)

    # 6. Serialise to the encoding implied by FORMAT + TRANSPARENT.
    if "jpeg" in p.fmt.lower() or "jpg" in p.fmt.lower():
        return "image/jpeg", image.tostring("jpeg90")
    encoding = "png32" if p.transparent else "png"
    return "image/png", image.tostring(encoding)


if __name__ == "__main__":
    params = GetMapParams(
        layers=["countries"],
        bbox=(-90.0, -180.0, 90.0, 180.0),   # 1.3.0 EPSG:4326 = lat,lon
        width=800,
        height=400,
        crs="EPSG:4326",
        version="1.3.0",
        transparent=True,
    )
    content_type, png = render_getmap(params, "ne_110m_admin_0_countries.shp")
    with open("getmap.png", "wb") as fh:
        fh.write(png)
    print(f"{content_type}: {len(png)} bytes")

Step-by-Step Walkthrough

Map construction and map.srs. mapnik.Map(p.width, p.height) fixes the output raster dimensions directly from the WIDTH and HEIGHT parameters, so the rendered image already matches what the client asked for and needs no resampling. Setting m.srs to the proj4 string for the requested CRS tells Mapnik which projection the BBOX coordinates and the output pixels live in. The PROJ4 lookup keeps this to a whitelist of supported codes; an unrecognised CRS should have been rejected during parameter validation with a ServiceException, not reached the renderer.

The axis swap in _bbox_to_xy. mapnik.Box2d takes (minx, miny, maxx, maxy) and always treats the first ordinate as x — longitude or easting. WMS 1.3.0, however, honours the OGC axis definition of EPSG:4326, which is latitude-first, so the client sends minLat,minLon,maxLat,maxLon. Feeding that straight into Box2d frames a degenerate or wrong extent and the render comes back empty. The helper swaps the pairs only for the version/CRS combinations that need it. The full set of codes that require swapping, and how to detect them programmatically, is documented in the SRS and Coordinate Reference System Handling guide.

Transparency and the background colour. m.background = mapnik.Color("rgba(0,0,0,0)") makes the canvas fully transparent so that when the client stacks this layer over a base map the areas outside your features show through. This only matters for formats with an alpha channel, so the code guards it on "png" being in FORMAT; a JPEG response is always opaque and falls back to a white fill.

Style, layer and datasource. build_style() assembles a Style containing one Rule with a PolygonSymbolizer (a translucent fill) and a LineSymbolizer (the outline). The style is registered on the map under the name "default" with m.append_style, and the Layer references it by that name through layer.styles.append("default"). The layer is bound to a mapnik.Shapefile datasource; in production this is where you would swap in a mapnik.PostGIS datasource pointed at the tables named by LAYERS. Because layer.srs is declared independently of m.srs, Mapnik reprojects the features on the fly when the data projection differs from the requested output projection.

Rendering and scale_factor. mapnik.render(m, image, scale_factor) draws the framed map into the mapnik.Image. The scale_factor argument multiplies every symbolizer dimension — line widths, marker sizes, text — which is how you produce high-DPI tiles without editing the style; a factor of 2.0 doubles pixel density. If you prefer to write straight to disk instead of holding bytes in memory, mapnik.render_to_file(m, "out.png", "png") is the equivalent one-shot call.

Serialisation to bytes. image.tostring(encoding) returns the encoded raster as a bytes object. The encoding is chosen from FORMAT: "jpeg90" for JPEG, "png32" for full 32-bit RGBA when transparency is requested, and plain "png" otherwise. Returning the matching Content-Type alongside the bytes lets the WSGI or ASGI layer set the header without re-inspecting the payload.

Verification

Render a single global frame and confirm the byte count and PNG signature from the shell:

python wms_getmap_mapnik.py
python - <<'PY'
data = open("getmap.png", "rb").read()
print("signature ok:", data[:8] == b"\x89PNG\r\n\x1a\n")
print("bytes:", len(data))
PY

Expected output:

image/png: 21874 bytes
signature ok: True
bytes: 21874

The exact byte count varies with the datasource and Mapnik build, but a valid PNG signature and a non-trivial size confirm the pipeline produced a real raster rather than an empty or corrupt buffer. To confirm the extent is correct rather than merely non-empty, request a small BBOX over a known feature and check that the visible geometry lands where you expect.

Gotchas & Edge Cases

Back to Understanding OGC Web Map Service Specifications

Related

• Understanding OGC Web Map Service Specifications — the full GetMap parameter contract and where rendering sits in the request lifecycle
• How to Parse OGC WMS GetCapabilities XML in Python — extracting the layer names and CRS lists that feed a GetMap request
• SRS and Coordinate Reference System Handling — axis-order rules and reprojection logic behind the BBOX swap