Introduction to imshow_hs

This notebook introduces the parameters of the main function imshow_hs of the graphics module.

Author: Joseph Barraud (@jobar8)

After importing a few modules, we will start by loading some magnetic data in a 2D array. Then we will display the data using a modified version of the pyplot.imshow function.

In [2]:
# import numpy and matplotlib
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
# Import local module
import graphics
In [3]:
# Check matplotlib version
matplotlib.__version__
Out[3]:
'2.0.0'

Loading data

OK, let's load some data.

I have downloaded some aeromagnetic data from the USGS website. The survey was flown in an area in southwest New Mexico. It is part of a larger set of surveys that cover the entire state. A description of the data can be found here.

The original file is in GXF format and can be found here.

I have already done the conversion to a numpy file, so loading it is as simple as:

In [4]:
# load magnetic data (USGS survey data)
mag = np.load(r'.\data\mag4017.npy')

Standard imshow

Before demonstrating how the modified version works, let's see first how it looks with the standard version of imshow.

In [17]:
# Note: the colormap is set to 'jet' to simulate the behaviour of pyplot prior to version 2.0.0
plt.imshow(mag,cmap='jet')
Out[17]:
<matplotlib.image.AxesImage at 0x240bfb452b0>

That does not look very good. It is not only because the pseudocolors are drawn with the infamous "jet" rainbow-like colormap, but also because the amplitude of the anomalies are unevenly distributed: a couple of anomalies dominate the map and the rest is displayed with just a few shades of yellow, cyan and green.

The other problem is that the map is upside-down as the original format puts the origin in the lower-left corner. This can be easily fixed by setting the extent of the grid with real-world coordinates. Let's try that with the one of the new matplotlib colormaps, "viridis" (this is now the default in matplotlib 2.0).

In [6]:
plt.imshow(mag,cmap='viridis',extent=(-308000.0, -183500.0, 0.0, 61500.0),origin='lower')
Out[6]:
<matplotlib.image.AxesImage at 0x240bd240668>

Modified colorbar with imshow_hs

OK. Time to try imshow_hs. It works like imshow but it has additional arguments, some of them are actually just calling other pyplot functions, like pyplot.colorbar. The plot can be added to an existing figure if a reference to existing axes is given. For example:

In [8]:
fig,ax = plt.subplots(figsize=(12,6))
graphics.imshow_hs(mag,ax,cmap='viridis',cmap_norm='nonorm',hs=False,colorbar=True,
                   extent=(-308000.0, -183500.0, 0.0, 61500.0),origin='lower')

There is one option for the colorbar that changes the position of the ticks and the labels. Instead of having them at regular intervals, they can be placed according to basic descriptive statistics: min, max, mean and standard deviation (sigma). To have ticks at mean+2*mean and mean-2*sigma, use:

In [9]:
fig,ax = plt.subplots(figsize=(12,6))
graphics.imshow_hs(mag,ax,cmap='viridis',cmap_norm='nonorm',hs=False,colorbar=True,
                   cb_ticks='stats',nSigma=2,
                   extent=(-308000.0, -183500.0, 0.0, 61500.0),origin='lower')

Histogram Equalization

The next option that imshow_hs introduces is the normalization of the colormap (norm_cmap). Normalization is already an option of imshow, and it remains one in imshow_hs. Its aim is to modify the distribution of intensities in the data in order to change the contrast or the brightness of the pseudocolored image (or exposure, using a photographic term that is also the name of a submodule in scikit-image that was the inspiration for this function).

Normalization of the colormap modifies the colormap instead of the data. It makes the changes clearly visible on the colorbar. There are two choices:

  • Histogram equalization (called with 'equalize' or 'equalization')
  • Autolevels (called with 'auto' or 'autolevels')

The equalization stretches the histogram and spreads it more evenly, boosting the intensity of some levels and reducing the intensity of some others. The effect may look dramatic but it can be useful to detect small variations in the data.

In [10]:
fig,ax = plt.subplots(figsize=(12,6))
graphics.imshow_hs(mag,ax,cmap='viridis',cmap_norm='equalize',hs=False,colorbar=True,
                   cb_ticks='stats',nSigma=2,
                   extent=(-308000.0, -183500.0, 0.0, 61500.0),origin='lower')

Another example of histogram equalization is shown below with the 'geosoft' colormap. This colormap (clra is the exact name) is popular in the community of geophysicists who interpret potential field data (from gravity and magnetic surveys essentially).

In [11]:
fig,ax = plt.subplots(figsize=(12,6))
graphics.imshow_hs(mag,ax,cmap='geosoft',cmap_norm='equalize',hs=False,colorbar=True,
                   cb_ticks='stats',nSigma=2,
                   extent=(-308000.0, -183500.0, 0.0, 61500.0),origin='lower')

The other method of normalization that is implemented in imshow_hs is 'autolevels'. It is actually a simple intensity rescaling with clipping of the extremes. It has two additional parameters, minPercent and maxPercent. It makes the colormap look brighter at high amplitude and darker for the lowest amplitudes.

In [12]:
fig,ax = plt.subplots(figsize=(12,6))
graphics.imshow_hs(mag,ax,cmap='viridis',cmap_norm='autolevels',hs=False,colorbar=True,
                   cb_ticks='stats',nSigma=2,minPercent=2,maxPercent=98,
                   extent=(-308000.0, -183500.0, 0.0, 61500.0),origin='lower')

Contours

Contour lines can be added with the contours option (True or False). The interval between the lines is automatically chosen if the levels argument is an integer. In this case, the number is the number of levels drawn on the map. If you want more control on the intervals, give it a sequence of increasing values. Like so:

In [13]:
fig,ax = plt.subplots(figsize=(12,6))
graphics.imshow_hs(mag,ax,cmap='coolwarm',cmap_norm='autolevels',hs=False,colorbar=True,
                   cb_ticks='stats',nSigma=2,minPercent=2,maxPercent=98,
                   contours=True,levels=[-600,-300,0,300,600],
                   extent=(-308000.0, -183500.0, 0.0, 61500.0),origin='lower')