Wed 22 February 2017

Filed under Cartopy

Tags python cartopy layout

Cartopy and Map Layout

So I was browsing an excellent blog, and came across this post The author discusses the various components you expect on a paper map (apart from just the spatial content), and proposes a number of layouts.

Now obviously, the article was intended mainly for users with a full-powered GIS. However a full-function GIS, even the freely available QGIS, is pretty large beast, and the learning curve to master it, and accomodate it into a workflow would be steep (or, at least, it has been for me, even given the comprehensive tutorial material available).

So I thought to myself, how would I implement one of the recommended layouts in Matplotlib and Cartopy. This post is the result.

Axes are your friend(s)

Excuses in advance to Matplotlib experts who read this, as I have probably got some of the key concepts askew (but that's all part of the learning). The main concept behind my solution is the concept that a figure can contain a number of Axes objects that co-exist. Each Axes object can be commanded to turn on or off the 'decorations' that surround the plotting content.

So the first thing is to define a function that will clear away all the decorations, leaving us with the equivalent of a blank space on the paper.

    def blank_axes(ax):
        blank_axes:  blank the extraneous spines and tick marks for an axes

        ax:  a matplotlib Axes object

        Output: None

        ax.tick_params(labelbottom='off', labeltop='off', labelleft='off', labelright='off' ,\
                        bottom='off', top='off', left='off', right='off' )
    #end blank_axes

The term spine is a little misleading (you or I have only one spine, but an Axes has four: top, bottom, left, right). On a spine there are tick marks, and labels for those tick marks. We clear them all away.

Draw the outer frame

    fig = plt.figure(figsize=(10,12))

    # ------------------------------- Surrounding frame ------------------------------
    # set up frame full height, full width of figure, this must be called first

    left = -0.05
    bottom = -0.05
    width = 1.1
    height = 1.05
    rect = [left,bottom,width,height]
    ax3 = plt.axes(rect)

    # turn on the spines we want, ie just the surrounding frame

    ax3.text(0.01,0.01,'© Don Cameron, 2017: '+
             'Map generated at '"%Y-%m-%d %H:%M:%S") + ' from ' +theNotebook, fontsize=8)

This outer frame is a little big larger than the usual square of side one, so the surround box doesn't touch any interior component. We set the enclosing spines visible, but supress tick marks, etc. Note that I have included a Copyright text, and some information to support reproducibility (the classic small print at the bottom of the document). I get the Jupyter Notebook name by some JavaScript hackery I got from Stack Exchange, but this could be a Python Module name, or anything you like.

Draw the spatial data

    # ---------------------------------  Main Map -------------------------------------
    # set up main map almost full height (allow room for title), right 80% of figure

    left = 0.2
    bottom = 0
    width = 0.8
    height = 0.90
    rect = [left,bottom,width,height]

    ax = plt.axes(rect, projection=ccrs.PlateCarree(), )
    ax.set_extent((150, 155, -30, -23))

    ax.coastlines(resolution='10m', zorder=2)

    #land polygons, including major islands, use cartopy default color
    ax.add_feature(BORDERS2_10m, edgecolor='grey')
    # stock image is good enough for example, but OCEAN_10m could be used, but very slow
    #       ax.add_feature(OCEAN_10m)

    ax.gridlines(draw_labels=True, xlocs=[150, 152, 154, 155])

I have allocated about 80% of the width of the map to the Cartopy spatial representation. The height is set to 90%, to allow a title at the top.

Because I want a high quality map, I have defined features using the Natural Earth 1:10Million scale data, as follows:

    BORDERS2_10m = cartopy.feature.NaturalEarthFeature('cultural', 'admin_1_states_provinces',
                                  '10m', edgecolor='black', facecolor='none')
    #"""country boundaries.""""

This gives us our state borders. We have facecolor='none', because we only want a border, not each state shaded.

    LAND_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',
    #"""land polygons, including major islands."""

Note that Cartopy defines a set of colors for land, sea, etc., which I am re-using.

    RIVERS_10m = cartopy.feature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m',
    #"""single-line drainages, including lake centerlines."""

Note that by default, Cartopy draws its tick labels outside the main map, so we have to adjust the surrounding Axes objects so they don't overlap them.

Scale Bar

Now come the hacky-est part of this post. Basemap has a nifty method to draw scale bars (drawmapscale()), but so far as I know, Cartopy has no such method. Maybe this is because Cartopy is oriented towards those who want to display quantitive or qualitative data in a spatial context, and they don't expect people to use their maps to measure distance. However, in my view, a scale bar is part of the spatial data that should be displayed. Just how big are those Texas county borders that we keep seeing in examples, color-coded for poverty, or some such sociological attribute?

The double-sad part is that I know on no easy way to draw a scale bar in one Axes object, that will exactly match the distance measurements on the main map. Maybe if I define a second Cartopy Axes, with lat/lon range adjusted for the different sizes of the Axes objects ... but that is a matter for more investigation.

What all this means is that I can't follow the original blog post that prompted this post, which suggested the scale bar as a separate object, off to the side. I am going to have to write it on the main map. Not too bad a compromise.

Define the scale bar parameters

    lon0, lon1, lat0, lat1 = ax.get_extent()

    # bar offset is how far from bottom left corner scale bar is (x,y) and how far up is scale bar text
    bar_offset = [0.05, 0.05, 0.07]
    bar_lon0 = lon0 + (lon1-lon0)*bar_offset[0]
    bar_lat0 = lat0 + (lat1-lat0)*bar_offset[1]

    text_lon0 = bar_lon0
    text_lat0 = lat0 + (lat1-lat0)*bar_offset[2]
    bar_tickmark = 20000 # metres
    bar_ticks = 5
    bar_alpha = 0.3

    bar_color = ['black', 'red']

My scale bar will be 5% up, 5% right from the lower left corner; made up of 5 sub-bars each 20 km long, alternating red and black, and quite transparent (in case they hide something important; alpha=0.3). The text will be up 7%, above the bar.

Draw the Scale Bar

    # draw a scale bar that is a set of colored line segments (bar_ticks of these), bar_tickmarks long
    for i in range(bar_ticks):
         #  90 degrees = direction of horizontal scale bar
        end_lat, end_lon = displace(bar_lat0, bar_lon0, 90, bar_tickmark)
        # capstyle must be set so line segments end square
        #TODO make transform match ax projection
        ax.plot([bar_lon0, end_lon], [bar_lat0, end_lat], color=bar_color[i%2], linewidth=20, 
                transform=ccrs.PlateCarree(), solid_capstyle='butt', alpha = bar_alpha)
        # start of next bar is end of last bar
        bar_lon0 = end_lon
        bar_lat0 = end_lat
    #end for

I have a helper function displace() that takes:

  • a lat/lon (in degrees - no love for radians here!)

  • a direction (degrees)

  • a distance in metres

and returns the lat/lon at the end of that distance in that direction. I draw a set of lines (NOT bars), but as a result, I have to set the appearence of the line segment end as 'butt', or else the line will be drawn with rounded ends that go past the point I want them to stop or start at.

Draw the scale bar text

    # highlight text with white background
    buffer = [patheffects.withStroke(linewidth=3, foreground="w")]
    # Plot the scalebar label
    units = 'km'
    #TODO make transform match ax projection
    t0 = ax.text(text_lon0, text_lat0, str(bar_ticks*bar_tickmark/1000) + ' ' + units, transform=ccrs.PlateCarree(),
        horizontalalignment='left', verticalalignment='bottom',
        path_effects=buffer, zorder=2)

We put a white background behind the black scale bar text.

Locating Map

    # ---------------------------------Locating Map ------------------------
    # set up index map 20% height, left 16% of figure
    left = 0
    bottom = 0
    width = 0.16
    height = 0.2
    rect = [left,bottom,width,height]

    ax2 = plt.axes(rect, projection=ccrs.PlateCarree(), )
    ax2.set_extent((110,160, -45, -10))
    #  ax2.set_global()  will show the whole world as context

    ax2.coastlines(resolution='110m', zorder=2)


    lon0,lon1,lat0,lat1 = ax.get_extent()
    box_x = [lon0, lon1, lon1, lon0, lon0]
    box_y = [lat0, lat0, lat1, lat1, lat0]

    plt.plot(box_x, box_y, color='red',  transform=ccrs.Geodetic())

We reserve a small area at the bottom left for a locating map. In my example, I have drawn Australia (well, some far flung islands like Christmas Island may be missing). A set_global() call will force the whole globe to be shown, but that is zoomed too far back for my case. I add coastlines, and shade the land and oceans with the default Cartopy resolution and colors. I draw gridlines, and then draw a red box that represents the main map bounds.

The Title

    # -------------------------------- Title -----------------------------
    # set up map title top 4% of figure, right 80% of figure

    left = 0.2
    bottom = 0.95
    width = 0.8
    height = 0.04
    rect = [left,bottom,width,height]
    ax6 = plt.axes(rect)
    ax6.text(0.5, 0.0,'Multi-Axes Map Example', ha='center', fontsize=20)

