How can I show a km ruler on a cartopy / matplotlib plot? How can I show a km ruler on a cartopy / matplotlib plot? python python

How can I show a km ruler on a cartopy / matplotlib plot?


With the addition of the geodesic module in CartoPy 0.15, we can now fairly easily compute exact lengths on a map. It was a bit tricky to figure out how to find two points on a straight-on-a-map line which are the right distance-on-a-sphere apart. Once the direction on the map is specified, I perform an exponential search to find a point far enough away. I then perform a binary search to find a point close enough to the desired distance.

The scale_bar function is simple enough, but it has a lot of options. The basic signature is scale_bar(ax, location, length). ax is any CartoPy axes, location is the position of the left-side of the bar in axes coordinates (so each coordinate is from 0 to 1), and length is the length of the bar in kilometres. Other lengths are supported like with the metres_per_unit and unit_name keyword arguments.

Extra keyword arguments (like color) are simply passed to text and plot. However, keyword arguments specific to text or plot (like family or path_effects) must be passed in as dictionaries through text_kwargs and plot_kwargs.

I've included examples of what I think are the common use cases.

Please share any questions, comments, or criticisms.

scalebar.py

import numpy as npimport cartopy.crs as ccrsimport cartopy.geodesic as cgeodef _axes_to_lonlat(ax, coords):    """(lon, lat) from axes coordinates."""    display = ax.transAxes.transform(coords)    data = ax.transData.inverted().transform(display)    lonlat = ccrs.PlateCarree().transform_point(*data, ax.projection)    return lonlatdef _upper_bound(start, direction, distance, dist_func):    """A point farther than distance from start, in the given direction.    It doesn't matter which coordinate system start is given in, as long    as dist_func takes points in that coordinate system.    Args:        start:     Starting point for the line.        direction  Nonzero (2, 1)-shaped array, a direction vector.        distance:  Positive distance to go past.        dist_func: A two-argument function which returns distance.    Returns:        Coordinates of a point (a (2, 1)-shaped NumPy array).    """    if distance <= 0:        raise ValueError(f"Minimum distance is not positive: {distance}")    if np.linalg.norm(direction) == 0:        raise ValueError("Direction vector must not be zero.")    # Exponential search until the distance between start and end is    # greater than the given limit.    length = 0.1    end = start + length * direction    while dist_func(start, end) < distance:        length *= 2        end = start + length * direction    return enddef _distance_along_line(start, end, distance, dist_func, tol):    """Point at a distance from start on the segment  from start to end.    It doesn't matter which coordinate system start is given in, as long    as dist_func takes points in that coordinate system.    Args:        start:     Starting point for the line.        end:       Outer bound on point's location.        distance:  Positive distance to travel.        dist_func: Two-argument function which returns distance.        tol:       Relative error in distance to allow.    Returns:        Coordinates of a point (a (2, 1)-shaped NumPy array).    """    initial_distance = dist_func(start, end)    if initial_distance < distance:        raise ValueError(f"End is closer to start ({initial_distance}) than "                         f"given distance ({distance}).")    if tol <= 0:        raise ValueError(f"Tolerance is not positive: {tol}")    # Binary search for a point at the given distance.    left = start    right = end    while not np.isclose(dist_func(start, right), distance, rtol=tol):        midpoint = (left + right) / 2        # If midpoint is too close, search in second half.        if dist_func(start, midpoint) < distance:            left = midpoint        # Otherwise the midpoint is too far, so search in first half.        else:            right = midpoint    return rightdef _point_along_line(ax, start, distance, angle=0, tol=0.01):    """Point at a given distance from start at a given angle.    Args:        ax:       CartoPy axes.        start:    Starting point for the line in axes coordinates.        distance: Positive physical distance to travel.        angle:    Anti-clockwise angle for the bar, in radians. Default: 0        tol:      Relative error in distance to allow. Default: 0.01    Returns:        Coordinates of a point (a (2, 1)-shaped NumPy array).    """    # Direction vector of the line in axes coordinates.    direction = np.array([np.cos(angle), np.sin(angle)])    geodesic = cgeo.Geodesic()    # Physical distance between points.    def dist_func(a_axes, b_axes):        a_phys = _axes_to_lonlat(ax, a_axes)        b_phys = _axes_to_lonlat(ax, b_axes)        # Geodesic().inverse returns a NumPy MemoryView like [[distance,        # start azimuth, end azimuth]].        return geodesic.inverse(a_phys, b_phys).base[0, 0]    end = _upper_bound(start, direction, distance, dist_func)    return _distance_along_line(start, end, distance, dist_func, tol)def scale_bar(ax, location, length, metres_per_unit=1000, unit_name='km',              tol=0.01, angle=0, color='black', linewidth=3, text_offset=0.005,              ha='center', va='bottom', plot_kwargs=None, text_kwargs=None,              **kwargs):    """Add a scale bar to CartoPy axes.    For angles between 0 and 90 the text and line may be plotted at    slightly different angles for unknown reasons. To work around this,    override the 'rotation' keyword argument with text_kwargs.    Args:        ax:              CartoPy axes.        location:        Position of left-side of bar in axes coordinates.        length:          Geodesic length of the scale bar.        metres_per_unit: Number of metres in the given unit. Default: 1000        unit_name:       Name of the given unit. Default: 'km'        tol:             Allowed relative error in length of bar. Default: 0.01        angle:           Anti-clockwise rotation of the bar.        color:           Color of the bar and text. Default: 'black'        linewidth:       Same argument as for plot.        text_offset:     Perpendicular offset for text in axes coordinates.                         Default: 0.005        ha:              Horizontal alignment. Default: 'center'        va:              Vertical alignment. Default: 'bottom'        **plot_kwargs:   Keyword arguments for plot, overridden by **kwargs.        **text_kwargs:   Keyword arguments for text, overridden by **kwargs.        **kwargs:        Keyword arguments for both plot and text.    """    # Setup kwargs, update plot_kwargs and text_kwargs.    if plot_kwargs is None:        plot_kwargs = {}    if text_kwargs is None:        text_kwargs = {}    plot_kwargs = {'linewidth': linewidth, 'color': color, **plot_kwargs,                   **kwargs}    text_kwargs = {'ha': ha, 'va': va, 'rotation': angle, 'color': color,                   **text_kwargs, **kwargs}    # Convert all units and types.    location = np.asarray(location)  # For vector addition.    length_metres = length * metres_per_unit    angle_rad = angle * np.pi / 180    # End-point of bar.    end = _point_along_line(ax, location, length_metres, angle=angle_rad,                            tol=tol)    # Coordinates are currently in axes coordinates, so use transAxes to    # put into data coordinates. *zip(a, b) produces a list of x-coords,    # then a list of y-coords.    ax.plot(*zip(location, end), transform=ax.transAxes, **plot_kwargs)    # Push text away from bar in the perpendicular direction.    midpoint = (location + end) / 2    offset = text_offset * np.array([-np.sin(angle_rad), np.cos(angle_rad)])    text_location = midpoint + offset    # 'rotation' keyword argument is in text_kwargs.    ax.text(*text_location, f"{length} {unit_name}", rotation_mode='anchor',            transform=ax.transAxes, **text_kwargs)

