If you are viewing this file in preview mode, some links won't work. Find the fully featured Jupyter Notebook file on the website of Prof. Jens Flemming at Zwickau University of Applied Sciences. This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Cartopy and Mplleaflet

Cartopy is a Python library for creating geographic maps. It tightly integrates into Matplotlib. Mplleaflet combines Matplotlib drawings with interactive maps.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import mplleaflet

The letters crs are the abbreviation of coordinate reference system. Cartopys crs module handles all coordinate transforms. In particular, it tells Matplotlib how to transform geographical coordinates to screen coordinates.

Cartopy basics

Creating a map requires only few steps:

• create a Matplotlib figure,
• create a GeoAxes object,
• add content to the map.

fig = plt.figure()
ax = fig.add_axes((0, 0, 1, 1), projection=ccrs.PlateCarree())

ax.coastlines()

plt.show()

ax is a Cartopy GeoAxes object, because the projection is a Cartopy projection. There are many different projection types available, see Cartopy projection list. Note that GeoAxes is derived from the usual Matplotlib Axes and thus provides very similar functionality.

Most projection types take arguments for specifing map center and some parameters. To get a map of a smaller geographic region use GeoAxes.set_extent. Detail level of coastlines can be controlled with resolution argument, which has to be 110m (default), 50m or 10m. Here is a map of Germany:

fig = plt.figure()
ax = fig.add_axes((0, 0, 1, 1), projection=ccrs.Orthographic(10.5, 51.25))
ax.set_extent([5.5, 15.5, 47, 55.5])

ax.coastlines(resolution='50m')

plt.show()

Adding more content

Cartopy gets its data from Natural Earth, which provides public domain map data. But other sources can be used, too. Map data is downloaded by Cartopy as needed and then reused if needed again. To add Natural Earth content to a map call GeoAxes.add_feature.

fig = plt.figure(figsize=(8,8))
ax = fig.add_axes((0, 0, 1, 1), projection=ccrs.Orthographic(10.5, 51.25))
ax.set_extent([5.5, 15.5, 47, 55.5])

ax.coastlines(resolution='50m')

ax.add_feature(
    cfeature.NaturalEarthFeature(
        name='land',
        scale='50m',
        category='physical'
    ),
    facecolor='#e0e0e0'
)
ax.add_feature(
    cfeature.NaturalEarthFeature(
        name='admin_0_boundary_lines_land',
        scale='50m',
        category='cultural'
    ),
    edgecolor='r',
    facecolor='#e0e0e0'
)

plt.show()

Available options for the name argument can be guessed from the download section of Natural Earth. Look at the filenames in the download links. Arguments other than name, resolution and category are passed on to Matplotlib.

Plotting onto the map

All Matplotlib plotting functions are available for extending the map. Coordinates can be provided with respect to arbitrary cartographic projections, if the correct transform is chosen. If latitude, longitude values are used, PlateCarree is the right choice. For all other projection types coordinates have to be provided in meters.

fig = plt.figure()
ax = fig.add_axes((0, 0, 1, 1), projection=ccrs.Orthographic(10.5, 51.25))
ax.set_extent([5.5, 15.5, 47, 55.5])

ax.coastlines(resolution='50m')

ax.plot(13.404, 52.518, 'ob', transform=ccrs.PlateCarree())
ax.text(13.404, 52.9, 'Berlin', ha='center', va='center', transform=ccrs.PlateCarree())

ax.plot(-200000, 0, 'or', transform=ccrs.Orthographic(10.5, 51.25))
ax.text(-200000, 80000, '200km west\nof map center', ha='center', va='center',
        transform=ccrs.Orthographic(10.5, 51.25))

plt.show()

Take care: On Cartopy maps plot does not connect two points by a straight line, but by the shortest path with respect to the chosen coordinate system.

fig = plt.figure()
ax = fig.add_axes((0, 0, 1, 1), projection=ccrs.PlateCarree())

ax.coastlines()

lat1 = 40
lon1 = -90
lat2 = 30
lon2 = 70

ax.plot([lon1, lon2], [lat1, lat2], '-ob', transform=ccrs.PlateCarree())
ax.plot([lon1, lon2], [lat1, lat2], '-r', transform=ccrs.Geodetic())

ax.set_global()

plt.show()

When plotting Matplotlib adjusts axes limits to the plotted objects. With GeoAxes.set_global we reset limits to show the whole map.

Interactive maps with Mplleaflet

The Python library Mplleaflet allows to show Matplotlib drawings on an interactive map based on Open Street Map or other map services. The result is a webpage, which can be embedded into an Jupyter notebook.

Use longitude and latitude values for plotting and then call mplleaflet.display to show the map inside the Jupyter notebook. Alternatively, call mplleaflet.show to open the map in new window.

fig = plt.figure(figsize=(10,6))
ax = fig.add_axes((0, 0, 1, 1))
ax.plot([12, 13, 13, 12, 12], [51, 51, 50, 50, 51], '-ob')
mplleaflet.display(fig=fig)

Exercises

1. Have a look at the documentation of GeoAxes.

2. Visualize mean traffic densities for 2018 on an interactive map of Germany. Here traffic density is number of vehicles per lane and hour. Indicate each measuring station by a colored dot with colors ranging from green (few traffic) over yellow to red (high traffic). Solution