The title is centered above the map. I could use some ornate font, but just went with the Cartopy defaults. Note that no Cartopy projection systems are used, just the 0.0-1.0 default Axes object units. We turn off all spines, etc.

The North Arrow

    # ---------------------------------North Arrow  ----------------------------
    left = 0
    bottom = 0.2
    width = 0.16
    height = 0.2
    rect = [left,bottom,width,height]
    rect = [left,bottom,width,height]
    ax4 = plt.axes(rect)

    # need a font that support enough Unicode to draw up arrow. need space after Unicode to allow wide char to be drawm?
    ax4.text(0.5, 0.0,u'\u25B2 \nN ', ha='center', fontsize=30, family='Arial', rotation = 0)

The North arrow lives above the locating map. The only gotchas of note here are:

  • You need to specify a font family that support the Unicode character used (the default font doesn't, in my case). I am uncertain as to where Matplotlib / Cartopy gets its fonts from, so there may be scope for a more florid North arrow.

  • In my case, the Unicode needed a following space before the training newline character, to be drawn properly. Why, I don't know.

In my case, my map is NOT rotated, but the text( rotation= ) parameter would allow for a rotated arrow, if needed.

the Legend

The legend is bit hacky as well. A simple call to the Matplotlib 'make me a legend' usually fails because under the covers, the necessary calls haven't been made to link a drawing object (like a line with linestyle and linecolor), to a name. So we have to do all that.

The trick is to create some drawing entities (patches and lines), get the handles to these, associate them with names, and then create a legend. It's not all bad, it does give you the chance to control the legend in fine detail. Spoiler alert, GeoPandas legends can be quite painful.

Create the Axes Object

    # ------------------------------------  Legend -------------------------------------

    # legends can be quite long, so set near top of map (0.4 - bottom + 0.5 height = 0.9 - near top)
    left = 0
    bottom = 0.4
    width = 0.16
    height = 0.5
    rect = [left,bottom,width,height]
    rect = [left,bottom,width,height]
    ax5 = plt.axes(rect)

Create Area Legend Entries

    # create an array of color patches and associated names for drawing in a legend
    # colors are the predefined colors for cartopy features (only for example, Cartopy names are unusual)
    colors = sorted(cartopy.feature.COLORS.keys())

    # handles is a list of patch handles
    handles = []
    # names is the list of corresponding labels to appear in the legend
    names = []

    # for each cartopy defined color, draw a patch, append handle to list, and append color name to names list
    for c in colors:
        patch = mpatches.Patch(color=cfeature.COLORS[c], label=c)
    #end for

Note that these Patches never appear in our diagram. Now Cartopy color names are not perfect, but it's only an example. If at some time Cartopy expend their color range, the code above should still work.

Create Line Legend Entries

    # do some example lines with colors
    river = mlines.Line2D([], [], color=cfeature.COLORS['water'], marker='',
                              markersize=15, label='river')
    coast = mlines.Line2D([], [], color='black', marker='',
                              markersize=15, label='coast')
    bdy  = mlines.Line2D([], [], color='grey', marker='',
                      markersize=15, label='state boundary')
    names.append('state boundary')

Note that here we have turned off all markers, and have invented our own labels to appear in the legend.

Create the Legend, and Display All!

    # create legend
    ax5.legend(handles, names)

It looks like by default, the legend() call turns the spines back on, but that is OK. We specify a legend title.

The Results

The end result of all this? See below.

Map Layout Results

In order to make sense of the different Axes objects, here they are again, with no content, and with the enclosing spines drawn. Some of these (like the main map) will expand as they are filled with content.

Map Layout Results

The Imports

Shown below are the imports for the Notebook from which this example was taken. Not all are used in the example fragments shown above.

    # all imports should go here

    import pandas as pd
    import sys
    import os
    import subprocess
    import datetime
    import platform
    import datetime
    import math

    import matplotlib.pyplot as plt
    #import seaborn as sb

    import cartopy
    import as ccrs
    from import OSM
    import cartopy.feature as cfeature
    from import shapereader
    from import StamenTerrain
    from import GoogleTiles
    from owslib.wmts import WebMapTileService

    from matplotlib.path import Path
    import matplotlib.patheffects as PathEffects
    from matplotlib import patheffects
    import matplotlib.patches as mpatches
    import matplotlib.lines as mlines

    import numpy as np
Comment Data Analysis Blog © Don Cameron Powered by Pelican and Twitter Bootstrap. Icons by Font Awesome and Font Awesome More