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")
Remember this figure?
Image(filename='Figures/etopo20.jpg',width=500)
We now have the tools to make it!
First, you should recognize the basic orthographic projection from last lecture.
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();
We already encountered the topographic data in the lecture on the hypsographic curve. Let's import them here and look at their shapes again:
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.
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.
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:
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]]
Now we want to make a plot of the data - right? There are many options to do this. One rather simple way would be to plot the x,y data points and use color to represent the value of the z. For this we have plt.pcolormesh( ).
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.
plt.pcolormesh(x,y,z,cmap=cm.jet)
plt.axis('equal');
We probably want a color bar too:
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.
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).
Another way to visualize data in a 2D grid would be to make a contour plot of it. For this we use the matplotlib.pyplot function contour or contourf.
help(plt.contour)
Help on function contour in module matplotlib.pyplot: contour(*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.
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:
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:
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.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
X, Y = np.meshgrid(elons,elats)
levels=np.arange(-10000,8000,500) # define contour intervals
m=ax.contourf(X, Y, etopo,levels,\
transform=ccrs.PlateCarree(),
cmap=cm.jet)
ax.set_global()
cbar=plt.colorbar(m) # put on a color bar of Elevations
Tada!!!
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:
Image(filename='Figures/components.png')
a) Lines of flux produced by a geocentric axial dipole. b) Lines of flux of the geomagnetic field of 2005. At point P the horizontal component of the field $B_H$, is directed toward the magnetic north. The vertical component $B_V$ is directed down and the field makes an angle $I$ with the horizontal, known as the inclination. c) Components of the geomagnetic field vector ${\bf B}$. The angle between the horizontal component (directed toward magnetic north and geographic north is the declination $D$.) [Modified from Ben-Yosef et al., 2008, doi:10.1016/j.jas.2008.05.016]
import mkigrf as mkigrf
help(mkigrf.doigrf)
Help on function doigrf in module mkigrf: doigrf(long, lat, date) Calculates the interpolated (<2015) or extrapolated (>2015) main field and secular variation coefficients and passes them to the Malin and Barraclough routine (function pmag.magsyn) to calculate the field from the coefficients. Parameters: ----------- lon : east longitude in degrees (0 to 360 or -180 to 180) lat : latitude in degrees (-90 to 90) date : Required date in years and decimals of a year (A.D.) Return ----------- x : north component of the magnetic field in nT y : east component of the magnetic field in nT z : downward component of the magnetic field in nT f : total magnetic field in nT By default, igrf12 coefficients are used between 1900 and 2020 from http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html. To check the results you can run the interactive program at the NGDC www.ngdc.noaa.gov/geomag-web
mkigrf.doigrf returns $x,y,z$ cartesian components of the magnetic field vector. But we want to plot the polar coordinates declination, inclination, and strength. So we need to convert from cartesian coordinates to polar coordinates. There is a handy function mkigrf.cart2dir() that will do this for you. For example, we find the declination for San Diego in 2017 like this:
help(mkigrf.cart2dir)
Help on function cart2dir in module mkigrf: cart2dir(x, y, z) Converts a direction in cartesian coordinates into declination, inclinations Parameters ---------- cart : input list of [x,y,z] or list of lists [[x1,y1,z1],[x2,y2,z2]...] Returns ------- direction_array : returns an array of [declination, inclination, intensity] Examples -------- >>> pmag.cart2dir([0,1,0]) array([ 90., 0., 1.])
San_lat=33
San_lon=243
x,y,z,f=mkigrf.doigrf(San_lon,San_lat,2018)
Dec,Inc,B=mkigrf.cart2dir(x,y,z)
print ('%7.1f'%(Dec)) # notice the formatting
11.5
We can use the tools in mkigrf to evaluate the magnetic field over the surface of the Earth and then contour it on, say, a Mollweide projection. There is another tool in mkigrf called magMap( ) that comes in handy now. It will generate the data to make a global map of the geomagnetic field which we can then plot on a projection of our choosing (say, Mollweide) using the tools in cartopy.
help(mkigrf.magMap)
Help on function magMap in module mkigrf: magMap(date, **kwargs) generates the data for a map of the magnetic field. Inputs: required: date = decimal year for evaluation (between 1900 and 2020) optional: lon_0 = desired zero longitude Returns: Bdec = declinations Binc = inclinations B = field strength (in microtesla) lons = array of longitudes lats = array of latitudes
We can call magMap for any date we like. We also have to specify the central longitude that we plan to use in the Mollweide projection. Here I use the Greenwich meridian (0).
date=2015 # let's do this for 2018
lon_0=0 # we can specify the grid spacing and the intended 0 longitude for the plot
Ds,Is,Bs,lons,lats=mkigrf.magMap(date,lon_0=lon_0)
And.... drum roll .... we can plot them on a contour map.
Here is the field strength data ($Bs$).
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=lon_0))
X, Y = np.meshgrid(lons,lats)
lincr=2
levmax=round(Bs.max())+lincr
levmin=round(Bs.min()-lincr)
levels=np.arange(levmin,levmax,lincr)
m=ax.contourf(X, Y, Bs, levels,transform=ccrs.PlateCarree(),cmap=cm.jet)
ax.set_global()
ax.coastlines()
cbar=plt.colorbar(m,orientation='horizontal') # put on a color bar of intensities
plt.title('Field strength ($\mu$T): '+str(date));
And the inclinations ($Is$):
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=lon_0))
X, Y = np.meshgrid(lons,lats)
lincr=1
levmax=Is.max()+lincr
levmin=round(Is.min()-lincr)
levels=np.arange(levmin,levmax,lincr)
ax.contourf(X,Y, Is, levels,transform=ccrs.PlateCarree(),cmap=cm.jet)
ax.contour(X,Y,Is,levels=np.arange(-80,90,10),transform=ccrs.PlateCarree(),colors='black')
ax.set_global()
ax.coastlines()
cbar=plt.colorbar(m,orientation='horizontal') # put on a color bar of intensities
plt.title('Inclination: '+str(date));
Okay - i love magnetic fields and I think these maps are cool.