<jupyter_start><jupyter_text>2 Visualising your data Scientific undertakings are rather pointless if your findings are not communicated to someone else (a colleague, the scientific community, the public, etc.)The gold standard of scientific communication is the ***figure*** - a presentation of your findings in pictorial form that maximizes reader comprehension of your ideas. In this notebook, we will consider a range of Python visualisations, from simple line plots, to complex dynamic figures. I will attempt (and fail) to be exhaustive, as I hope this notebook can serve as a sort of reference when you come to do your own figures (*"what was the command for a scatter plot again?"*)"***Why have you tried to cram so much into this notebook?***"*Think about learning Python as like learning a new language. Once you have complete mastery of the vocabulary (and the rules for using it) you have total flexibility in how you express yourself... BUT you have to **know** a word exists before you can **use** it. **Python is similar - you have total control over all aspects of plotting, but first you need to know the commands that do the different things. Hence, my attempt to be exhaustive here.*  2.1D - line plots, histograms, subplots The workhorse of Python plotting is the [matplotlib](https://matplotlib.org/index.html) module (a particularly good reference is the [gallery](https://matplotlib.org/gallery.html) - choose your plot type by **thumbnail** sample and click to get the **source code**).We can make Python plots **display** in the notebook by calling a special IPython [magic](http://ipython.readthedocs.io/en/stable/interactive/python-ipython-diff.html?highlight=magicsmagics) function.<jupyter_code>%matplotlib inline<jupyter_output><empty_output><jupyter_text>2.1.1 The basics Now, let's **import** the matplotlib module and create an **empty figure and axes**.<jupyter_code>from matplotlib import pyplot as plt      
f,ax = plt.subplots(1,1)                  # creates a 'matrix' of subplots, in this case, 1 row and 1 column

# CHANGE THIS PLOT to create two subplots, one on top of the other<jupyter_output><empty_output><jupyter_text>The `subplots` function returns two outputs: a 'figure' object, which we have assigned to the variable `f` and a set of 'axes' objects (in this case, just one), which we have assigned to the variable `ax`.By default, both axes plot from 0 to 1, although this will be automatically adjusted as we add lines, points, data, etc. We can **resize** the figure and **save a copy** of it to the file system.<jupyter_code>f,ax = plt.subplots(1,1)                           # 1 row of subplots, 1 column of subplots
f.set_size_inches(3,8)                             # 3 inches wide, 8 inches high
plt.savefig('my_empty_figure.png', dpi = 300)      # dpi is the resolution, 300 is good for many applications

# **to do**
# CHANGE THIS PLOT so that the figure is short and wide<jupyter_output><empty_output><jupyter_text>***Verify that the file `my_empty_figure.png` has been created and corresponds to the awkwardly tall axis above.******Modify the code above to produce two SQUARE subplots, side-by-side.*** You will need to modify BOTH the lines `f,ax = plt.subplots(1,1)` and `f.set_size_inches(3,8)`.  2.1.2 Line plots Let's make some simple line plots.<jupyter_code>import numpy as np 
f,ax = plt.subplots(1,1,figsize=(8,4))            # can set figsize here directly

# get some topograpy data
x,z = np.genfromtxt('../data/topo.csv',delimiter=',',skip_header=1).T

ax.plot(x,z)                                      # the 'plot' command is actually a 'method' called on the 'ax' object<jupyter_output><empty_output><jupyter_text>Want to change the ***color*** of the line?<jupyter_code>f,ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(x,z,'r')                                 # r = red, k = black, b = blue, g = green, c = cyan, m = magenta, y = ??  

# **to do**
# CHANGE THIS PLOT to plot a green line<jupyter_output><empty_output><jupyter_text>Want to change the ***style*** of the line? And add a marker?<jupyter_code>f,ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(x,z,'r--s')                             # -- = dashed, : = dotted, -. = dash-dot, - = solid
                                                # s = square, ^ = up triangle, o = circle, . = dot, * = star, p = pentagon
# **to do**
# CHANGE THIS PLOT to plot a blue dotted line with circle markers<jupyter_output><empty_output><jupyter_text>*"That plot: 3/10. It doesn't have axes labels or a title."*<jupyter_code>f,ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(x,z,'r--s')  

ax.set_xlabel('x')                                  # add a basic label to the x-axis
ax.set_ylabel('z', size = 12)                       # add a large label to the y-axis
ax.set_title(r'topography along profile $\alpha$', size = 14)      # add a title using LaTeX math mode

# **to do**
# CHANGE THIS PLOT so that:
# - the xlabel reads "along strike distance / km"
# - the ylabel reads "elevation / m"<jupyter_output><empty_output><jupyter_text>*"6/10, the plot shouldn't exceed the limits of x axis, it's unseemly."*<jupyter_code>f,ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(x,z,'r--s')      
ax.set_xlabel('x')                                  
ax.set_ylabel('z', size = 12)                       
ax.set_title(r'topography along profile $\alpha$', size = 14)       

# set plot limits in x direction to '0' on the left and 'the last element of the array x' on the right
ax.set_xlim([0, x[-1]])      

# **to do**
# CHANGE THIS PLOT so that the ylimits are set to 300 to 600<jupyter_output><empty_output><jupyter_text>*"8/10 - I'd like a legend and a text annotation for the subplot."* Of course!<jupyter_code>f,ax = plt.subplots(1,1,figsize=(8,4))
ax.plot(x,z,'r--s', label = 'profile $\\alpha$')                # add a label when you call the plot command      
ax.set_xlabel('x')                                  
ax.set_ylabel('z', size = 12)                       
ax.set_title(r'topography along profile $\alpha$', size = 14)    
ax.set_xlim([0, x[-1]])     

# add a legend
ax.legend(loc = 1)                                     # 1 = top right, 2 = top left, 0 = default, etc ...
# add a text annotation in the bottom left (this is great when you're doing subfigures)
ax.text(0.05, 0.05, "A", ha = 'left', va = 'bottom', transform=ax.transAxes, size = 15)

# **to do**
# CHANGE THIS PLOT so that:
# - the legend label reads "topography"
# - the legend appears in the bottom righthand corner
# - the text annotation is in the top lefthand corner and reads 'I'<jupyter_output><empty_output><jupyter_text>*"9/10 - Can you put more data on there? And another subplot?"*<jupyter_code># we'll add a second subplot (1 row, 2 columns) and 'unpack' the axes as 'ax1' and 'ax2'
f,(ax1,ax2) = plt.subplots(1,2,figsize=(10,4))
ax1.plot(x,z,'r--s', label = 'profile $\\alpha$')    

xv = np.linspace(x[0],x[-1],101)
ax1.plot(xv, 110*np.cos(xv/x[-1]*2*np.pi)+420, 'g--', label='not a great model')       # transform and plot the data on the first axis    

ax2.plot(x, z/1.e3, 'k-')       # plot a different transformation on the second axis

# let's set the labels for all axes simultaneously using a FOR loop
for ax in [ax1,ax2]:
    ax.set_xlabel('x')                                
    ax.set_ylabel('z', size = 12)            

# add a legend for the first axis 
ax1.legend(loc = 1)       
# add a text annotation in the two plots
ax1.text(0.05, 0.05, "A", ha = 'left', va = 'bottom', transform=ax1.transAxes, size = 15)
ax2.text(0.05, 0.05, "B", ha = 'left', va = 'bottom', transform=ax2.transAxes, size = 15)

# **to do**
# CHANGE THIS PLOT so that:
# - the 2nd axis plots the line y propto b(1-exp(-a(x-xm)**2))
# - the 2nd axis has a legend, your choice of labels for the two lines<jupyter_output><empty_output><jupyter_text>***10/10***As you can see, plotting commands are additive and, depending on the complexity of your figure, will rapidly compound into a quite sizeable Python script.(*This code also implemented in script: [module 2 - visualisation/line_plot.py](line_plot.py) - try running it at the command line* `ipython line_plot.py`)  Replicate the plot below (exercise)<jupyter_code># **to do**
# write python commands to produce the same figure as above<jupyter_output><empty_output><jupyter_text>2.1.3 Histograms These are an okay way to present the **'frequency'** of binned values within a dataset. First, let's generate some **random** 'data'.<jupyter_code>N = 100                       # number of data points to create
obs = np.random.randn(N)      # generate normally distributed data
mean, std = [2.5, 0.3]        # mean and standard deviation to transform the data
obs = obs*std + mean          # transform the data
print(obs)<jupyter_output><empty_output><jupyter_text>Before we plot a histogram, let's look at some useful functions for ***summarising*** the data.<jupyter_code># the basics
print('min =', np.min(obs))
print('max =', np.max(obs))
print('mean =', np.mean(obs))
print('std =', np.std(obs))

# more fancy
print('LQ =', np.percentile(obs,25))
print('median =', np.percentile(obs,50))
print('UQ =', np.percentile(obs,75))
print('count = ', len(obs))
print('sum =', np.sum(obs))<jupyter_output><empty_output><jupyter_text>Let's plot a **histogram**: there are two steps and some **decision-making**.1. Place the data into **bins**. This requires an **expert decision** on how many bins to use (rule-of-thumb, for $N$ data, use $\sqrt{N}/2$ bins)2. Plot the histogram using the `bar` method.<jupyter_code># create the "bin" as a vector of evenly spaced points
Nbins = 8                                                         # alternatively, Nbins = int(np.sqrt(N)/2)
bin_edges = np.linspace(np.min(obs), np.max(obs), Nbins+1)        # create the bins, a vector of bin edges
h,e = np.histogram(obs, bin_edges)                                # sort the data into their bins
print(e)                                                          # bin edges
print(h)                                                          # bin heights

# create the histogram
f,ax = plt.subplots(1,1)
ax.bar(e[:-1], height = h, width = e[1]-e[0], color = [0, 0.5, 0], edgecolor = 'k')
# Above, I have specified the bar 'color' as an [r g b] vector - this gives more flexibility in the type of colors that can
# be used. Experiment with changing the three values (make sure your values are between 0 and 1 though!)

# **to do**
# CHANGE THIS PLOT so that:
# - there are 20 bins
# - there are 1000 data points<jupyter_output><empty_output><jupyter_text>***Is the shape of the plot as you expect? Why is it not perfectly normally distributed?******What happens as you increase the number of bins up or down? Describe the trade-offs.***(*This code also implemented in script: [2_visualisation/histogram.py](histogram.py) - try running it at the command line* `ipython histogram.py`)  &lt;homework&gt; The task below requires you to selectively copy-paste from the code above and modify for your purposes.<jupyter_code># let's read in the production well data from python101.ipynb
time, mf = np.genfromtxt('../data/PW1.dat',delimiter=',',skip_header=True).T

# **to do**
# - create a series of bin_edges using the min and max of the data 'mf'
# - plot a histogram of 'mf'
# - extra: 'detrend' the data first by subtracting the mean of 'mf'<jupyter_output><empty_output><jupyter_text>&lt;/homework&gt;  2.2D - scatter plots, contour plots, colour maps, cartopy, KDE, images Much of what we covered in the previous section - axis labels, legends, line styles - applies to the plots in this section. For brevity though, I won't include these commands.  2.2.1 Scatter plots This is a favourite of mine for plotting earthquake locations, because we can modify the **size** and **colour** of the marker to reflect earthquake magnitude and time.Let's have a look .***First, ensure you have executed the code below Cell 1.5.2 in the previous notebook - this will download some earthquake data***.<jupyter_code># load the earthquake data downloaded previously
data = np.genfromtxt('../1_data/nd_eqs.txt', delimiter = ',', skip_header=1)

# we will pick out four columns: date (col 1), lat (col 4), lon (col 5), magnitude (col 7)
date_all = data[:,0]                # extract all the rows, first column
lat_all = data[:,3]/180.*np.pi      # same, fourth column (coverted from degrees to radians)
lon_all = data[:,4]/180.*np.pi      # etc
mag_all = data[:,6]

# convert lat-lon to approximate x-y
r = 6371                        # radius of earth (km)
lat_avg = np.mean(lat_all)
x_all = r*lon_all*np.cos(lat_avg)
y_all = r*lat_all<jupyter_output><empty_output><jupyter_text>Let's just do a really basic x-y scatter plot.<jupyter_code>f,ax = plt.subplots(1,1,figsize=(6,6))                # this command should be familiar by now...
ax.scatter(x_all,y_all)<jupyter_output><empty_output><jupyter_text>Looks like most of the events are clustered in a single region $410<x<470$ and $5800<y<5950$. Let's just focus on those events.To do that, we're going to need to **filter** our four data arrays - `date`, `x`, `y`, and `mag` - to retain *only events in the range we're interested in*.<jupyter_code># use the 'where' function and conditions to find the data indices we're interested in 
    # (while we're at it, get rid of the events smaller than M 0)
    # & = 'and' but when comparing arrays
inds = np.where((435 < x_all) & (x_all < 475) & (5900 < y_all) & (y_all < 5950) & (mag_all > 0.0))      


# 'slice' the arrays to keep only the values satisfying the conditions above
date = date_all[inds]
x = x_all[inds]
y = y_all[inds]
mag = mag_all[inds]

# plot the clipped data
f,ax = plt.subplots(1,1,figsize=(6,6))   
ax.scatter(x,y)
ax.set_aspect('equal', adjustable='box')    # this command avoids weird stretching of the axes which distorts the data

# **to do**
# CHANGE THIS PLOT so that only earthquakes larger than magnitude 1.0 are shown<jupyter_output><empty_output><jupyter_text>Let's now modify the **size** of the markers so that they reflect the magnitude of the earthquake.<jupyter_code>f,ax = plt.subplots(1,1,figsize=(6,6))   

# convert magnitude to a size between 1 and 100 (largest event is a 100 point circle)
    # first, rescale magnitudes to range between 0 and 1
s = (mag - np.min(mag))/(np.max(mag) - np.min(mag))
    # second, rescale to range between 1 and 100 (this will be the size of the marker)
s = s*(50-2)+1
ax.scatter(x,y,s)         # passing a third argument now, marker size
ax.set_aspect('equal', adjustable='box')    

# **to do**
# CHANGE THIS PLOT to implement a DIFFERENT scaling of the sizes<jupyter_output><empty_output><jupyter_text>Finally, let's distinguish the events based on the time they occurred. However, before we can do that, we'll need to **convert** the 'datetime' - YYYYMMDD - that we read from the text file from a float to a string, and then **reinterpret** this as a 'decimal year'. For this task, we will use the Python [datetime](https://docs.python.org/3/library/datetime.html) module.<jupyter_code># import the datetime module
from datetime import datetime

# convert the dates from float -> integer -> string
t0 = datetime.strptime('19910101', '%Y%m%d')           # reference time 1 Jan 1991
times = []                                             # an empty list for storing each date as it is calculated
for each_date in date:
    # first, convert from float -> integer -> string    (converting straight to a string leaves an awkward decimal point)
    str_date = str(int(each_date))
    # interpret each datestring 
    t = datetime.strptime(str_date, '%Y%m%d')
    # take the difference between the datestring and the reference time
    dt = t - t0
    # find the total seconds and convert the date to decimal years
    times.append(dt.total_seconds()/(3600*24*365.25)+1991)<jupyter_output><empty_output><jupyter_text>Now let's plot events with colours **corresponding** with earthquake time (and add a colour bar for good measure)<jupyter_code>f,ax = plt.subplots(1,1,figsize=(6,6))   

CS = ax.scatter(x,y,s, c = times)            # use the new 'times' array to indicate to scale a color map
ax.set_aspect('equal', adjustable='box')    
plt.colorbar(CS, ax = ax)                   # create a colorbar and attach it to the 'ax' axis<jupyter_output><empty_output><jupyter_text>Not happy with this particular colour scheme? Let's change it (here is the [full selection](https://matplotlib.org/examples/color/colormaps_reference.html)).***For super-users, consider using a [colour-blind friendly colourmap](http://colorbrewer2.org/type=sequential&scheme=BuGn&n=3).***<jupyter_code>f,ax = plt.subplots(1,1,figsize=(6,6))   
import matplotlib
coolwarm = matplotlib.cm.get_cmap('coolwarm_r')     # import a new colormap - this one is 'coolwarm', reversed it with '_r'

CS = ax.scatter(x,y,s, c= times, cmap=coolwarm)     # pass the colormap in here
ax.set_aspect('equal', adjustable='box')    
plt.colorbar(CS, ax = ax)               

# **to do**
# CHANGE THIS PLOT to use a sequential colormap (your choice)<jupyter_output><empty_output><jupyter_text>Phew - you're getting pretty expert at this!**But we can do more...**Let's say we have an **image** that outlines the reservoir in which these events are occurring (we do, this image is in the same folder as this notebook, `groningen_reservoir.png`). We can load the image into Python, plot it, and then superimpose the earthquakes over the top.<jupyter_code>f,ax = plt.subplots(1,1,figsize=(6,6))   

im = plt.imread('groningen_reservoir.png')                  # use the imread function to read an image file
                                                            # 'im' is a MxNx3 array of image data (each pixel has three 
                                                            # components, corresponding to an RGB color) 

implot = ax.imshow(im, extent=[430, 475, 5902, 5950])     # plot the image and 'stretch' it to the given x and y limits

# then, let's plot the earthquakes over the top
CS = ax.scatter(x,y,s, c= times, cmap=coolwarm)     
ax.set_aspect('equal', adjustable='box')    
plt.colorbar(CS, ax = ax)<jupyter_output><empty_output><jupyter_text>Starting to look pretty good right? You get the idea how a complex figure can be **built up** by just starting with some simple commands and then introducing **incremental improvements**.(*This code also implemented in script: [2_visualisation/scatter_plot.py](scatter_plot.py) - try running it at the command line* `ipython scatter_plot.py`)In the next section, we'll consider a **variation** on this figure, mainly as a vehicle to introduce a few concepts that are important for 2D plotting.  &lt;homework&gt; Read in a set of porosity-permeability-rock type data and use a `scatter` plot to illustrate the poroperm relationship for each rock (color by rock type).<jupyter_code># read in the poroperm data
phi,k,rt = np.genfromtxt('../data/poroperm.csv',delimiter=',',skip_header=1).T

# **to do**
# uncomment and complete the commands below to scatter plot permeability vs. porosity
# color the markers by rock type 'rt'

#f,ax =???
#ax.set_yscale('log')                       # use a log axes for permeability
#ax.scatter(???<jupyter_output><empty_output><jupyter_text>&lt;/homework&gt;  2.2.2 Contour plots Let's say that instead of plotting the **individual** earthquakes themselves, we'd rather some visualisation of earthquake **density**. One way to go about this is the **histogram route** - divide the area (this time 2D instead of 1D) into **rectangular bins**, **count** the earthquakes in each, and then **assign a colour** corresponding to the number of earthquakes.I have done this in the cell below<jupyter_code>f,ax = plt.subplots(1,1)                
f.set_size_inches(10,10)
im = plt.imread('groningen_reservoir.png')                  
implot = ax.imshow(im, extent=[430, 475, 5902, 5950])     

# use histogram2d to add up EQs in discrete bins
xbin_edges = np.linspace(430,475,int((475-430)/1)+1)                    # divide x limits into 1 km bins
ybin_edges = np.linspace(5899, 5950, int((5950-5899)/1)+1)              # same for y
dens, xe, ye = np.histogram2d(x,y, [xbin_edges,ybin_edges])        # create the counting array
dens = np.fliplr(dens)                                             # flip the matrix columns, left to right

# create an empty RGBA figure (3 color components, 1 'alpha' component = transparency)
    # create the 'transparency layer'
alpha_dens = 0.*dens                       # all bins initially transparent
alpha_dens[np.where(dens>0)] = 1     # any bin with data is not transparent
dens_scale = dens/np.max(dens)*255         # scale to range [0,255] required for images
im_dens = np.array([dens_scale, dens_scale, dens_scale, alpha_dens])
implot = ax.imshow(im_dens.T, extent=[430, 475, 5899, 5950])    

ax.set_aspect('equal', adjustable='box')<jupyter_output><empty_output><jupyter_text>However, we can **do better** using something called a [Kernel Density Estimator](http://www.mglerner.com/blog/?p=28). This applies some smoothing to the histogram in an effort to combat arbitrary bin choice.I will use an implementation of [KDE](http://scikit-learn.org/stable/modules/density.htmlkernel-density-estimation) provided in a module called [scikit-learn](http://scikit-learn.org/stable/). This module provides a bunch of machine learning tools.<jupyter_code>f,ax = plt.subplots(1,1)                
f.set_size_inches(10,10)            
YlOrBr = matplotlib.cm.get_cmap('CMRmap_r')

# import the KDE
from sklearn.neighbors import KernelDensity
data = np.vstack([x,y]).T                    # format the data that we will use to train the KDE (the EQ locations)

# Set up a grid for eventual CONTOUR PLOTTING of the KDE
x0,x1,y0,y1 = [430, 475, 5899, 5950]                   # grid limits
xg = np.linspace(x0,x1,101)                            # 101 points in the x dir
yg = np.linspace(y0,y1,101)                            # same in the y
X, Y = np.meshgrid(xg,yg)                              # create a meshgrid, 2D arrays of x and y coords
xy = np.vstack([X.ravel(), Y.ravel()]).T               # combine the grid data together

# create the KDE - I will use a Gaussian kernal (that is, normally distributed location errors) with a bandwidth
# of 1 km (roughly the event location error)
kde = KernelDensity(bandwidth=1, kernel='gaussian', algorithm='ball_tree')   # set KDE properties
kde.fit(data)                                                                # train it with the EQ data
# now we need to evaluate the KDE at our grid locations
Z = np.exp(kde.score_samples(xy))
Z = Z.reshape(X.shape)
# plot contours of the density
levels = np.linspace(0, Z.max(), 11)                # the levels at which contours will be drawn
CS = ax.contourf(X, Y, Z, levels=levels, cmap=YlOrBr, zorder = 1)     # the (filled) contour method
plt.colorbar(CS,ax=ax)

# plot the reservoir outline over the top
im = plt.imread('groningen_reservoir.png')  
im[:,:,-1] = 1.*(im[:,:,0]<1.e-6)
implot = ax.imshow(im, extent=[430, 475, 5899, 5950], zorder=2)<jupyter_output><empty_output><jupyter_text>&lt;neat&gt; Cartopy The final aspect of 2D plotting to show is a handy little package called [cartopy](http://scitools.org.uk/cartopy/docs/latest/index.html). You'll have to install this package yourself (a good skill to have).Assuming you have installed Anaconda, go to the command line and type `conda install -c conda-forge cartopy`. Anaconda will automatically download and install cartopy and any packages it relies on. **Once you have done this, execute the cell below to plot a map of the world.**<jupyter_code># make cartopy available
import cartopy.crs as ccrs

# here, I am positioning the axes myself - this gives you much more flexibility in constructing figures with multiple subplots
ax1 = plt.axes([0.10, 0.15, 0.35, 0.7], projection = ccrs.PlateCarree())    # note the different projections      
ax2 = plt.axes([0.50, 0.15, 0.35, 0.7], projection = ccrs.Mollweide())          
    
# adding simple coastlines to each plot
for ax in [ax1,ax2]:
    ax.coastlines()

f = plt.gcf()
f.set_size_inches(20,20)<jupyter_output><empty_output><jupyter_text>I can zoom the plot in to show everyone's **second favourite** continent after [Zealandia](https://en.wikipedia.org/wiki/Zealandia), adding a little colour for splash.<jupyter_code># code credit: <NAME>
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

figsize = [7.48, 9.06]             # sets the dimensions of the figure to a full page
text_size = 10.                    # text size for axis ticks and annotations

# here we create the figure and set some its dimensions
    # fiddling with margins, extent is not strictly necessary, but does demonstrate the control you have
fig = plt.figure(figsize=figsize)
margin_x = [0.1, 0.02]
margin_y = [0.05, 0.02]
ax = plt.axes([margin_x[0], margin_y[0], 1-sum(margin_x), 1-sum(margin_y)], projection=ccrs.PlateCarree())
x_extent = [112, 154]              # longitude plot limits
y_center = -25                     # center of lattitude axis (limits computed to be consistent with lon extent) 
y_delta = figsize[1]/figsize[0]*(x_extent[1]-x_extent[0])
y_extent = [y_center-y_delta/2, y_center+y_delta/2]
ax.set_extent(x_extent+y_extent, ccrs.PlateCarree())

# add lines
ax.coastlines(resolution='50m')

# add coloured land and ocean
land_50m = cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='k',facecolor=cfeature.COLORS['land'])
ax.add_feature(land_50m, linewidth=.01)
ocean_50m = cfeature.NaturalEarthFeature('physical', 'ocean', '50m', edgecolor='face',facecolor=cfeature.COLORS['water'])
ax.add_feature(ocean_50m);<jupyter_output><empty_output><jupyter_text>We can do more!<jupyter_code># code credit: <NAME>
fig = plt.figure(figsize=figsize)
ax = plt.axes([margin_x[0], margin_y[0], 1-sum(margin_x), 1-sum(margin_y)], projection=ccrs.PlateCarree())
ax.set_extent(x_extent+y_extent, ccrs.PlateCarree())
ax.coastlines(resolution='50m')
ax.add_feature(ocean_50m)
ax.add_feature(land_50m, linewidth=.01)

# add state provinces 
states_provinces = cfeature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none', edgecolor='k')
ax.add_feature(states_provinces, linewidth=.5)

# add locations of two Enhanced Geothermal System projects
hb = {'lon': 140.45, 'lat':-27.495, 'txt':'Habanero'}
pl = {'lon': 139.7, 'lat':-30.2, 'txt':'Paralana'}
for EGS in [hb, pl]:
    # we'll use 'scatter' to add the points
    ax.scatter(EGS['lon'], EGS['lat'], marker = 's', c='r', zorder = 5, s = 3, lw=.1, edgecolors='k')
    # and 'text' to annotate them
    ax.text(EGS['lon']-0.5, EGS['lat'], EGS['txt'], fontsize = text_size, va='center', ha='right')

# here are some things we can do to pretty up the axes (again, nice to have control of, but not strictly necessary)
for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, alpha=0., linestyle = '--')
gl.xlabels_top = False
gl.ylabels_right = False
gl.xlines = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': text_size}
gl.ylabel_style = {'size': text_size}
ax.tick_params(direction='out', length=6, width=2)<jupyter_output><empty_output><jupyter_text>For more resources on the use of Cartopy, refer to the [gallery](http://scitools.org.uk/cartopy/docs/latest/gallery.html) and the [Natural Earth](http://www.naturalearthdata.com/) database.  &lt;/neat&gt;   2.3D - interactive view, meshed surfaces Plotting in 3D comes up less often than the plot types we've looked at so far. Perhaps that is due to difficulties finding a good **view angle** when looking at 3D data.However, there are times when 3D plots provide a really powerful visual aide. I'll just show some simple examples here.<jupyter_code># import tools for 3D axes
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib import pyplot as plt
import numpy as np

# create a grid
xg = np.linspace(0,1,31)              # evenly spaced grid points
yg = np.linspace(0,1,31)
ymin,ymax = [0.15,0.85]               # create a smaller subgrid in the y-dir for coloring
i1 = np.argmin(abs(yg-ymin))
i2 = np.argmin(abs(yg-ymax))
yg2 = yg[i1:i2+1]                     # subsample y coords
[X,Y] = np.meshgrid(xg,yg)            # create the two mesh grids
[X2,Y2] = np.meshgrid(xg,yg2)

# create a custom surface
    # parameters
xm = np.mean(xg)*0.8
ym = np.mean(yg)*1.2
sx = 0.02*3.
sy = 0.04*3.
    # function defining the surface in terms of x, y and parameters
def r(X,Y): 
    return (5-np.exp(-((X-xm)**2/sx+(Y-ym)**2/sy)))*(1-(X/4)**2)*(1+(Y/4)**2)

# create a figure with a 3D projection
fig = plt.figure(figsize=[15,8])
ax = fig.add_subplot(111, projection='3d')

# plot the function as a wireframe over the large grid
ax.plot_wireframe(X, Y, r(X,Y), lw = 0.5, color = 'k')
    # shade part of the wireframe according to the function value
CS = ax.plot_surface(X2, Y2, r(X2,Y2), rstride=1, cstride=1,cmap=cm.Oranges, lw = 0.5)
plt.colorbar(CS, ax=ax)<jupyter_output><empty_output><jupyter_text>It looks fine, but are we happy with that **view angle**? It's the **default** that Python chose - out of all possible angles, how can we expect Python to choose the **best**?It would be nicer if we could **modify it in real-time**...<jupyter_code># we'll add this IPython magic command
%matplotlib notebook

fig = plt.figure(figsize=[10,5])
ax = fig.add_subplot(111, projection='3d')

ax.plot_wireframe(X, Y, r(X,Y), lw = 0.5, color = 'k')
CS = ax.plot_surface(X2, Y2, r(X2,Y2), rstride=1, cstride=1,cmap=cm.Oranges, lw = 0.5)
plt.colorbar(CS, ax=ax)

# and also a show() command
plt.show()

# execute this cell and click and drag the subplot below until you find a view you are happy with
# HINT: depending on your Python version, you'll need to execute the cell twice OR restart the Python kernel (hotkey: 0 0) 
# and then execute the cell above and then this one (I don't know why this is the case...)<jupyter_output><empty_output><jupyter_text>(*This code also implemented in script: [2_visualisation/mesh_plot.py](mesh_plot.py) - try running it at the command line* `ipython mesh_plot.py`)  2.4D - interactive plots, movies Because the fourth dimension is **time**, right?Anyway, this section is about **moving** plots. There are two flavours we'll cover: interactive Jupyter Notebook plots with slider bars, and movies.  2.4.1 Interact Interactive plots are **excellent** in Jupyter notebooks. The `interact` method essentially **wraps around** any plot written as a **function**, adding slider bars, radio buttons, etc., that can be used to **update** plotting parameters in **real-time**. Let's look at the example of contamination due to a chemical spill into a groundwater aquifer with local flow velocity, $v$.\begin{equation}C(x,t) = \frac{M}{2\sqrt{\pi D t}}\exp\left(-\frac{(x-vt)^2}{4Dt}\right),\end{equation}where $C$ is concentration, $M$ is the mass of spilled contaminant, and $D$ is the diffusion coefficient.<jupyter_code>%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
# plot contamination down stream after two years
f,ax = plt.subplots(1,1)
x = np.linspace(0,10,101)                       # a vector of distance values to plot over
M = 1000.                                       # mass of chemical spilled
D = 0.1                                         # aquifer diffusivity
v = 1.                                          # groundwater velocity
t = 2.                                          # time at which to plot solution

# compute chemical concentration
C = M/(2*np.sqrt(np.pi*D*t))*np.exp(-(x-v*t)**2/(4*D*t))
ax.plot(x,C, 'r-')
ax.set_xlim([0, 10])
ax.set_ylim([0, 1.e3])
ax.set_xlabel('distance downstream [km]')
ax.set_ylabel('concentration [kg/m^3]')<jupyter_output><empty_output><jupyter_text>Let's turn it into a **function**. This is less difficult that you might think.<jupyter_code># a function that plots downstream contamination after two years, with flow velocity v as an input
def plot_contamination(v):
    f,ax = plt.subplots(1,1)
    x = np.linspace(0,10,101)
    M = 1000.
    D = 0.1
    t = 2.
    C = M/(2*np.sqrt(np.pi*D*t))*np.exp(-(x-v*t)**2/(4*D*t))
    ax.plot(x,C, 'r-')
    ax.set_xlim([0, 10])
    ax.set_ylim([0, 1.e3])
    ax.set_xlabel('distance downstream [km]')
    ax.set_ylabel('concentration [kg/m^3]')
    
# call the function with v = 1.5 to produce the same plot as above
plot_contamination(v=1.5)

# change the value of v and rerun the cell - does the plot update in a manner that you expect?<jupyter_output><empty_output><jupyter_text>Okay, now for the **magic**. We'll use `interact` from the ipywidgets module to turn this into an **interactive graph**.<jupyter_code># first, let's get the interact function
from ipywidgets import interact

# then its as simple as one line
interact(plot_contamination, v = (0.2, 4.8, 0.2))

# play with the slider bar below and see how the graph is redrawn in real time, each time with an updated value v
# the first input to 'interact' is the function to be executed
# the second input says 'create a slider for the parameter v, let it vary between 0.2 and 4.8 in increments of 0.2)<jupyter_output><empty_output><jupyter_text>&lt;homework&gt;<jupyter_code># CHANGE THE INTERACTIVE PLOT ABOVE to make D the slider parameter, varying between 0.02 and 0.98 in increments of 0.06

# **change the function so D is the input**
def plot_contamination(v):
    f,ax = plt.subplots(1,1)
    x = np.linspace(0,10,101)
    M = 1000.
    D = 0.1
    t = 2.
    C = M/(2*np.sqrt(np.pi*D*t))*np.exp(-(x-v*t)**2/(4*D*t))
    ax.plot(x,C, 'r-')
    ax.set_xlim([0, 10])
    ax.set_ylim([0, 1.e3])
    ax.set_xlabel('distance downstream [km]')
    ax.set_ylabel('concentration [kg/m^3]')
    
# **change the interact call**
interact();
# **to do**
# MAKE ANOTHER CHANGE so that BOTH v and D are slider parameters
# **your code here**<jupyter_output><empty_output><jupyter_text>&lt;/homework&gt; We'll come back to the interact function in the **next notebook** when we look at fitting and calibrating models.  2.4.2 Movies When you watch a movie, what does it actually comprise? A **sequence of still images**, flashed in front of you in sequence (usually about 24 each second).We shall take a similar approach to creating movies in Python. You can already create still images (that's what this whole notebook is about) so now we just need to extend this idea to creating a **series** of still images, each one different from the previous.Finally, we will need a way to **turn those images into a movie**. The example I'll show is for stress changes on a fault as a rupture propagates along it (an earthquake). Let's create a 5 second movie, in which the position of the rupture front, `a`, goes from 0.25 to 10 km.First, create a function that draws and saves a copy of the image.<jupyter_code>def plot_frame(i,a):
    ''' Plots a single frame of the fault rupture, at position a.
    
        i = frame number
        a = position of rupture front
        
        Stress changes are computed according to linear elastic fracture mechanics.
        s = s0 for x>>a         (far field stress)
        s = s0-ds for x<a       (stress drop inside rupture)
        s ~ K/sqrt(x) for x>a   (decaying stresses with singularity at crack tip)
    '''
    f,ax = plt.subplots(1,1)
    x = np.linspace(0,10.,1001)
    
    # calculate the stress field
    s0 = 25.                      # pre-earthquake stress
    ds = 3.                       # stress drop
    K = s0*np.sqrt(np.pi*a)       # stress intensity factor
    y = s0*np.ones(len(x))        # background stress vector
    y[np.where(x<a)] = s0-ds      # stress behind the rupture
    # stress increase in front of the rupture
    y[np.where(x>=a)] = s0+K/np.sqrt(2*np.pi*(x[np.where(x>=a)]-a)*1.e3) 
    
    # plot the stress changes
    ax.plot(x,y,'k-')                  # plot stress changes
    ax.plot([a,a], [s0-ds, 30], 'k:')  # stress singularity
    ax.set_xlim([x[0], x[-1]])
    ax.set_ylim([21, 30])
    ax.set_xlabel('position along fault / km')
    ax.set_ylabel('stress / MPa')
    
    # save the frame in a separate folder and name it using the index i
    plt.savefig('all_frames/frame{:04d}.png'.format(i), dpi=300)
    plt.close(f)                   # it will be important to close all these figures when we're looping
    
# this part should be familiar!
import os
if not os.path.isdir('all_frames'):
    os.makedirs('all_frames')

plot_frame(0,1.5)<jupyter_output><empty_output><jupyter_text>**Verify that a file called `frame0000.png` has been created in a directory `all_frames`.**   Now, let's create all the frames we'll need. We'll use a for loop to call `plot_frame` a bunch of times.<jupyter_code>FPS = 10           # ten frames per second, let's not be greedy
secs = 5           # total seconds
Nframes = secs*FPS  # total frames

# each frame has a different value a, lets create a vector of a values
avals = np.linspace(0.25, 10.0, Nframes)

# then loop over each avalue and create the frame
for i,a in enumerate(avals):               # enumerate is a handy function that both loops over the values in an array
    plot_frame(i,a)                        # AND gives you a corresponding index<jupyter_output><empty_output><jupyter_text>Once this cell has finished running, there should be a large number of figures in the directory `all_frames`. Last step is to 'glue' them all together. We'll use an external command line tool called [ffmpeg](https://www.ffmpeg.org/documentation.html). (Note, if you downloaded this course from GitHub, you'll need to download a copy of ffmpeg.exe and drop it into the `module 2 - visualisation` folder.)This isn't a course on use of ffmpeg. I have provided the Python code to call out to ffmpeg to create a video, but if you want to get fancy you'll need to do your own Googling.**Execute the cell below to create a short movie. (This won't work if you're on Microsoft Azure)**<jupyter_code># we'll use os.system - a handy command that let's us implement command line calls without returning to the command line
os.system('ffmpeg -framerate {:d} -i all_frames/frame%04d.png earthquake_movie.mp4'.format(FPS))<jupyter_output><empty_output><jupyter_text>**Verify that the movie you've created corresponds to the plotted frames.** Last job - tidy your bedroom!<jupyter_code>from glob import glob
fls = glob('all_frames/*.png')
for fl in fls: os.remove(fl)
os.rmdir('all_frames')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11682605324373223, 'fraction_numerical': 0.023546135952442492, 'mean_word_length': 3.587502964192554, 'pattern_counts': {'":': 0, '<': 170, '<?xml version=': 0, '>': 170, 'https://': 6, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 41, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '39', 'dir': 'jupyter-structured-clean-dedup', 'id': '335121', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 13176, 'n_tokens_neox': 12245, 'n_words': 4450, 'path': '2_visualisation/Visualisation.ipynb', 'repo_name': 'ddempsey/python_for_geoscientists', 'size': '38708'}