Making a basic map with Mapnik - Part 2: Projections

Introduction

My previous blog article described how to create a basic map showing a flight plan route with mapnik. In this part, I describe how to set the projection of the map.

Map Projections

The Earth is approximately spherical (actually an oblate spheroid, a sphere which is slightly flattened at the poles and bulging at the equator) but maps are two-dimensional and planar, so creating a two-dimensional map requires a mathematical transformation from geographical coordinates (latitude/longitude) to projected coordinates, typically Cartesian (x,y) coordinates.

The maps we have been generating so far have been using a projection whereby latitude and longitude coordinate values are effectively treated as if they were on a flat plane with Cartesian axes with a constant scaling between degrees longitude/latitude and x, y value (a cylindrical projection). With this projection, areas appear to be laterally stretched at latitudes far from the Equator since lines of longitude (meridians), converge toward the poles whereas lines of latitude remain parallel (which is why they're called parallels).

There are many other ways to project coordinates from an ellipsoid onto a plane, some more useful than others. The Mercator projection (a cylindrical projection, a slight variant of which is used by most web maps) was useful to navigators in the past because a straight line drawn on a Mercator projection map approximates a line which one would travel if one followed a constant compass bearing relative to true north (a rhumb line). Aeronautical charts use a Lambert conformal conic projection (LCC), on which a straight line approximates a Great Circle (the shortest path between two points on the surface of an ellipsoid), but along which a true north compass bearing varies. In the contiguous USA (CONUS), VFR sectional charts and terminal area charts use the LCC projection with standard parallels at 33°N and 45°N.

Expressing Coordinate Reference Systems

The "mapping" between the geographic coordinates on the earth (or other celestial body) and Cartesian coordinates on a map is defined by a Coordinate Reference System (CRS) or Spatial Reference System (SRS). There are a few standards for expressing CRS/SRS to geographic information systems; for example in the Well-Known Text (WKT) format WKT2:2019/ISO-19162:2019 [1], by referencing a projection logged in the EPSG (European Petroleum Survey Group) registry of geometric parameters [2], or as a "proj string" used by coordinate transformation software proj [3]. See footnote [A].

With mapnik, one specifies the map projection as the "srs" attribute of the map, either in the XML style file as an attribute of the <map> tag, or programatically as an attribute of a map object, using a proj string (mapnik itself uses the proj library). For example, in XML:

<Map srs="+proj=latlong +datum=WGS84">

or in a python program using the mapnik API:

map.srs = 'init=EPSG:4326'

Note from the example above that EPSG registry entries can be referenced directly but see footnote [B].

An Example

In my sample application I am plotting flight plans. Given the a standard Mercator projection, the flight plan for a flight from Washington National airport (KDCA) to Chicago O'Hare (KORD) appears as shown below:

Flight plan route KDCA to KORD: Mercator projection

Note that the Mason-Dixon line along Pennsylvania's southern border with Maryland and West Virginia (roughly along a parallel of 39°43′20″N but actually composed of multiple segments) is a straight line and the Great Lakes appear stretched laterally.

Since this is an aeronautical-type application, the LCC projection used in CONUS VFR sectional charts would seem an appropriate candidate, with standard parallels of 33°N and 45°N. We further specify the central meridian of the projection as in the middle of the route (i.e. in the middle of the meridians of KCDA and KORD, which I reckon to be about 82°28′23″W). This can be specified to mapnik as the following proj string

+proj=lcc +lon_0=-82.473 +lat_1=33 +lat_2=45

where lcc specifies Lambert Conformal Conical projection, +lon_0 specifies the central median, and +lat_1 and +lat_2 give the standard parallels. The result is as follows:

Flight plan route KDCA to KORD: Lambert Conformal Conic projection

The Mason-Dixon line is now curved and the Great Lakes better approximate their shape as seen on a globe.

A gotcha - bounding boxes

When creating a map, we need to specify the area that will be shown. This can be done by specifying the centre point of the map and a scale factor, or by specifying a bounding box containing the area we want to show. We usually want to specify the map centre or bounding box in geographic coordinates, but with mapnik we must specify these in projected coordinates.

The proj utility provided with the proj library can be used to convert between geographic and projected coordinates on the command line or in a script. See the man page. In a python program, the pyproj module can be used (of course, the proj library must be installed as a prerequisite). The following python program demonstrates the use of the pyproj module to transform a bounding box in geographical coordinates to projected coordinates and set the map bounds.

import mapnik
import pyproj

# Target SRS as a proj string
t_srs = '+proj=lcc +lon_0=-82.473 +lat_1=33 +lat_2=45'

m = mapnik.Map(500,500)                 # Create a 500x500 pixel map
mapnik.load_map(m, 'stylesheet.xml')    # Load map stylesheet
m.srs = t_srs                           # Set the SRS

# Compute and set the bounding box in
# projected coordinates. In geographic
# coordinates, let the lower left corner be
# minlat, minlon and the upper right corner
# be at maxlat, maxlon

transformer = \
    pyproj.Transformer.from_crs(
        'EPSG:4326',    # From CRS: WGS-84
        t_srs,          # Target CRS
        always_xy=True) # geographic coordinates in (x,y) order

bbox_xmin, bbox_ymin = transformer.transform(minlon, minlat)
bbox_xmax, bbox_ymax = transformer.transform(maxlon, maxlat)
m.zoom_to_bbox(
    mapnik.Box2D(bbox_xmin, bbox_ymin, bbox_xmax, bbox_ymax))

mapnik.render_to_file(m, 'output.png)

Coordinate Transformations between Data Sets

We have just talked about projecting to LCC map coordinates assuming the data sources are latitude/longitude coordinates, but each mapnik data source for mapnik might its own CRS. Indeed, there are several latitude/longitude systems based on different datums (data?). Mapnik assumes WGS-84 by default, which has now become a global de facto standard since it is used by the GPS (Global Positioning System) satellite navigation system but before satellite navigation, different countries tended to use different datums. If data sources are not in WGS-84 coordinates, the CRS can be specified to mapnik which will perform the necessary transformations from the source's CRS to the target map CRS.

Footnotes

[A]If you know the EPSG number of a CRS, you can get details such as its scope, intended usage, bounds and definition in various formats via the web service epsg.io [4]; e.g. the EPSG:4326 reference (WGS-84) can be accessed from https://epsg.io/4326.
[B]Being able to reference the EPSG registry entry is very convenient. The proj library (used either directly or indirectly through an application such as mapnik) contains EPSG registry information but as the registry changes with time, it can become out of date and newer entries may be "missing". Further, it can take time for entries to be added to the registry; EPSG was slow to add the "web Mercator" projection, for example, despite its widespread use in web mapping. The ability to use a "proj" string is therefore still useful, as well as enabling you to use a custom projection.

References

[1]Open Geospatial Consortium, Geographic information --- Well-known text representation of coordinate reference systems. http://docs.opengeospatial.org/is/18-010r7/18-010r7.html
[2]EPSG Geodetic Parameter Dataset https://epsg.org/
[3]PROJ https://proj.org/
[4]EPSG.io https://epsg.io/

social