demo.py

import cartopy.crs as ccrsimport cartopy.feature as cfeatureimport matplotlib.pyplot as pltfrom scalebar import scale_barfig = plt.figure(1, figsize=(10, 10))ax = fig.add_subplot(111, projection=ccrs.Mercator())ax.set_extent([-180, 180, -85, 85])ax.coastlines(facecolor='black')ax.add_feature(cfeature.LAND)# Standard 6,000 km scale bar.scale_bar(ax, (0.65, 0.4), 6_000)# Length of the bar reflects its position on the map.scale_bar(ax, (0.55, 0.7), 6_000, color='green')# Bar can be placed at any angle. Any units can be used.scale_bar(ax, (0.4, 0.4), 3_000, metres_per_unit=1609, angle=-90,          unit_name='mi', color='red')# Text and line can be styled separately. Keywords are simply passed to# text or plot.text_kwargs = dict(family='serif', size='xx-large', color='red')plot_kwargs = dict(linestyle='dashed', color='blue')scale_bar(ax, (0.05, 0.3), 6_000, text_kwargs=text_kwargs,          plot_kwargs=plot_kwargs)# Angles between 0 and 90 may result in the text and line plotted at# slightly different angles for an unknown reason.scale_bar(ax, (0.45, 0.15), 5_000, color='purple', angle=45, text_offset=0)# To get around this override the text's angle and fiddle manually.scale_bar(ax, (0.55, 0.15), 5_000, color='orange', angle=45, text_offset=0,          text_kwargs={'rotation': 41})plt.show()

Several different scale bars on one world map.


Here's the Cartopy scale bar function I wrote for my own use which uses simpler version of pp-mo's answer:Edit: modified code to create a new projection centred so that the scale bar is parallel to the axies for many coordinate systems, including some orthographic and larger maps, and removing the need to specify a utm system.Also added code to calculate a scale bar length if one wasn't specified.

import cartopy.crs as ccrsimport numpy as npdef scale_bar(ax, length=None, location=(0.5, 0.05), linewidth=3):    """    ax is the axes to draw the scalebar on.    length is the length of the scalebar in km.    location is center of the scalebar in axis coordinates.    (ie. 0.5 is the middle of the plot)    linewidth is the thickness of the scalebar.    """    #Get the limits of the axis in lat long    llx0, llx1, lly0, lly1 = ax.get_extent(ccrs.PlateCarree())    #Make tmc horizontally centred on the middle of the map,    #vertically at scale bar location    sbllx = (llx1 + llx0) / 2    sblly = lly0 + (lly1 - lly0) * location[1]    tmc = ccrs.TransverseMercator(sbllx, sblly)    #Get the extent of the plotted area in coordinates in metres    x0, x1, y0, y1 = ax.get_extent(tmc)    #Turn the specified scalebar location into coordinates in metres    sbx = x0 + (x1 - x0) * location[0]    sby = y0 + (y1 - y0) * location[1]    #Calculate a scale bar length if none has been given    #(Theres probably a more pythonic way of rounding the number but this works)    if not length:         length = (x1 - x0) / 5000 #in km        ndim = int(np.floor(np.log10(length))) #number of digits in number        length = round(length, -ndim) #round to 1sf        #Returns numbers starting with the list        def scale_number(x):            if str(x)[0] in ['1', '2', '5']: return int(x)                    else: return scale_number(x - 10 ** ndim)        length = scale_number(length)     #Generate the x coordinate for the ends of the scalebar    bar_xs = [sbx - length * 500, sbx + length * 500]    #Plot the scalebar    ax.plot(bar_xs, [sby, sby], transform=tmc, color='k', linewidth=linewidth)    #Plot the scalebar label    ax.text(sbx, sby, str(length) + ' km', transform=tmc,            horizontalalignment='center', verticalalignment='bottom')

