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 ,
by referencing a projection logged in the EPSG (European Petroleum Survey
Group) registry of geometric parameters , or as a "proj string"
used by coordinate transformation software proj .
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:
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:
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)
References