import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image
import pandas as pd
from matplotlib import cm
%matplotlib inline
Contour plots are really just a way to visualize something that is inherently 3D on a 2D surface. Think about our topographic map - the contour intervals are elevations and our brains can reconstruct the 3D world by looking at the contours on the map. But with computers we can visualize the 3D world in a more realistic manner. There are many 3D plotting packages that apply many different approaches to plot in 3D. For example, mplot3d, is a 3D toolkit of matplotlib that uses the same logic as for "regular" matplotlib. For more on this module, see:
http://matplotlib.sourceforge.net/mpl_toolkits/mplot3d/index.html
But for more 3D horsepower, there is a module called mlab, which is part of the mayavi package. See:
http://github.enthought.com/mayavi/mayavi/mlab.html
And then there is Mayavi itself, which comes with the Canopy Python Edition. This was way beyond what I know, but if you are curious, check out this website:
http://github.enthought.com/mayavi/mayavi/examples.html
Until now, we have used %matplotlib inline to plot our matplotlib figures within the jupyter notebook. In order to plot a 3D figure within the jupyter notebook and allow interaction with a 3D plot, we must use this alternative: %matplotlib notebook
Remember this plot from Lecture 18?
Image(filename='Figures/earthquakes_depth.png')
We used color to represent the depth of each earthquake . You can see that there are increasingly deep earthquakes as you move away from subduction zones- sort of.
It would be more instructive to view the depth of the earthquakes as points in a 3D plot. So, let's read-in the data, filter for a lat/lon box from 35$^{\circ}$S to 20$^{\circ}$N and 170$^{\circ}$-190$^{\circ}$E (the Marianas trench), and plot the data in 3D.
To begin, let's filter the data by the desired bounds and make an array of longitudes,x, an array of latitudes,y,and an array of depths,z.
# read in the data from the Lecture 18 as a Pandas DataFrame:
eq_data=pd.read_csv('Datasets/EarthquakeLocations/last5Years.csv',skiprows=1)
# define some boundaries for our box
lat_min,lon_min,lat_max,lon_max=-35,175,-15,190
# use Pandas filtering to fish out lat/lon in this range
box=eq_data[(eq_data.longitude.values%360<lon_max)&(eq_data.longitude.values%360>=lon_min)]
box=box[(box.latitude.values<lat_max)&(box.latitude.values>=lat_min)]
# and export them to NumPy arrays
x=box.longitude.values%360
y=box.latitude.values
z=-box.depth.values
We can use the Axes3D class from the mplot3d, which is another tool in the matplotlib toolkit (mpl_toolkits). [Remember that is where we got Basemap.]
We will import the Axes3D module.
But before we can use it, we have to make the notebook ready for 3D plots. One somewhat annoying thing about using the 3D options in a Jupyter notebook is that there are two different magic commands that are used. You are already familiar with the %matplotlib inline command that allows us to plot "regular" matplotlib plots in the notebook. But now we need a NEW one: %matplotlib notebook strictly for plotting 3D objects. The trouble comes when you try to do both in one notebook, because you have to always use the magic correct command before you can see your figure. So we import the Axes3D module and call the magic command:
%matplotlib notebook
# import the module
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook
The trick to using Axes3D is to first make a figure object using plt.figure( ) (called fig below) and then use the figure method fig.add_subplot( ) to make an Axis instance (here called ax). Finally we set the keyword projection of the add_sublot call to '3d'.
fig=plt.figure(1,(7,7)) # we need to make a figure object
ax = fig.add_subplot(111, projection='3d'); # so we can do this.
This makes an empty set of axes. But you can twirl them around! [Okay, I thought that was cool.]
Now we can decorate the plot with some data.
fig=plt.figure(2,(7,7)) # let's make a new one.
ax = fig.add_subplot(111, projection='3d') # so we can do this.
ax.scatter(x,y,z,c='r',marker='o')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_zlabel('Depth (km)');
Try twirling the figure around! Can you see the slab?
You can imagine other enhancements, for example setting the size of the points to be proportional to the size of the earthquake.
In Lecture 19, we learned how to make 2D color contour plots, but 3D is much more fun, so let's try to make a 3D version of something geophysical, namely the gravity anomaly of a buried sphere.
Let's choose a sphere with a radius $R$ of 2 m (whose volume is ${{4\pi}\over{3}}R^3$), that is buried $z=$ 3 m deep. The sphere's density ($\Delta \rho$) is 500 kg/m$^3$, which is, for this purpose, much higher than the surrounding material. Oh and we'll need the universal gravitational constant $G$, which is 6.67x 10$^{-11}$ Nm$^2$/kg$^2$.
The formula for gravitational attraction of such a body is:
$$g= {{4\pi}\over{3}} R^3 {{G \Delta \rho}\over{(h^2+z^2})},$$where $h$ is the horizontal distance from the mass.
The units (from dimensional analysis remembering that Newtons are kg $\cdot$ m $\cdot$ s$^2$) are m $\cdot$ s$^{-2}$. 1 Gal (for Galileo) = 1 cm $\cdot$s$^{-2}$, so to convert $g$ to the units of microgals, we multiply the above by 10$^{8}$.
We can write our equation as a lambda function.
gravity= lambda G,R,drho,h,z : (1e8*G*4.*np.pi*R**3.*drho)/(3.*(h**2+z**2)) # gravitational attraction.
And set up the arrays the same way as we did in the previous lectures as a color contour on a 2D surface.
from matplotlib import tri
x=np.arange(-6.,6.,0.1) # range of x values
y=np.arange(-6,6.,0.1) # range of y values
X,Y=np.meshgrid(x,y) # make a meshgrid
To get the gravity array $g$, we need $z,G,R$ and $\Delta \rho$ and $h$ the horizontal distance from ground zero of the buried sphere, which is given by good old Pythagorous as:
$$ h=\sqrt {x^2 + y^2}. $$# define the other parameters
z=3.
G=6.67e-11 # grav constant in Nm^2/kg^2 (SI)
R=2. # radius in meters
drho=500 # density contrast in kg/m^3
h=np.sqrt(X**2+Y**2) # get the horizontal distance from ground zero for x,y
# and make the g array
g=gravity(G,R,drho,h,z) # multipy by a million to get the units reasonable for the plots.
We want to make the plot of the gravitational attraction first using our old friend from the Lecture 19, plt.pcolormesh( ).
# put this back to 'normal' for use in the notebook using the correct Jupyter magic command
%matplotlib inline
plt.figure(3)
plt.pcolormesh(x,y,g,cmap=cm.jet,shading='gouraud')
plt.axis('equal')
plt.axis('off');
Well that IS pretty. But now we want to do this in 3D.
There are a few ways to plot this in 3D: wireframes and surfaces. Let's start with a wireframe plot. For this we need to use matplotlib Axes3D.
There is an Axes3D method called plot_wireframe( ) which will do the trick. It is in many ways similar to the pcolormesh( ), but instead of using color to represent values, it uses a 3D projection.
Let's try the plot_wireframe( ) function on our gravity data.
In the following plot, we create an Axes3D instance called ax from the figure object, fig. By setting $projection='3d'$, we get a 3D projection.
Remember that because we were just plotting in 2D (%matplotlib inline), we have to call the 3D version (%matplotlib notebook).
# makes a 3D plot in a notebook
%matplotlib notebook
fig=plt.figure(4,(8,5)) # make a figure object
ax=fig.add_subplot(111,projection='3d') # give it the powers of an Axes3D object
surf=ax.plot_wireframe(X,Y,g) # use plot_wireframe to do the plot.
ax.set_xlabel('X') # and label the axes
ax.set_ylabel('Y')
ax.set_zlabel('Z');
Try twirling it around - pretty cool, huh?
There is a related method plot_surface( ) which is also nice. We then call plot_surface( ) on the Axes3D instance and set labels to the 3 axis. Note that if we assign the surface instance to a variable, e.g., surf, we can do other things to it, like add a color bar. Finally, please admire the use of a colormap (cmap=cm.rainbow).
This works pretty much like plot_wireframe( ) but this one we can put a color bar on too.
# make the notebook ready for interactive 3D plots again
%matplotlib notebook
fig=plt.figure(4,(8,5)) # make a figure object
ax=fig.add_subplot(111,projection='3d') # get a 3d plot object
surf=ax.plot_surface(X,Y,g,cmap=cm.jet) # use plot_surface to do the plot using a color map
ax.set_xlabel('X') # and label the axes
ax.set_ylabel('Y')
ax.set_zlabel('Z')
bar=fig.colorbar(surf) # put on the color bar
bar.set_label('Gravity in $\mu$gal'); #label the color bar
OR - we can plot the data as a 3D color contour plot. For this, we can use ax.contour( ).
%matplotlib notebook
fig=plt.figure(5,(8,5)) # make a figure object
ax=fig.add_subplot(111,projection='3d') # give it the powers of an Axes3D object
surf=ax.contour(X,Y,g,cmap=cm.jet) # use contour to do the plot using a color map
ax.set_xlabel('X') # and label the axes
ax.set_ylabel('Y')
ax.set_zlabel('Z');
And ax.contour3D() does the same thing. Here, I've added a color bar.
# make the notebook ready for interactive 3D plots again
%matplotlib notebook
fig=plt.figure(4,(8,5)) # make a figure object
ax=fig.add_subplot(111,projection='3d') # get a 3d plot object
surf=ax.contour3D(X,Y,g,cmap=cm.jet) # use plot_surface to do the plot using a color map
ax.set_xlabel('X') # and label the axes
ax.set_ylabel('Y')
ax.set_zlabel('Z')
bar=fig.colorbar(surf) # put on the color bar
bar.set_label('Gravity in $\mu$gal'); #label the color bar
We can even make a wedding cake setting the extend3d keyword of the ax.contour( ) function to True.
%matplotlib notebook
fig=plt.figure(5,(8,5)) # make a figure object
ax=fig.add_subplot(111,projection='3d') # give it the powers of an Axes3D object
surf=ax.contour(X,Y,g,extend3d=True,cmap=cm.jet) # use contour to do the plot using a color map
Or to really get into it, we can plot contours onto projected planes using the ax.contour( ) method. The keyword offset sets where each plane will be put and the zdir keyword says which axis to project. You just have to play with this a while, to see how it works.
%matplotlib notebook
fig=plt.figure(6,(8,5)) # make a figure object
ax=fig.add_subplot(111,projection='3d') # give it the powers of an Axes3D object
surf=ax.plot_surface(X,Y,g,alpha=0.3)
cset=ax.contour(X,Y,g,zdir='z',offset=1,cmap=cm.jet)
cset=ax.contour(X,Y,g,zdir='x',offset=-7,cmap=cm.jet)
cset=ax.contour(X,Y,g,zdir='y',offset=7,cmap=cm.jet)
ax.set_xlabel('X') # and label the axes
ax.set_ylabel('Y')
ax.set_zlabel('Z');
Before we leave the gravity anomaly problem, consider one more way to plot the data. Gravity data are inherently vectors, so we could plot them as arrows. This can be done using the plt.quiver( ) method.
help(plt.quiver)
Help on function quiver in module matplotlib.pyplot: quiver(*args, **kw) Plot a 2-D field of arrows. Call signatures:: quiver(U, V, **kw) quiver(U, V, C, **kw) quiver(X, Y, U, V, **kw) quiver(X, Y, U, V, C, **kw) *U* and *V* are the arrow data, *X* and *Y* set the location of the arrows, and *C* sets the color of the arrows. These arguments may be 1-D or 2-D arrays or sequences. If *X* and *Y* are absent, they will be generated as a uniform grid. If *U* and *V* are 2-D arrays and *X* and *Y* are 1-D, and if ``len(X)`` and ``len(Y)`` match the column and row dimensions of *U*, then *X* and *Y* will be expanded with :func:`numpy.meshgrid`. The default settings auto-scales the length of the arrows to a reasonable size. To change this behavior see the *scale* and *scale_units* kwargs. The defaults give a slightly swept-back arrow; to make the head a triangle, make *headaxislength* the same as *headlength*. To make the arrow more pointed, reduce *headwidth* or increase *headlength* and *headaxislength*. To make the head smaller relative to the shaft, scale down all the head parameters. You will probably do best to leave minshaft alone. *linewidths* and *edgecolors* can be used to customize the arrow outlines. Parameters ---------- X : 1D or 2D array, sequence, optional The x coordinates of the arrow locations Y : 1D or 2D array, sequence, optional The y coordinates of the arrow locations U : 1D or 2D array or masked array, sequence The x components of the arrow vectors V : 1D or 2D array or masked array, sequence The y components of the arrow vectors C : 1D or 2D array, sequence, optional The arrow colors units : [ 'width' | 'height' | 'dots' | 'inches' | 'x' | 'y' | 'xy' ] The arrow dimensions (except for *length*) are measured in multiples of this unit. 'width' or 'height': the width or height of the axis 'dots' or 'inches': pixels or inches, based on the figure dpi 'x', 'y', or 'xy': respectively *X*, *Y*, or :math:`\sqrt{X^2 + Y^2}` in data units The arrows scale differently depending on the units. For 'x' or 'y', the arrows get larger as one zooms in; for other units, the arrow size is independent of the zoom state. For 'width or 'height', the arrow size increases with the width and height of the axes, respectively, when the window is resized; for 'dots' or 'inches', resizing does not change the arrows. angles : [ 'uv' | 'xy' ], array, optional Method for determining the angle of the arrows. Default is 'uv'. 'uv': the arrow axis aspect ratio is 1 so that if *U*==*V* the orientation of the arrow on the plot is 45 degrees counter-clockwise from the horizontal axis (positive to the right). 'xy': arrows point from (x,y) to (x+u, y+v). Use this for plotting a gradient field, for example. Alternatively, arbitrary angles may be specified as an array of values in degrees, counter-clockwise from the horizontal axis. Note: inverting a data axis will correspondingly invert the arrows only with ``angles='xy'``. scale : None, float, optional Number of data units per arrow length unit, e.g., m/s per plot width; a smaller scale parameter makes the arrow longer. Default is *None*. If *None*, a simple autoscaling algorithm is used, based on the average vector length and the number of vectors. The arrow length unit is given by the *scale_units* parameter scale_units : [ 'width' | 'height' | 'dots' | 'inches' | 'x' | 'y' | 'xy' ], None, optional If the *scale* kwarg is *None*, the arrow length unit. Default is *None*. e.g. *scale_units* is 'inches', *scale* is 2.0, and ``(u,v) = (1,0)``, then the vector will be 0.5 inches long. If *scale_units* is 'width'/'height', then the vector will be half the width/height of the axes. If *scale_units* is 'x' then the vector will be 0.5 x-axis units. To plot vectors in the x-y plane, with u and v having the same units as x and y, use ``angles='xy', scale_units='xy', scale=1``. width : scalar, optional Shaft width in arrow units; default depends on choice of units, above, and number of vectors; a typical starting value is about 0.005 times the width of the plot. headwidth : scalar, optional Head width as multiple of shaft width, default is 3 headlength : scalar, optional Head length as multiple of shaft width, default is 5 headaxislength : scalar, optional Head length at shaft intersection, default is 4.5 minshaft : scalar, optional Length below which arrow scales, in units of head length. Do not set this to less than 1, or small arrows will look terrible! Default is 1 minlength : scalar, optional Minimum length as a multiple of shaft width; if an arrow length is less than this, plot a dot (hexagon) of this diameter instead. Default is 1. pivot : [ 'tail' | 'mid' | 'middle' | 'tip' ], optional The part of the arrow that is at the grid point; the arrow rotates about this point, hence the name *pivot*. color : [ color | color sequence ], optional This is a synonym for the :class:`~matplotlib.collections.PolyCollection` facecolor kwarg. If *C* has been set, *color* has no effect. Notes ----- Additional :class:`~matplotlib.collections.PolyCollection` keyword arguments: 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 -------- quiverkey : Add a key to a quiver plot
So we need the x,y coordinates of the gravity vector and plot the arrows at the evaluation points.
# we need to redo what we already did. but at lower resolution
%matplotlib inline
z=3.
G=6.67e-11 # grav constant in Nm^2/kg^2 (SI)
R=2. # radius in meters
drho=500 # density contrast in kg/m^3
h=np.sqrt(X**2+Y**2) # get the horizontal distance from ground zero for x,y
# and make the g array
g=gravity(G,R,drho,h,z) # multipy by a million to get the units reasonable for the plots.
x=np.arange(-6,6.5,1) # range of x values
y=np.arange(-6.,6.5,1) # range of y values
X,Y=np.meshgrid(x,y) # make a meshgrid
h=np.sqrt(X**2+Y**2) # get the horizontal distance from ground zero for x,y
g=gravity(G,R,drho,h,z) # re-use our lambda function
# now make the horizontal projections along X and Y
U=-X*g
V=-Y*g
# and plot
plt.quiver(x,y,U,V)
plt.axis('equal');
There is a lovely plotting package with lots more options (although not 3D). One of them is a nicer quiver plot option, plotly.create_quiver. Let's take a look:
import plotly.figure_factory as ff
%matplotlib inline
fig=ff.create_quiver(X,Y,U,V)
fig.show() # for some reason, this plotting option requires this line.
Let's return for a moment to the topic of Lecture 17 where we plotted the data with seaborn.pairplot which showed some interesting relationships. Now we can plot the different components on a 3D plot. We can use a new trick to create a discrete color map based on an input color map, too. Here's a function discrete_cmap that does it:
def discrete_cmap(N, base_cmap=None):
"""Create an N-bin discrete colormap from the specified input map"""
# Note that if base_cmap is a string or None, you can simply do
# return plt.cm.get_cmap(base_cmap, N)
# The following works for string, None, or a colormap instance:
base = plt.cm.get_cmap(base_cmap)
color_list = base(np.linspace(0, 1, N))
cmap_name = base.name + str(N)
return color_list
So we read in the data from Lecture 17. We can use some concepts from Lecture 24 and apply a PCA approach to split the data into the principal components.
MantleArray=pd.read_csv('Datasets/GeoRoc/MantleArray_OIB.csv')
from sklearn.decomposition import PCA
X_PCA=np.array([MantleArray['EPSILON_ND'],MantleArray['SR87_SR86'],MantleArray['PB206_PB204'],\
MantleArray['EPSILON_HF']])
X_PCA=X_PCA.T
pca = PCA(3)
pca.fit(X_PCA)
Y_pca = pca.transform(X_PCA)
MantleArray['PCA_1']=Y_pca[:,0]
MantleArray['PCA_2']=Y_pca[:,1]
MantleArray['PCA_3']=Y_pca[:,2]
Now let's plot the end members in 3D!
%matplotlib notebook
fig=plt.figure(2,(7,7)) # let's make a new one.
ax=fig.add_subplot(111,projection='3d') # so we can do this.
locations=MantleArray['LOCATION'].unique() # get a list of unique locations
colors=discrete_cmap(len(locations),base_cmap='nipy_spectral') # get the color map
# assign the colors
MantleArray['COLOR']=[[0,0,0,0] for n in list(range(len(MantleArray)))]
# step through the locations and plot the data using the specified colors
for n in list(range(len(locations))):
newslice=MantleArray.loc[MantleArray.LOCATION==locations[n],'COLOR']
colorlist=pd.Series([colors[n] for m in list(range(len(newslice)))],index=newslice.index)
MantleArray.loc[MantleArray.LOCATION==locations[n],'COLOR']=colorlist
ax.scatter(MantleArray['PCA_1'],MantleArray['PCA_2'],MantleArray['PCA_3'],c=MantleArray['COLOR'].values);
This lecture was just a sample of what is possible with 3D plotting. Here are more examples that may inspire you. Check out this website: https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html for hints on how to do more.
Image(filename='Figures/mplot3d-examples.png')