It has some limitations but is relatively simple so I hope you could see how to change it if you want something different.

Example usage:

import matplotlib.pyplot as pltax = plt.axes(projection=ccrs.Mercator())plt.title('Cyprus')ax.set_extent([31, 35.5, 34, 36], ccrs.Geodetic())ax.coastlines(resolution='10m')scale_bar(ax, 100)plt.show()

enter image description here


Here's a refined version of @Siyh's answer which adds:

  • automatic UTM zone selection
  • a buffer behind the text/bar so it shows against the background
  • a North arrow

Notes:

  • if you don't use UTM for you axes, the bar will be drawn crooked
  • the North arrow assumes north is up

Code:

import osimport cartopy.crs as ccrsfrom math import floorimport matplotlib.pyplot as pltfrom matplotlib import patheffectsimport matplotlibif os.name == 'nt':    matplotlib.rc('font', family='Arial')else:  # might need tweaking, must support black triangle for N arrow    matplotlib.rc('font', family='DejaVu Sans')def utm_from_lon(lon):    """    utm_from_lon - UTM zone for a longitude    Not right for some polar regions (Norway, Svalbard, Antartica)    :param float lon: longitude    :return: UTM zone number    :rtype: int    """    return floor( ( lon + 180 ) / 6) + 1def scale_bar(ax, proj, length, location=(0.5, 0.05), linewidth=3,              units='km', m_per_unit=1000):    """    http://stackoverflow.com/a/35705477/1072212    ax is the axes to draw the scalebar on.    proj is the projection the axes are in    location is center of the scalebar in axis coordinates ie. 0.5 is the middle of the plot    length is the length of the scalebar in km.    linewidth is the thickness of the scalebar.    units is the name of the unit    m_per_unit is the number of meters in a unit    """    # find lat/lon center to find best UTM zone    x0, x1, y0, y1 = ax.get_extent(proj.as_geodetic())    # Projection in metres    utm = ccrs.UTM(utm_from_lon((x0+x1)/2))    # Get the extent of the plotted area in coordinates in metres    x0, x1, y0, y1 = ax.get_extent(utm)    # Turn the specified scalebar location into coordinates in metres    sbcx, sbcy = x0 + (x1 - x0) * location[0], y0 + (y1 - y0) * location[1]    # Generate the x coordinate for the ends of the scalebar    bar_xs = [sbcx - length * m_per_unit/2, sbcx + length * m_per_unit/2]    # buffer for scalebar    buffer = [patheffects.withStroke(linewidth=5, foreground="w")]    # Plot the scalebar with buffer    ax.plot(bar_xs, [sbcy, sbcy], transform=utm, color='k',        linewidth=linewidth, path_effects=buffer)    # buffer for text    buffer = [patheffects.withStroke(linewidth=3, foreground="w")]    # Plot the scalebar label    t0 = ax.text(sbcx, sbcy, str(length) + ' ' + units, transform=utm,        horizontalalignment='center', verticalalignment='bottom',        path_effects=buffer, zorder=2)    left = x0+(x1-x0)*0.05    # Plot the N arrow    t1 = ax.text(left, sbcy, u'\u25B2\nN', transform=utm,        horizontalalignment='center', verticalalignment='bottom',        path_effects=buffer, zorder=2)    # Plot the scalebar without buffer, in case covered by text buffer    ax.plot(bar_xs, [sbcy, sbcy], transform=utm, color='k',        linewidth=linewidth, zorder=3)if __name__ == '__main__':    ax = plt.axes(projection=ccrs.Mercator())    plt.title('Cyprus')    ax.set_extent([31, 35.5, 34, 36], ccrs.Geodetic())    ax.stock_img()    ax.coastlines(resolution='10m')    scale_bar(ax, ccrs.Mercator(), 100)  # 100 km scale bar    # or to use m instead of km    # scale_bar(ax, ccrs.Mercator(), 100000, m_per_unit=1, units='m')    # or to use miles instead of km    # scale_bar(ax, ccrs.Mercator(), 60, m_per_unit=1609.34, units='miles')    plt.show()

Demo image, map of Cyprus with scalebar