Colormap example for bathymetry
Manipulating colormaps for bathymetric data
The following notebook gives one example of how to make a decent looking chart with filled bathymetric contours using matplotlib - decent being in the eye of the beholder. This builds on discussion from the 27 Oct 2014 PyHOGS meeting, earlier meetings on colormaps with matplotlib, and on the 10 Nov 2014 PyHOGS meeting.
Aim
Using Matlab it's hard to directly align desired contour intervals with the desired color, and this becomes very tricky if your contour intervals are not evenly spaced. It turns out there are easy ways to do this with matplotlib.
Colormap
I'd like a colormap that goes from white (shallow water) to dark blue (deep water), that has more colors for shallow depths for resolving large bathymetric changes close to land, and that has a small discrete number of colors. Few discrete colors helps the viewer to determine which isobath corresponds to what depth.
Nonuniform contour intervals
There is an easy way to implement this with built-in scaling functions that scale the data to the colormap. Although one can always scale the data oneself before plotting, e.g. by taking its log, it is more convenient in python to use the built in normalization methods instead.
Generally speaking, this is related to how the colormap is mapped to the data.
Normally (at left), this is done automatically, or by matching the [0,1] range of the colormap to the specified [vmin,vmax] range of the data (here, I've used the [-8000,0] range of the bathymetric dataset I plot below). Uniform intervals in the data are mapped to uniform intervals in the colormap. What I'm trying to do (at right) is to align non-uniform levels in the data to uniform levels in the colormap, to compress/extend the colormap to better resolve certain ranges of data.
Spacing of colorbar ticks
When using plt.colorbar(), the contour intervals can be shown two ways, either with uniform spacing on the plotted colorbar (the default)
colorbar(pc, spacing='uniform')
or by placing the colorbar ticks relative to their respective numerical value
colorbar(spacing='proportional')
Implementation
Import and define functions
%matplotlib inline
from __future__ import division # to remove integer division (e.g. 2/5=0)
from __future__ import unicode_literals # to remove unicode/ascii distinction in python2
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy.io as spio
import matplotlib
# Use custom colormap function from Earle
def custom_div_cmap(numcolors=11, name='custom_div_cmap',
mincol='blue', midcol='white', maxcol='red'):
""" Create a custom diverging colormap with three colors
Default is blue to white to red with 11 colors. Colors can be specified
in any way understandable by matplotlib.colors.ColorConverter.to_rgb()
"""
from matplotlib.colors import LinearSegmentedColormap
cmap = LinearSegmentedColormap.from_list(name=name,
colors =[mincol, midcol, maxcol],
N=numcolors)
return cmap
# make negative contours, normally dashed by default, be solid
matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
Load in bathymetric data from the western subtropical Pacific, off Luzon Island
# load bathymetry data from matlab data file
BB = spio.loadmat('/Users/zszuts/python/bathy.mat');
blon = BB['blon']
blat = BB['blat']
bdepth = BB['bdepth']
Blon,Blat = np.meshgrid(blon,blat)
Setup contour intervals and colormaps
blevels = [-8000, -7000, -6000, -5000, -4000, -3000, -2000, -1500, -1000, -500, -200, -100, 0]
N = len(blevels)-1
cmap2 = custom_div_cmap(N, mincol='DarkBlue', midcol='CornflowerBlue' ,maxcol='w')
cmap2.set_over('0.7') # set positive values (land) as light gray
Simple attempt: the poor result of basic plotting
Now start plotting. First do a naive plot to show why the straightforward approach doesn't give a good looking example.
fig, ax = plt.subplots(1,1, figsize=(10,6))
ax.set_aspect(1/np.cos(np.average(blat)*np.pi/180)) # set the aspect ratio for a local cartesian grid
pc = plt.contourf(Blon,Blat,bdepth, vmin=-8000, vmax=0, levels=blevels, cmap=cmap2, extend='both')
plt.colorbar(pc, ticks=blevels, spacing='uniform')
Comments on number of colors in the colormap
My preference is for a small number of colors in the colormap, since this facilitates determining the value each contour stands for without requiring inline labels. When a colormap is used that has a small number of discrete levels, multiple isobaths can end up with the same color, e.g. the example above has white for all contours shallower than -500 m.
The linear mapping from the data range [vmin,vmax] to the colormap range [0,1] gives a precise decimal value. The associated color is found by finding the nearest available color. If there are many levels in the colormap (effectively 30-40 or higher, I find), then the color is interpolated well and the colormap appears (quasi) continuous. If there are only a few levels in the colormap, then the nearest one is found even if it is not unique.
This is related to the need for divergent color maps (for instance blue-white-red) to have an odd number of colors so that the central color white is explicitly resolved - see Colormap Examples for more discussion of this.
Second attempt: use normalization to meet my aims
Now make a non-linear colormap using BoundaryNorm() and make it pretty. I make two colorbars to show the difference between uniform and proportional spacing. The proportionally-spaced colorbar requires the removal of one tick level to avoid overprinting of the tick labels.
blevels = [-8000, -7000, -6000, -5000, -4000, -3000, -2000, -1500, -1000, -500, -200, -100, 0]
N = len(blevels)-1
cmap2 = custom_div_cmap(N, mincol='DarkBlue', midcol='CornflowerBlue' ,maxcol='w')
cmap2.set_over('0.7') # light gray
from matplotlib.colors import BoundaryNorm
bnorm = BoundaryNorm(blevels, ncolors=N, clip=False)
fig, ax = plt.subplots(1,1, figsize=(10,6))
ax.set_aspect(1/np.cos(np.average(blat)*np.pi/180)) # set the aspect ratio for a local cartesian grid
pc = plt.contourf(Blon,Blat,bdepth, vmin=-8000, vmax=0, levels=blevels, norm=bnorm, cmap=cmap2, extend='both')
plt.colorbar(pc, ticks=blevels[:-2]+[blevels[-1]], spacing='proportional')
plt.colorbar(pc, ticks=blevels, spacing='uniform')
pc = plt.contour(Blon,Blat,bdepth, levels=[-1000], colors='k')
Another way to plot the land is to extract the zero isobath and plot it as a fill object. This provides more control over the linestyle at the expense of more complex code. (For some reason, if I use the pc created by the first call to plt.contour with levels=blevels, the zero contours are not continuous. This results in the multiple fill objects, such that the center of the big island of Luzon is not filled as desired.)
fig, ax = plt.subplots(1,1, figsize=(10,6))
ax.set_aspect(1/np.cos(np.average(blat)*np.pi/180)) # set the aspect ratio for a local cartesian grid
pc = plt.contourf(Blon,Blat,bdepth, vmin=-8000, vmax=0, levels=blevels, norm=bnorm, cmap=cmap2, extend='min')
plt.colorbar(pc, ticks=blevels, spacing='uniform')
# extract the coastlines and plot separately with lines around them
pc = plt.contour(Blon,Blat,bdepth, levels=[0], colors='0.7', linewidth=0)
b0s = pc.collections[0].get_paths()
for i in range(len(b0s)):
b0x = b0s[i].vertices[:,0]
b0y = b0s[i].vertices[:,1]
plt.fill(b0x,b0y, '0.7', linestyle='solid') # plot land as filled gray
# plot specific isobath(s)
pc = plt.contour(Blon,Blat,bdepth, levels=[-1000], colors='k')
An easier way with almost identical outcome would be to add another contour object ontop of the existing plot to define the linestyle, while using the set_over switch as above. The only difference/benefit to the use of plt.fill as above is avoiding an upward arrow on the colorbar, such as made by defining set_over.
Additional scaling options
There are other ways to automatically scale data when plotting with colormaps in python. Additional norms are available by setting the options
norm=LogNorm()
norm=PowerNorm()
norm=SymLogNorm()
These scalings need different types of input, as can be found in their documentation, for LogNorm, PowerNorm, and SymLogNorm. Symlognorm is a positive and negative lognorm connected by a linear scale across zero, which performs something trivially that is otherwise very difficult to do (I know from experience)! These scalings are not suitable for plotting bathymetry, however, since bathymetry data is negative and spans 0 (and since we don't care to emphasize topography on land).
For all of these complex normalizations, the option clip=True/False determines if data values outside the range of [vmin,vmax] is clipped to 0 or 1, respectively. Note that the documentation advises:
"Clipping silently defeats the purpose of setting the over, under, and masked colors in the colormap, so it is likely to lead to surprises; therefore the default is clip = False."