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]: