Joseph Barraud
This notebook introduces the various possibilities of grid display with interpies
. The function show
combines several functions of matplotlib
so that relatively complex maps can be created with only one line of code.
The range of values of geophysical data, typically from gravity and magnetic surveys, is often dominated by extreme and localised anomalies. Displaying this sort of data can be a challenge. For example, using a linear scale would achieve a true representation of the data content; and large highs and lows would be clearly visible. However, smaller anomalies would also certainly disappear in a large extent of uniform colour.
Two of the most common ways to solve this problem are hillshading, and data normalisation. Hillshading works with derivatives so will naturally remove large trends and highlight small gradients. Normalisation modifies the data to rebalance the distribution of anomalies: the simplest way to do that is to clip highs and lows to reduce the range of values. There are more sophisticated ways to normalise the data, for example using a logarithmic scale. See the matplotlib documentation for more information.
The way data values are turned into colours is called colormapping. As explained in the matplotlib documentation, this involves two steps: (1) the data is normalised to the 0-1 interval, (2) the 0-1 interval is mapped to a sequence of colours (the "colormap"). In this approach, normalisation happens in the first step, and the second step is purely a linear match between two sequences.
In interpies
, I have experimented a different approach: instead of modifying the data, the colormap is modified to visually achieve the same result as in the first approach. The two steps of colormapping remain linear but the colormap has changed, for example to rebalance an uneven data distribution.
The advantage is that the colorbar works as a visual indicator of the amount of normalisation that has been applied to the data.
Let's look at a few examples of how this works in practice.
In [1]:
import sys
import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline
In [2]:
import interpies
The example data come from an aeromagnetic survey collected for the USGS over the town of Blanca in Colorado (Bankey and Grauch, 2004). The survey lines, the gridded data and the report are available from this page.
The grid can be loaded with the function open
.
In [3]:
inFile = '../data/brtpgrd.gxf'
grid1 = interpies.open(inFile)
Once loaded in interpies
as a Grid object, gridded data can be displayed with the show()
method.
In [4]:
ax = grid1.show()
The default parameters are the ones from the imshow_hs
function in the graphics
module.
imshow_hs(source, ax=None, cmap='geosoft', cmap_norm='equalize', hs=True,
zf=10, azdeg=45, altdeg=45, dx=1, dy=1, hs_contrast=1.5, cmap_brightness=1.0,
blend_mode='alpha', alpha=0.7, contours=False, colorbar=True,
cb_contours=False, cb_ticks='linear', std_range=1, figsize=(8, 8),
title=None, **kwargs)
The default normalisation is equalisation. It can have quite a pronounced effect on the result but it is generally well suited for grids with an uneven data distribution. After all, that is what the algorithm of histogram equalisation is about: redistributing values to obtain a flat histogram, therefore maximising the contrast.
In [5]:
ax = grid1.show(cmap_norm='none')
Notice how the colorbar has changed between the two plots. The second one is the way it would normally look like with the pyplot.imshow()
function. In the previous plot, the colorbar clearly indicates that a large part of the data range is represented by only one colour, so it's an honest way to remind the reader that the colours have been "manipulated" somehow to highlight some of the smaller anomalies.
In [6]:
ax = grid1.show(cmap_norm='none', hs=False)
In [7]:
ax = grid1.show(cmap='coolwarm')
Here is parula
, the MATLAB default.
In [8]:
ax = grid1.show(cmap='parula')
Let's try viridis
, the matplotlib default. Let's also remove the normalisation.
In [9]:
ax = grid1.show(cmap='viridis', cmap_norm='none')
I have also included the cet_rainbow
, which a rainbow colormap designed by Peter Kovesi. This particular version is actually called rainbow_bgyr_35_85_c73
and was downloaded from the colorcet project.
In [10]:
ax = grid1.show(cmap='cet_rainbow')
Autolevels normalisation works by clipping the data intelligently. The upper and lower boundaries are calculated as percentiles (by default the 98th and 2nd ones, respectively) so they should be well adapted to the input values. Thinking of the distribution of the data in an histogram, this normalisation effectively clips the tails of the typical bell-shaped curve that characterises most datasets. As a result, the histogram is stretched, improving the contrast of the image.
As before with equalisation, the autolevels are calculated from the data but applied to the colormap. So the colorbar on the side still shows the full range of values and makes it clear that a lot of values are represented with the same colour.
In [11]:
ax = grid1.show(cmap_norm='auto')
Let's try tighter clipping boundaries by changing the percentiles.
In [12]:
ax = grid1.show(cmap_norm='auto', min_percent=10, max_percent=90)
Here is the same with a spectral colormap and no hillshading. Weird looking but effective at isolating extremes.
In [13]:
ax = grid1.show(hs=False, cmap='nipy_spectral', cmap_norm='autolevels', min_percent=10, max_percent=90)
In [14]:
ax = grid1.show(cmap_norm='none', vmin=-600, vmax=-200)
All the different types of normalisation from matplotlib
are also available in interpies
, including BoundaryNorm
. This allows you to create a colormap composed of discrete colours.
Here is an example using boundaries that are equally spaced from the mean. The graphics.stats_boundaries
function provides directly the necessary list of values with an interval equal to a multiple of the standard deviation. The function also adds the minimum and maximum values to the list.
So here is an example in which a four-colour colormap is created: the central value is the mean (-418), then on both sides there are two intervals equal to the standard deviation, and finally two colours at each end to represent values up to the maximum and minimum of the range.
In [15]:
import matplotlib.colors as mcolors
from interpies import graphics
bounds = graphics.stats_boundaries(grid1.data, std_range=1, step=0.5)
norm = mcolors.BoundaryNorm(boundaries=bounds, ncolors=256)
ax = grid1.show(cmap='Reds', cmap_norm='none', norm=norm)
Despite the appearance of the colorbar, only the four intervals around the mean are equal, the two extreme ones are much bigger.
One can also increase the number of intervals and decrease the size of each interval. In the following cell, the cb_ticks
parameter is a trick to ensure meaningful statistics are displayed on the colorbar.
In [16]:
bounds = graphics.stats_boundaries(grid1.data, std_range=2, step=0.5)
norm = mcolors.BoundaryNorm(boundaries=bounds, ncolors=256)
ax = grid1.show(cmap='Reds', cmap_norm='none', norm=norm, cb_ticks='stats')
In [17]:
ax = grid1.show(cmap_norm='equalize', cb_ticks='stats', std_range=2)
In [18]:
ax = grid1.show(cmap='YlGnBu_r', cmap_norm='autolevels')
The cmap_brightness
parameter aims at correcting that effect by increasing the brightness of the colormap. A value greater than 1.0 will increase the brightness, while a value smaller than 1.0 will decrease it.
In [19]:
ax = grid1.show(cmap='YlGnBu_r', cmap_norm='autolevels', cmap_brightness=1.5)
Here is a darker version.
In [20]:
ax = grid1.show(cmap='YlGnBu_r', cmap_norm='autolevels', cmap_brightness=0.75)
In [21]:
ax = grid1.show(alpha=0.0)
In [22]:
ax = grid1.show(alpha=0.0, hs_contrast=0.5)
Here is an example with high contrast.
In [23]:
ax = grid1.show(cmap='inferno', cmap_norm='autolevels', min_percent=1, max_percent=99,
cmap_brightness=1.2, hs_contrast=3.0, alpha=0.6)
Hillshading has been designed originally for terrain data. So when visualising gravity and magnetic data, the Z unit is artificially assumed to be metres, as if it is topography. Therefore, appropriate scaling of Z values is generally necessary to obtain a satisfying shading. The result also depends on the cell size.
The zf
parameter (or Z factor) controls the vertical exaggeration. The default value is 10.
Here is the version with a smaller vertical exaggeration.
In [24]:
ax = grid1.show(alpha=0.0, zf=1, hs_contrast=1)
This is what happens when zf
increases.
In [25]:
ax = grid1.show(alpha=0.0, zf=1000, hs_contrast=1)
In [26]:
ax = grid1.show(contours=True)
The contours levels can also be shown on the colorbar with the cb_contours
boolean. This example shows how the contour levels are equally spaced starting at zero.
In [27]:
ax = grid1.show(contours=12, cb_contours=True)
There is also an argument (ct_colors
) to change the colour of the contours. Here I am making them slightly transparent using a RGBA tuple.
In [28]:
ax = grid1.show(cmap='Spectral_r', cmap_norm='auto', cb_ticks='stats', contours=24, ct_colors=[(0, 0, 0, 0.6)])
The keyword ct_cmap
controls the colormap of the contours, to be used together with ct_colors=None
.
In [29]:
ax = grid1.show(alpha=0.0, zf=0.5, hs_contrast=0.1, contours=True, ct_colors=None, ct_cmap='magma')
A final example with white contours and no hillshade.
In [30]:
ax = grid1.show(hs=False, cmap='seismic', cmap_norm='auto', contours=32, ct_colors='w', cb_ticks='stats')