In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import Image
import pandas as pd
from matplotlib import cm


import cartopy
import cartopy.crs as ccrs
from cartopy import config
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy import feature as cfeature
from cartopy.feature import NaturalEarthFeature, LAND, COASTLINE, OCEAN, LAKES, BORDERS
import matplotlib.ticker as mticker
import warnings 
# this will allow us to use the code in peace :) 
warnings.filterwarnings("ignore")

Lecture 19:

  • Learn about gridding and contouring with cartopy

Gridding

Remember this figure?

In [2]:
Image(filename='Figures/etopo20.jpg',width=500)
Out[2]:

We now have the tools to make it!

First, you should recognize the basic orthographic projection from last lecture.

In [3]:
plt.figure(1,(6,6))
ax = plt.axes(projection=ccrs.Orthographic(-75, 42))
San_lat=33
San_lon=-117%360  # takes the west longitude and converts to 0=>360
gl=ax.gridlines(crs=ccrs.PlateCarree(),color='black',linewidth=1,linestyle='dotted')
#gl.xlabels_top = False
gl.ylocator=mticker.FixedLocator(np.arange(-80,81,20))
gl.xlocator=mticker.FixedLocator(np.arange(-180,181,60));
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
ax.set_global()
ax.coastlines();

The topographic data we already encountered in the lecture on the hypsographic curve. Let's import them again here and look at their shapes again:

In [4]:
etopo=np.loadtxt('Datasets/Etopo/etopo20data.gz')
elons=np.loadtxt('Datasets/Etopo/etopo20lons.gz')
elats=np.loadtxt('Datasets/Etopo/etopo20lats.gz')
print (etopo.shape)
print (elons.shape)
print (elats.shape)
(540, 1081)
(1081,)
(540,)

elats is a 1D array with 540 latitudinal bins (this is a long and skinny array- 540 x 1).

elons is a 1D array with 1081 longitudinal bins (this as a fat and wide array- 1 x 1081).
And etopo is a 2D array with 540 rows and 1081 columns (540 x 1080).

So etopo has an elevation for each lat/lon cell.

In order to plot the elevation data onto a lat/lon grid, we have to first make a grid out of the 1D arrays of elats/ and elons.

We use the numpy function meshgrid for this.

In [5]:
help (np.meshgrid)
Help on function meshgrid in module numpy.lib.function_base:

meshgrid(*xi, **kwargs)
    Return coordinate matrices from coordinate vectors.
    
    Make N-D coordinate arrays for vectorized evaluations of
    N-D scalar/vector fields over N-D grids, given
    one-dimensional coordinate arrays x1, x2,..., xn.
    
    .. versionchanged:: 1.9
       1-D and 0-D cases are allowed.
    
    Parameters
    ----------
    x1, x2,..., xn : array_like
        1-D arrays representing the coordinates of a grid.
    indexing : {'xy', 'ij'}, optional
        Cartesian ('xy', default) or matrix ('ij') indexing of output.
        See Notes for more details.
    
        .. versionadded:: 1.7.0
    sparse : bool, optional
        If True a sparse grid is returned in order to conserve memory.
        Default is False.
    
        .. versionadded:: 1.7.0
    copy : bool, optional
        If False, a view into the original arrays are returned in order to
        conserve memory.  Default is True.  Please note that
        ``sparse=False, copy=False`` will likely return non-contiguous
        arrays.  Furthermore, more than one element of a broadcast array
        may refer to a single memory location.  If you need to write to the
        arrays, make copies first.
    
        .. versionadded:: 1.7.0
    
    Returns
    -------
    X1, X2,..., XN : ndarray
        For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` ,
        return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij'
        or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy'
        with the elements of `xi` repeated to fill the matrix along
        the first dimension for `x1`, the second for `x2` and so on.
    
    Notes
    -----
    This function supports both indexing conventions through the indexing
    keyword argument.  Giving the string 'ij' returns a meshgrid with
    matrix indexing, while 'xy' returns a meshgrid with Cartesian indexing.
    In the 2-D case with inputs of length M and N, the outputs are of shape
    (N, M) for 'xy' indexing and (M, N) for 'ij' indexing.  In the 3-D case
    with inputs of length M, N and P, outputs are of shape (N, M, P) for
    'xy' indexing and (M, N, P) for 'ij' indexing.  The difference is
    illustrated by the following code snippet::
    
        xv, yv = np.meshgrid(x, y, sparse=False, indexing='ij')
        for i in range(nx):
            for j in range(ny):
                # treat xv[i,j], yv[i,j]
    
        xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
        for i in range(nx):
            for j in range(ny):
                # treat xv[j,i], yv[j,i]
    
    In the 1-D and 0-D case, the indexing and sparse keywords have no effect.
    
    See Also
    --------
    index_tricks.mgrid : Construct a multi-dimensional "meshgrid"
                     using indexing notation.
    index_tricks.ogrid : Construct an open multi-dimensional "meshgrid"
                     using indexing notation.
    
    Examples
    --------
    >>> nx, ny = (3, 2)
    >>> x = np.linspace(0, 1, nx)
    >>> y = np.linspace(0, 1, ny)
    >>> xv, yv = np.meshgrid(x, y)
    >>> xv
    array([[ 0. ,  0.5,  1. ],
           [ 0. ,  0.5,  1. ]])
    >>> yv
    array([[ 0.,  0.,  0.],
           [ 1.,  1.,  1.]])
    >>> xv, yv = np.meshgrid(x, y, sparse=True)  # make sparse output arrays
    >>> xv
    array([[ 0. ,  0.5,  1. ]])
    >>> yv
    array([[ 0.],
           [ 1.]])
    
    `meshgrid` is very useful to evaluate functions on a grid.
    
    >>> x = np.arange(-5, 5, 0.1)
    >>> y = np.arange(-5, 5, 0.1)
    >>> xx, yy = np.meshgrid(x, y, sparse=True)
    >>> z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2)
    >>> h = plt.contourf(x,y,z)

Let's try out meshgrid on some smaller arrays first, to wrap our mind around it, so to speak.

In [6]:
x = np.arange(-3.2, 3.3, 0.1) # make a 1D array from -3.3 to 3.2 with a spacing of .1
y = np.arange(-3.2, 3.3, 0.1) # ditto
xx, yy = np.meshgrid(x, y) # make a meshgrid
print (x.shape)
print (y.shape)
print (xx.shape)
print (xx[0:5])
(65,)
(65,)
(65, 65)
[[-3.20000000e+00 -3.10000000e+00 -3.00000000e+00 -2.90000000e+00
  -2.80000000e+00 -2.70000000e+00 -2.60000000e+00 -2.50000000e+00
  -2.40000000e+00 -2.30000000e+00 -2.20000000e+00 -2.10000000e+00
  -2.00000000e+00 -1.90000000e+00 -1.80000000e+00 -1.70000000e+00
  -1.60000000e+00 -1.50000000e+00 -1.40000000e+00 -1.30000000e+00
  -1.20000000e+00 -1.10000000e+00 -1.00000000e+00 -9.00000000e-01
  -8.00000000e-01 -7.00000000e-01 -6.00000000e-01 -5.00000000e-01
  -4.00000000e-01 -3.00000000e-01 -2.00000000e-01 -1.00000000e-01
   2.66453526e-15  1.00000000e-01  2.00000000e-01  3.00000000e-01
   4.00000000e-01  5.00000000e-01  6.00000000e-01  7.00000000e-01
   8.00000000e-01  9.00000000e-01  1.00000000e+00  1.10000000e+00
   1.20000000e+00  1.30000000e+00  1.40000000e+00  1.50000000e+00
   1.60000000e+00  1.70000000e+00  1.80000000e+00  1.90000000e+00
   2.00000000e+00  2.10000000e+00  2.20000000e+00  2.30000000e+00
   2.40000000e+00  2.50000000e+00  2.60000000e+00  2.70000000e+00
   2.80000000e+00  2.90000000e+00  3.00000000e+00  3.10000000e+00
   3.20000000e+00]
 [-3.20000000e+00 -3.10000000e+00 -3.00000000e+00 -2.90000000e+00
  -2.80000000e+00 -2.70000000e+00 -2.60000000e+00 -2.50000000e+00
  -2.40000000e+00 -2.30000000e+00 -2.20000000e+00 -2.10000000e+00
  -2.00000000e+00 -1.90000000e+00 -1.80000000e+00 -1.70000000e+00
  -1.60000000e+00 -1.50000000e+00 -1.40000000e+00 -1.30000000e+00
  -1.20000000e+00 -1.10000000e+00 -1.00000000e+00 -9.00000000e-01
  -8.00000000e-01 -7.00000000e-01 -6.00000000e-01 -5.00000000e-01
  -4.00000000e-01 -3.00000000e-01 -2.00000000e-01 -1.00000000e-01
   2.66453526e-15  1.00000000e-01  2.00000000e-01  3.00000000e-01
   4.00000000e-01  5.00000000e-01  6.00000000e-01  7.00000000e-01
   8.00000000e-01  9.00000000e-01  1.00000000e+00  1.10000000e+00
   1.20000000e+00  1.30000000e+00  1.40000000e+00  1.50000000e+00
   1.60000000e+00  1.70000000e+00  1.80000000e+00  1.90000000e+00
   2.00000000e+00  2.10000000e+00  2.20000000e+00  2.30000000e+00
   2.40000000e+00  2.50000000e+00  2.60000000e+00  2.70000000e+00
   2.80000000e+00  2.90000000e+00  3.00000000e+00  3.10000000e+00
   3.20000000e+00]
 [-3.20000000e+00 -3.10000000e+00 -3.00000000e+00 -2.90000000e+00
  -2.80000000e+00 -2.70000000e+00 -2.60000000e+00 -2.50000000e+00
  -2.40000000e+00 -2.30000000e+00 -2.20000000e+00 -2.10000000e+00
  -2.00000000e+00 -1.90000000e+00 -1.80000000e+00 -1.70000000e+00
  -1.60000000e+00 -1.50000000e+00 -1.40000000e+00 -1.30000000e+00
  -1.20000000e+00 -1.10000000e+00 -1.00000000e+00 -9.00000000e-01
  -8.00000000e-01 -7.00000000e-01 -6.00000000e-01 -5.00000000e-01
  -4.00000000e-01 -3.00000000e-01 -2.00000000e-01 -1.00000000e-01
   2.66453526e-15  1.00000000e-01  2.00000000e-01  3.00000000e-01
   4.00000000e-01  5.00000000e-01  6.00000000e-01  7.00000000e-01
   8.00000000e-01  9.00000000e-01  1.00000000e+00  1.10000000e+00
   1.20000000e+00  1.30000000e+00  1.40000000e+00  1.50000000e+00
   1.60000000e+00  1.70000000e+00  1.80000000e+00  1.90000000e+00
   2.00000000e+00  2.10000000e+00  2.20000000e+00  2.30000000e+00
   2.40000000e+00  2.50000000e+00  2.60000000e+00  2.70000000e+00
   2.80000000e+00  2.90000000e+00  3.00000000e+00  3.10000000e+00
   3.20000000e+00]
 [-3.20000000e+00 -3.10000000e+00 -3.00000000e+00 -2.90000000e+00
  -2.80000000e+00 -2.70000000e+00 -2.60000000e+00 -2.50000000e+00
  -2.40000000e+00 -2.30000000e+00 -2.20000000e+00 -2.10000000e+00
  -2.00000000e+00 -1.90000000e+00 -1.80000000e+00 -1.70000000e+00
  -1.60000000e+00 -1.50000000e+00 -1.40000000e+00 -1.30000000e+00
  -1.20000000e+00 -1.10000000e+00 -1.00000000e+00 -9.00000000e-01
  -8.00000000e-01 -7.00000000e-01 -6.00000000e-01 -5.00000000e-01
  -4.00000000e-01 -3.00000000e-01 -2.00000000e-01 -1.00000000e-01
   2.66453526e-15  1.00000000e-01  2.00000000e-01  3.00000000e-01
   4.00000000e-01  5.00000000e-01  6.00000000e-01  7.00000000e-01
   8.00000000e-01  9.00000000e-01  1.00000000e+00  1.10000000e+00
   1.20000000e+00  1.30000000e+00  1.40000000e+00  1.50000000e+00
   1.60000000e+00  1.70000000e+00  1.80000000e+00  1.90000000e+00
   2.00000000e+00  2.10000000e+00  2.20000000e+00  2.30000000e+00
   2.40000000e+00  2.50000000e+00  2.60000000e+00  2.70000000e+00
   2.80000000e+00  2.90000000e+00  3.00000000e+00  3.10000000e+00
   3.20000000e+00]
 [-3.20000000e+00 -3.10000000e+00 -3.00000000e+00 -2.90000000e+00
  -2.80000000e+00 -2.70000000e+00 -2.60000000e+00 -2.50000000e+00
  -2.40000000e+00 -2.30000000e+00 -2.20000000e+00 -2.10000000e+00
  -2.00000000e+00 -1.90000000e+00 -1.80000000e+00 -1.70000000e+00
  -1.60000000e+00 -1.50000000e+00 -1.40000000e+00 -1.30000000e+00
  -1.20000000e+00 -1.10000000e+00 -1.00000000e+00 -9.00000000e-01
  -8.00000000e-01 -7.00000000e-01 -6.00000000e-01 -5.00000000e-01
  -4.00000000e-01 -3.00000000e-01 -2.00000000e-01 -1.00000000e-01
   2.66453526e-15  1.00000000e-01  2.00000000e-01  3.00000000e-01
   4.00000000e-01  5.00000000e-01  6.00000000e-01  7.00000000e-01
   8.00000000e-01  9.00000000e-01  1.00000000e+00  1.10000000e+00
   1.20000000e+00  1.30000000e+00  1.40000000e+00  1.50000000e+00
   1.60000000e+00  1.70000000e+00  1.80000000e+00  1.90000000e+00
   2.00000000e+00  2.10000000e+00  2.20000000e+00  2.30000000e+00
   2.40000000e+00  2.50000000e+00  2.60000000e+00  2.70000000e+00
   2.80000000e+00  2.90000000e+00  3.00000000e+00  3.10000000e+00
   3.20000000e+00]]

Now let's put some 'data' on this grid:

In [7]:
z = np.sin(xx**2 + yy**2) / (xx**2 + yy**2)
print (z.shape)
print (z[0:5])
(65, 65)
[[ 0.04874129  0.04240357  0.01978166 -0.01062917 -0.03848542 -0.05525412
  -0.05655279 -0.04273656 -0.01799396  0.01139893  0.03895874  0.05947836
   0.06985375  0.0692455   0.05873339  0.04068985  0.01808671 -0.0061086
  -0.02936306 -0.04981395 -0.0663341  -0.07847277 -0.08632353 -0.09036374
  -0.0912986  -0.08992935 -0.08705429 -0.08340291 -0.0795987  -0.07614349
  -0.07341601 -0.0716779  -0.07108182 -0.0716779  -0.07341601 -0.07614349
  -0.0795987  -0.08340291 -0.08705429 -0.08992935 -0.0912986  -0.09036374
  -0.08632353 -0.07847277 -0.0663341  -0.04981395 -0.02936306 -0.0061086
   0.01808671  0.04068985  0.05873339  0.0692455   0.06985375  0.05947836
   0.03895874  0.01139893 -0.01799396 -0.04273656 -0.05655279 -0.05525412
  -0.03848542 -0.01062917  0.01978166  0.04240357  0.04874129]
 [ 0.04240357  0.01883608 -0.01274966 -0.04093405 -0.05646843 -0.05497775
  -0.0375519  -0.00954929  0.0215723   0.04851553  0.06584536  0.07083765
   0.06350006  0.04599635  0.02177758 -0.00530575 -0.03172333 -0.05472825
  -0.07255873 -0.08442648 -0.09036374 -0.09100059 -0.08733107 -0.0805067
  -0.0716779  -0.06188818 -0.05201685 -0.04275981 -0.03463722 -0.02801656
  -0.02314226 -0.02016469 -0.01916387 -0.02016469 -0.02314226 -0.02801656
  -0.03463722 -0.04275981 -0.05201685 -0.06188818 -0.0716779  -0.0805067
  -0.08733107 -0.09100059 -0.09036374 -0.08442648 -0.07255873 -0.05472825
  -0.03172333 -0.00530575  0.02177758  0.04599635  0.06350006  0.07083765
   0.06584536  0.04851553  0.0215723  -0.00954929 -0.0375519  -0.05497775
  -0.05646843 -0.04093405 -0.01274966  0.01883608  0.04240357]
 [ 0.01978166 -0.01274966 -0.04172151 -0.0569443  -0.05375768 -0.03374625
  -0.00330033  0.02899162  0.05502907  0.06916331  0.06908742  0.05570922
   0.03232054  0.00345881 -0.02619341 -0.05264649 -0.07308788 -0.08602738
  -0.09118315 -0.08921182 -0.08137909 -0.0692441  -0.05440211 -0.03830427
  -0.02215398 -0.00686784  0.00691588  0.01879886  0.02856932  0.03614516
   0.04152139  0.04472661  0.04579094  0.04472661  0.04152139  0.03614516
   0.02856932  0.01879886  0.00691588 -0.00686784 -0.02215398 -0.03830427
  -0.05440211 -0.0692441  -0.08137909 -0.08921182 -0.09118315 -0.08602738
  -0.07308788 -0.05264649 -0.02619341  0.00345881  0.03232054  0.05570922
   0.06908742  0.06916331  0.05502907  0.02899162 -0.00330033 -0.03374625
  -0.05375768 -0.0569443  -0.04172151 -0.01274966  0.01978166]
 [-0.01062917 -0.04093405 -0.0569443  -0.05330574 -0.03174658  0.00050721
   0.0337764   0.05910014  0.0705336   0.06612808  0.04766875  0.01957246
  -0.01254908 -0.04322709 -0.0681028  -0.08442648 -0.09112789 -0.08857611
  -0.07817039 -0.06188818 -0.04188051 -0.02016469  0.0015704   0.02205529
   0.04044927  0.05629448  0.06943935  0.07995113  0.08803094  0.09393966
   0.09793814  0.10024263  0.10099446  0.10024263  0.09793814  0.09393966
   0.08803094  0.07995113  0.06943935  0.05629448  0.04044927  0.02205529
   0.0015704  -0.02016469 -0.04188051 -0.06188818 -0.07817039 -0.08857611
  -0.09112789 -0.08442648 -0.0681028  -0.04322709 -0.01254908  0.01957246
   0.04766875  0.06612808  0.0705336   0.05910014  0.0337764   0.00050721
  -0.03174658 -0.05330574 -0.0569443  -0.04093405 -0.01062917]
 [-0.03848542 -0.05646843 -0.05375768 -0.03174658  0.00178314  0.03610834
   0.06128707  0.07089339  0.06317366  0.04068985  0.00894204 -0.0253975
  -0.05609466 -0.07847277 -0.08993125 -0.08992935 -0.0795987  -0.06117283
  -0.03739583 -0.01102077  0.01554663  0.04044927  0.06244505  0.08087369
   0.09556213  0.10670186  0.11472324  0.12018254  0.12366978  0.12573922
   0.1268606   0.12738753  0.12753855  0.12738753  0.1268606   0.12573922
   0.12366978  0.12018254  0.11472324  0.10670186  0.09556213  0.08087369
   0.06244505  0.04044927  0.01554663 -0.01102077 -0.03739583 -0.06117283
  -0.0795987  -0.08992935 -0.08993125 -0.07847277 -0.05609466 -0.0253975
   0.00894204  0.04068985  0.06317366  0.07089339  0.06128707  0.03610834
   0.00178314 -0.03174658 -0.05375768 -0.05646843 -0.03848542]]

Colormesh

Now we want to make a plot of the data - right? There are many options to do this. The simplest is to plot the x,y data points and use color to represent the value of the z. For this we have plt.pcolormesh( ).

In [8]:
help(plt.pcolormesh)
Help on function pcolormesh in module matplotlib.pyplot:

pcolormesh(*args, **kwargs)
    Create a pseudocolor plot with a non-regular rectangular grid.
    
    Call signature::
    
        pcolor([X, Y,] C, **kwargs)
    
    *X* and *Y* can be used to specify the corners of the quadrilaterals.
    
    .. note::
    
       ``pcolormesh()`` is similar to :func:`~Axes.pcolor`. It's much
       faster and preferred in most cases. For a detailed discussion on
       the differences see
       :ref:`Differences between pcolor() and pcolormesh()
       <differences-pcolor-pcolormesh>`.
    
    Parameters
    ----------
    C : array_like
        A scalar 2-D array. The values will be color-mapped.
    
    X, Y : array_like, optional
        The coordinates of the quadrilateral corners. The quadrilateral
        for ``C[i,j]`` has corners at::
    
            (X[i+1, j], Y[i+1, j])          (X[i+1, j+1], Y[i+1, j+1])
                                  +--------+
                                  | C[i,j] |
                                  +--------+
                (X[i, j], Y[i, j])          (X[i, j+1], Y[i, j+1]),
    
        Note that the column index corresponds to the
        x-coordinate, and the row index corresponds to y. For
        details, see the :ref:`Notes <axes-pcolormesh-grid-orientation>`
        section below.
    
        The dimensions of *X* and *Y* should be one greater than those of
        *C*. Alternatively, *X*, *Y* and *C* may have equal dimensions, in
        which case the last row and column of *C* will be ignored.
    
        If *X* and/or *Y* are 1-D arrays or column vectors they will be
        expanded as needed into the appropriate 2-D arrays, making a
        rectangular grid.
    
    cmap : str or `~matplotlib.colors.Colormap`, optional
        A Colormap instance or registered colormap name. The colormap
        maps the *C* values to colors. Defaults to :rc:`image.cmap`.
    
    norm : `~matplotlib.colors.Normalize`, optional
        The Normalize instance scales the data values to the canonical
        colormap range [0, 1] for mapping to colors. By default, the data
        range is mapped to the colorbar range using linear scaling.
    
    vmin, vmax : scalar, optional, default: None
        The colorbar range. If *None*, suitable min/max values are
        automatically chosen by the `~.Normalize` instance (defaults to
        the respective min/max values of *C* in case of the default linear
        scaling).
    
    edgecolors : {'none', None, 'face', color, color sequence}, optional
        The color of the edges. Defaults to 'none'. Possible values:
    
        - 'none' or '': No edge.
        - *None*: :rc:`patch.edgecolor` will be used. Note that currently
          :rc:`patch.force_edgecolor` has to be True for this to work.
        - 'face': Use the adjacent face color.
        - An mpl color or sequence of colors will set the edge color.
    
        The singular form *edgecolor* works as an alias.
    
    alpha : scalar, optional, default: None
        The alpha blending value, between 0 (transparent) and 1 (opaque).
    
    shading : {'flat', 'gouraud'}, optional
        The fill style, Possible values:
    
        - 'flat': A solid color is used for each quad. The color of the
          quad (i, j), (i+1, j), (i, j+1), (i+1, j+1) is given by
          ``C[i,j]``.
        - 'gouraud': Each quad will be Gouraud shaded: The color of the
          corners (i', j') are given by ``C[i',j']``. The color values of
          the area in between is interpolated from the corner values.
          When Gouraud shading is used, *edgecolors* is ignored.
    
    snap : bool, optional, default: False
        Whether to snap the mesh to pixel boundaries.
    
    Returns
    -------
    mesh : `matplotlib.collections.QuadMesh`
    
    Other Parameters
    ----------------
    **kwargs
        Additionally, the following arguments are allowed. They are passed
        along to the `~matplotlib.collections.QuadMesh` constructor:
    
      agg_filter: a filter function, which takes a (m, n, 3) float array and a dpi value, and returns a (m, n, 3) array 
      alpha: float or None 
      animated: bool 
      antialiased or antialiaseds: Boolean or sequence of booleans 
      array: ndarray
      capstyle: unknown
      clim: a length 2 sequence of floats; may be overridden in methods that have ``vmin`` and ``vmax`` kwargs. 
      clip_box: a `.Bbox` instance 
      clip_on: bool 
      clip_path: [(`~matplotlib.path.Path`, `.Transform`) | `.Patch` | None] 
      cmap: a colormap or registered colormap name 
      color: matplotlib color arg or sequence of rgba tuples
      contains: a callable function 
      edgecolor or edgecolors: matplotlib color spec or sequence of specs 
      facecolor or facecolors: matplotlib color spec or sequence of specs 
      figure: a `.Figure` instance 
      gid: an id string 
      hatch: [ '/' | '\\' | '|' | '-' | '+' | 'x' | 'o' | 'O' | '.' | '*' ] 
      joinstyle: unknown
      label: object 
      linestyle or dashes or linestyles: ['solid' | 'dashed', 'dashdot', 'dotted' | (offset, on-off-dash-seq) | ``'-'`` | ``'--'`` | ``'-.'`` | ``':'`` | ``'None'`` | ``' '`` | ``''``]
      linewidth or linewidths or lw: float or sequence of floats 
      norm: `.Normalize`
      offset_position: [ 'screen' | 'data' ] 
      offsets: float or sequence of floats 
      path_effects: `.AbstractPathEffect` 
      picker: [None | bool | float | callable] 
      pickradius: float distance in points
      rasterized: bool or None 
      sketch_params: (scale: float, length: float, randomness: float) 
      snap: bool or None 
      transform: `.Transform` 
      url: a url string 
      urls: List[str] or None 
      visible: bool 
      zorder: float 
    
    
    See Also
    --------
    pcolor : An alternative implementation with slightly different
        features. For a detailed discussion on the differences see
        :ref:`Differences between pcolor() and pcolormesh()
        <differences-pcolor-pcolormesh>`.
    imshow : If *X* and *Y* are each equidistant, `~.Axes.imshow` can be a
        faster alternative.
    
    Notes
    -----
    
    **Masked arrays**
    
    *C* may be a masked array. If ``C[i, j]`` is masked, the corresponding
    quadrilateral will be transparent. Masking of *X* and *Y* is not
    supported. Use `~.Axes.pcolor` if you need this functionality.
    
    .. _axes-pcolormesh-grid-orientation:
    
    **Grid orientation**
    
    The grid orientation follows the standard matrix convention: An array
    *C* with shape (nrows, ncolumns) is plotted with the column number as
    *X* and the row number as *Y*.
    
    .. _differences-pcolor-pcolormesh:
    
    **Differences between pcolor() and pcolormesh()**
    
    Both methods are used to create a pseudocolor plot of a 2-D array
    using quadrilaterals.
    
    The main difference lies in the created object and internal data
    handling:
    While `~.Axes.pcolor` returns a `.PolyCollection`, `~.Axes.pcolormesh`
    returns a `.QuadMesh`. The latter is more specialized for the given
    purpose and thus is faster. It should almost always be preferred.
    
    There is also a slight difference in the handling of masked arrays.
    Both `~.Axes.pcolor` and `~.Axes.pcolormesh` support masked arrays
    for *C*. However, only `~.Axes.pcolor` supports masked arrays for *X*
    and *Y*. The reason lies in the internal handling of the masked values.
    `~.Axes.pcolor` leaves out the respective polygons from the
    PolyCollection. `~.Axes.pcolormesh` sets the facecolor of the masked
    elements to transparent. You can see the difference when using
    edgecolors. While all edges are drawn irrespective of masking in a
    QuadMesh, the edge between two adjacent masked quadrilaterals in
    `~.Axes.pcolor` is not drawn as the corresponding polygons do not
    exist in the PolyCollection.
    
    Another difference is the support of Gouraud shading in
    `~.Axes.pcolormesh`, which is not available with `~.Axes.pcolor`.
    
    .. note::
        In addition to the above described arguments, this function can take a
        **data** keyword argument. If such a **data** argument is given, the
        following arguments are replaced by **data[<arg>]**:
    
        * All positional and all keyword arguments.

In [9]:
plt.pcolormesh(x,y,z,cmap=cm.jet)
plt.axis('equal');

We probably want a color bar too:

In [10]:
plt.pcolormesh(x,y,z,cmap=cm.jet)
plt.axis('equal')
plt.colorbar();

And the plot looks a little 'pixilated' (rough). We can smooth it with the keyword shading.

In [11]:
plt.pcolormesh(x,y,z,cmap=cm.jet,shading='gouraud')
plt.axis('equal')
plt.colorbar();

Very nice. A related technique to plotting on a mesh is contouring. Contouring breaks the data down into regions separated by equal values (contours).

Contouring

One way to visualize data in a 2D grid is to make a contour plot of it. For this we use the matplotlib.pyplot function contour or contourf.

In [12]:
help(plt.contourf)
Help on function contourf in module matplotlib.pyplot:

contourf(*args, **kwargs)
    Plot contours.
    
    Call signature::
    
        contour([X, Y,] Z, [levels], **kwargs)
    
    :func:`~matplotlib.pyplot.contour` and
    :func:`~matplotlib.pyplot.contourf` draw contour lines and
    filled contours, respectively.  Except as noted, function
    signatures and return values are the same for both versions.
    
    
    Parameters
    ----------
    X, Y : array-like, optional
        The coordinates of the values in *Z*.
    
        *X* and *Y* must both be 2-D with the same shape as *Z* (e.g.
        created via :func:`numpy.meshgrid`), or they must both be 1-D such
        that ``len(X) == M`` is the number of columns in *Z* and
        ``len(Y) == N`` is the number of rows in *Z*.
    
        If not given, they are assumed to be integer indices, i.e.
        ``X = range(M)``, ``Y = range(N)``.
    
    Z : array-like(N, M)
        The height values over which the contour is drawn.
    
    levels : int or array-like, optional
        Determines the number and positions of the contour lines / regions.
    
        If an int *n*, use *n* data intervals; i.e. draw *n+1* contour
        lines. The level heights are automatically chosen.
    
        If array-like, draw contour lines at the specified levels.
        The values must be in increasing order.
    
    Returns
    -------
    :class:`~matplotlib.contour.QuadContourSet`
    
    Other Parameters
    ----------------
    corner_mask : bool, optional
        Enable/disable corner masking, which only has an effect if *Z* is
        a masked array.  If ``False``, any quad touching a masked point is
        masked out.  If ``True``, only the triangular corners of quads
        nearest those points are always masked out, other triangular
        corners comprising three unmasked points are contoured as usual.
    
        Defaults to ``rcParams['contour.corner_mask']``, which defaults to
        ``True``.
    
    colors : color string or sequence of colors, optional
        The colors of the levels, i.e. the lines for `.contour` and the
        areas for `.contourf`.
    
        The sequence is cycled for the levels in ascending order. If the
        sequence is shorter than the number of levels, it's repeated.
    
        As a shortcut, single color strings may be used in place of
        one-element lists, i.e. ``'red'`` instead of ``['red']`` to color
        all levels with the same color. This shortcut does only work for
        color strings, not for other ways of specifying colors.
    
        By default (value *None*), the colormap specified by *cmap*
        will be used.
    
    alpha : float, optional
        The alpha blending value, between 0 (transparent) and 1 (opaque).
    
    cmap : str or `.Colormap`, optional
        A `.Colormap` instance or registered colormap name. The colormap
        maps the level values to colors.
        Defaults to :rc:`image.cmap`.
    
        If given, *colors* take precedence over *cmap*.
    
    norm : `~matplotlib.colors.Normalize`, optional
        If a colormap is used, the `.Normalize` instance scales the level
        values to the canonical colormap range [0, 1] for mapping to
        colors. If not given, the default linear scaling is used.
    
    vmin, vmax : float, optional
        If not *None*, either or both of these values will be supplied to
        the `.Normalize` instance, overriding the default color scaling
        based on *levels*.
    
    origin : {*None*, 'upper', 'lower', 'image'}, optional
        Determines the orientation and exact position of *Z* by specifying
        the position of ``Z[0, 0]``.  This is only relevant, if *X*, *Y*
        are not given.
    
        - *None*: ``Z[0, 0]`` is at X=0, Y=0 in the lower left corner.
        - 'lower': ``Z[0, 0]`` is at X=0.5, Y=0.5 in the lower left corner.
        - 'upper': ``Z[0, 0]`` is at X=N+0.5, Y=0.5 in the upper left
          corner.
        - 'image': Use the value from :rc:`image.origin`. Note: The value
          *None* in the rcParam is currently handled as 'lower'.
    
    extent : (x0, x1, y0, y1), optional
        If *origin* is not *None*, then *extent* is interpreted as
        in :func:`matplotlib.pyplot.imshow`: it gives the outer
        pixel boundaries. In this case, the position of Z[0,0]
        is the center of the pixel, not a corner. If *origin* is
        *None*, then (*x0*, *y0*) is the position of Z[0,0], and
        (*x1*, *y1*) is the position of Z[-1,-1].
    
        This keyword is not active if *X* and *Y* are specified in
        the call to contour.
    
    locator : ticker.Locator subclass, optional
        The locator is used to determine the contour levels if they
        are not given explicitly via *levels*.
        Defaults to `~.ticker.MaxNLocator`.
    
    extend : {'neither', 'both', 'min', 'max'}, optional
        Unless this is 'neither', contour levels are automatically
        added to one or both ends of the range so that all data
        are included. These added ranges are then mapped to the
        special colormap values which default to the ends of the
        colormap range, but can be set via
        :meth:`matplotlib.colors.Colormap.set_under` and
        :meth:`matplotlib.colors.Colormap.set_over` methods.
    
    xunits, yunits : registered units, optional
        Override axis units by specifying an instance of a
        :class:`matplotlib.units.ConversionInterface`.
    
    antialiased : bool, optinal
        Enable antialiasing, overriding the defaults.  For
        filled contours, the default is *True*.  For line contours,
        it is taken from :rc:`lines.antialiased`.
    
    Nchunk : int >= 0, optional
        If 0, no subdivision of the domain.  Specify a positive integer to
        divide the domain into subdomains of *nchunk* by *nchunk* quads.
        Chunking reduces the maximum length of polygons generated by the
        contouring algorithm which reduces the rendering workload passed
        on to the backend and also requires slightly less RAM.  It can
        however introduce rendering artifacts at chunk boundaries depending
        on the backend, the *antialiased* flag and value of *alpha*.
    
    linewidths : float or sequence of float, optional
        *Only applies to* `.contour`.
    
        The line width of the contour lines.
    
        If a number, all levels will be plotted with this linewidth.
    
        If a sequence, the levels in ascending order will be plotted with
        the linewidths in the order specified.
    
        Defaults to :rc:`lines.linewidth`.
    
    linestyles : {*None*, 'solid', 'dashed', 'dashdot', 'dotted'}, optional
        *Only applies to* `.contour`.
    
        If *linestyles* is *None*, the default is 'solid' unless the lines
        are monochrome.  In that case, negative contours will take their
        linestyle from :rc:`contour.negative_linestyle` setting.
    
        *linestyles* can also be an iterable of the above strings
        specifying a set of linestyles to be used. If this
        iterable is shorter than the number of contour levels
        it will be repeated as necessary.
    
    hatches : List[str], optional
        *Only applies to* `.contourf`.
    
        A list of cross hatch patterns to use on the filled areas.
        If None, no hatching will be added to the contour.
        Hatching is supported in the PostScript, PDF, SVG and Agg
        backends only.
    
    
    Notes
    -----
    1. :func:`~matplotlib.pyplot.contourf` differs from the MATLAB
       version in that it does not draw the polygon edges.
       To draw edges, add line contours with
       calls to :func:`~matplotlib.pyplot.contour`.
    
    2. contourf fills intervals that are closed at the top; that
       is, for boundaries *z1* and *z2*, the filled region is::
    
          z1 < Z <= z2
    
       There is one exception: if the lowest boundary coincides with
       the minimum value of the *Z* array, then that minimum value
       will be included in the lowest interval.

In [13]:
h = plt.contour(x,y,z,cmap=cm.jet) # plot the contours
plt.axis('equal'); 
# this makes the axes "square" so we get a circle and not a squashed ellipse

Contour plots are all well and good, but what we were after was a continuous gradation of color, not the contour lines, so for that we use the function contourf instead:

In [14]:
fig = plt.contourf(x,y,z,cmap=cm.jet)
plt.axis('equal') # this makes the axes square
plt.axis('off') # and this turns off the frame, ticks and labels
plt.colorbar();

Note the difference between the colormesh (continuous gradation) and contouf (constant color between contours).

So, where were we? Oh yes, we wanted to plot the contoured elevation data onto an orthographic projection. First, we will generate a meshgrid using np.meshgrid from numpy and the method contourf of our map object ax defined above. The method ax.contourf works like matplotlib's:

In [15]:
plt.figure(1,(6,6))
ax = plt.axes(projection=ccrs.Orthographic(-75, 42))
gl=ax.gridlines(crs=ccrs.PlateCarree(),color='black',linewidth=1,linestyle='dotted')
#gl.xlabels_top = False
gl.ylocator=mticker.FixedLocator(np.arange(-80,81,20))
gl.xlocator=mticker.FixedLocator(np.arange(-180,181,60));
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
xx, yy = np.meshgrid(elons,elats)
levels=np.arange(-10000,8000,500) # define contour intervals
m=ax.contourf(xx, yy, etopo,levels,\
                transform=ccrs.PlateCarree(),
                cmap=cm.jet)
ax.set_global()
cbar=plt.colorbar(m) # put on a color bar of Elevations

Tada!!!

Application to the Earth's magnetic field

Now we know how to grid and contour data, we can try to plot other features of our planet, like the magnetic field strength or inclination (dip below or above the horizontal). I prepared a module named mkigrf for your mapping pleasure. It includes a few functions that evaluate the elements of the magnetic field for any date between 1900 and 2020 including total field strength, direction, and radial field strength. Up until 2015, the data come from the International Geomagnetic Reference Field (IGRF) model and after that, the data are extrapolated. To learn more, check out this website: https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html.

You can use the function mkigrf.doigrf to find out the magnetic declination at a particular time and place. You could then use this function to set your compass. Handy for us geology types...

But first, I should explain what the magnetic vector is. As with all vectors, it has both direction and length. We can express the vector in terms of cartesian coordinates (say, North, East and Down) or these polar coordinates:

  • declination: the angle of the magnetic direction in the horizontal plane with respect to the North pole
  • inclination: the angle of the magnetic direction in the vertical plane with respect to the horizontal
  • intensity: the strength of the field, usually in units of tesla (either nano or micro). Tesla is is magnetic induction and is usually represented by the letter B.
In [16]:
Image(filename='Figures/components.png')
Out[16]: