#!/usr/bin/env python # coding: utf-8 # In[1]: import numpy as np import matplotlib.pyplot as plt get_ipython().run_line_magic('matplotlib', 'inline') from IPython.display import Image import pandas as pd # ### Lecture 18: # # - start to make some basic maps using **Cartopy**. Yippee (we love maps). # # # # # # ## Introduction to maps # # Maps in Python used to be plotted using the tools in **matplotlib**'s **Basemap** module. But **Basemap** is being deprecated, so we are switching to a new, actively developed toolkit called **cartopy**. # # The first thing we have to do is import a bunch of stuff from **cartopy** and from matplotlib (assuming you have **cartopy** installed - if not, install at least cartopy=0.17.0) by typing this on your command line: # # conda install cartopy=0.17.0 # # Here's what we need to make maps with **cartopy**: # In[2]: 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 # There are many different types of maps used in the Earth Sciences. A map tries to represent something that is essentially 3D (a globe) onto a 2D medium (paper or a computer screen). So all maps, except those at the smallest scale, will increasingly distort the area as the scale increases because the Earth is not 2-dimensional. [No the Earth is not flat! https://en.wikipedia.org/wiki/Modern_flat_Earth_societies] # # When we choose a map projection, we seek the one that distorts the least for our purpose. # # Here we will use a few popular projections to make maps on both the global and local scale. Let's start with the global projections. # # ### Mercator Projection # The one probably everyone is familiar with is the standard Mercator Projection: # In[3]: Image(filename='Figures/mercator.jpg',width=500) # In **cartopy** maps are instances of the **Axes** figure class. They have many methods, for example outlining continents, national boundaries, and state boundaries and plotting geospatial data, such as sampling locations, earthquakes, and much else. # # Before we plot can plot a set of coordinates, they must be transformed from latitudes and longitudes to map coordinates and then plotted like anything else in **matplotlib**. # # To make a map instance, we use the **projection** keyword in the **plt.axes( )** object (normally **ax**). # To set up the map, you need to know what type of projection you want (there are a growing number to choose from) and depending on the projection type, you need to know the map boundaries and other particulars. # # We are starting with the Mercator projection. To do this, we set the projection keyword to '**ccrs.Mercator**'. # # # In[4]: help(ccrs.Mercator) # So we need: **central_longitude, min_latitude**, and **max_latitude** # # # We make the map instance with a call like: # # **ax = plt.axes(projection=ccrs.Mercator(central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None))** # # # We also have to put something on the map, so let's draw the coastlines by using the **coastlines()** method: # # **ax.coastlines()** # # Here's our basic Mercator map of the whole world: # # # In[5]: ax = plt.axes(projection=ccrs.Mercator( central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None)) ax.coastlines(); # We usually want lines of latitude (parallels) and longitude (meridians) on there as well, so we use the **gridlines( )** method with the keyword **crs=ccrs.PlateCarree( )** to make a **globe** object (usually called $g$ or $g1$. # In[6]: help(ax.gridlines) # Now we stick the latitudes and longitudes onto our globe object ($g1$ in the following example). We will need the **matplotlib.ticker** module (already imported as **mticker**) which has the useful attributes **mticker.FixedLocator( )** method which allows us to make an create a special object of latitudes or longitudes with the desired spacing (using our old friend **np.arange( )**). These can be attached to the figure object using figure attributes **g1.ylocator** and **g1.xlocator** for latitudes and longitudes respectively. (You'll see how this works soon, so be a little patient!). # # The steps are: # # 1) Make the figure object ($ax$) with the desired projection using **plt.axes( )**. # # 2) Put on the coastlines # # 3) Make the globe object ($g1$) using **ax.gridlines( )** # # 4) Put on the latitude and longitudes by setting the **g1.ylocator** and **g1.xlocator** attributes to the desired **mticker.FixedLocator(ARRAY)** object. # # Here we go: # In[7]: # make the figure object ax = plt.axes(projection=ccrs.Mercator( central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None)) # put on the coastlines ax.coastlines() # make the globe object gl=ax.gridlines(crs=ccrs.PlateCarree(),linewidth=2,linestyle='dotted') # put on the latitudes gl.ylocator=mticker.FixedLocator(np.arange(-90,91,30)) # put on the longitudes gl.xlocator=mticker.FixedLocator(np.arange(-180.,181.,60.)); # # And surely you want the lat/long labels to show, so you can modify your commands like this: # In[8]: ax = plt.axes(projection=ccrs.Mercator( central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None)) ax.coastlines() g1=ax.gridlines(crs=ccrs.PlateCarree(),linewidth=2,linestyle='dotted',draw_labels=True) g1.ylocator=mticker.FixedLocator(np.arange(-90,91,30)) g1.xlocator=mticker.FixedLocator(np.arange(-180.,181,60.)) g1.xformatter = LONGITUDE_FORMATTER g1.yformatter = LATITUDE_FORMATTER # It seems that you can either have a gap in the latitude lines (as in this example), or a doubled 180$^{\circ}$W/180$^{\circ}$E. This is a "feature". But I did find a workaround for this special case. You draw the gridlines separately from the labels. We don't get the nice degree labels, but at least we don't get the doubled 180 tick mark! # # In[9]: ax = plt.axes(projection=ccrs.Mercator( central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None)) ax.coastlines() ylabels=np.arange(-90,91,30) xlabels=np.arange(-180.,181,60.) g1=ax.gridlines(xlocs=xlabels,ylocs=ylabels, crs=ccrs.PlateCarree(),linewidth=2,linestyle='dotted',draw_labels=False) xlabels=np.arange(-180.,180,60.) glabels = ax.gridlines(xlocs=xlabels, ylocs=ylabels, draw_labels=True, alpha=0) # To turn labels off for the top, bottom, left or right sides, we can use attributes of the globe object like this: # # In[10]: ax = plt.axes(projection=ccrs.Mercator( central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None)) ax.coastlines() ylabels=np.arange(-90,91,30) xlabels=np.arange(-180.,181,60.) g1=ax.gridlines(xlocs=xlabels,ylocs=ylabels, crs=ccrs.PlateCarree(),linewidth=2,linestyle='dotted',draw_labels=False) xlabels=np.arange(-180.,180,60.) glabels = ax.gridlines(xlocs=xlabels, ylocs=ylabels, draw_labels=True, alpha=0) glabels.xlabels_top = False glabels.ylabels_right = False # # # Now maybe you want some color? You can set the color of the continents with: # # **ax.add_feature(LAND,color='orange')** # # and color in the oceans with: # # **ax.add_feature(OCEAN,color='lightblue')** # # # In[11]: ax = plt.axes(projection=ccrs.Mercator( central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None)) ylabels=np.arange(-90,91,30) xlabels=np.arange(-180.,181,60.) g1=ax.gridlines(xlocs=xlabels,ylocs=ylabels, crs=ccrs.PlateCarree(), linewidth=2,linestyle='dotted',draw_labels=False) xlabels=np.arange(-180.,180,60.) glabels = ax.gridlines(xlocs=xlabels, ylocs=ylabels, draw_labels=True, alpha=0) glabels.xlabels_top = False ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.coastlines(); # this has to come last, or it is buried by the ocean, land colors # The lakes are missing! You can change lake colors using the **LAKES** attribute if you so desire. You can set the facecolor and edgecolor separately. # # You also want to put something ON the map? No problem. Just take your latitude and longitude (or arrays of them) and convert them to map coordinates using the **transform=ccrs.Geodetic( )** keyword in the **plt.plot( )** command. These you plot just like any ordinary **matplotlib** plot, except you need a **ax.set_global( )** command after the plot, or you get a map near San Diego (for some reason...). # # Let's plot the position of San Diego as a big white star with the marker="\*", markersize=20, and color = 'white' arguments. # # In[12]: ax = plt.axes(projection=ccrs.Mercator( central_longitude=180.0, min_latitude=-70.0, max_latitude=70.0, globe=None)) ylabels=np.arange(-90,91,30) xlabels=np.arange(-180.,181,60.) g1=ax.gridlines(xlocs=xlabels,ylocs=ylabels, crs=ccrs.PlateCarree(), linewidth=2,linestyle='dotted',draw_labels=False) xlabels=np.arange(-180.,180,60.) glabels = ax.gridlines(xlocs=xlabels, ylocs=ylabels, draw_labels=True, alpha=0) glabels.xlabels_top = False ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,facecolor='lightblue',edgecolor='black') ax.coastlines() San_lat=33 San_lon=-117%360 # takes the west longitude and converts to 0=>360 ax.plot([San_lon],[San_lat],marker='*',color='white',\ markersize=20,transform=ccrs.Geodetic(),markeredgecolor='black') ax.set_global() # we need this or weird things happen ax.coastlines(); # ### Orthographic Projection # # The Mercator is a nice classical map, but it sure does distort the map at high latitudes. Think back to the lecture on the hypsometric curves... # # Another type of map projection is the orthographic projection which is much less distorted. The downside to this projection is that you cannot see the whole globe at once. To create an orthographic map, you initialize a map instance with the projection set to ccrs.Orthographic and a tuple of the central longitude and latitude. # In[13]: ax = plt.axes(projection=ccrs.Orthographic(-75, 42)) ax.coastlines(); # And finish the map as before. # In[14]: ax = plt.axes(projection=ccrs.Orthographic(-75, 42)) San_lat=33 San_lon=-117%360 # takes the west longitude and converts to 0=>360 using modular syntax 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.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,facecolor='lightblue',edgecolor='black') ax.plot([San_lon],[San_lat],marker='*',color='white',\ markersize=20,transform=ccrs.Geodetic(),markeredgecolor='black') ax.set_global() ax.coastlines(); # ### Mollweide projection # # One more global scale example of a map projection is the Hammer projection (one of my favorites). Unfortunately, **cartopy** does not have a Hammer projection (yet), so we will learn about a Mollweide projection instead. And really it is pretty similar, so no worries. The Mollweide projection is always a global map centered on the equator, so all you need to specify is the central longitude (with the **central_longitude** keyword). Note that any central_longitude but 180 is super slow for some reason... # In[15]: #no Hammer in cartopy but yes mollweide ax = plt.axes(projection=ccrs.Mollweide(central_longitude=180)) 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(-90,91,30)) gl.xlocator=mticker.FixedLocator(np.arange(0,400,30)); gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,facecolor='lightblue',edgecolor='black') ax.plot([San_lon],[San_lat],marker='*',color='white',\ markersize=20,transform=ccrs.Geodetic(),markeredgecolor='black') ax.set_global() ax.coastlines(); # ### Lambert conformal # # The maps we've explored so far are well and good for global scale problems, for example plotting the locations of earthquakes around the globe, but not so great for more local problems, like a map of sampling sites. For this we need a smaller scale map and the Lambert confomal conic projection is a popular choice. For this we must specify the central latitude and longitude and the map boundaries (with **ax.set_extent( )**. # # Here's what we know so far: # In[16]: # cannot get labels for this projection yet proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33) ax = plt.axes(projection=proj) ax.set_extent([-130, -70, 20, 52], crs=ccrs.PlateCarree()) ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,facecolor='lightblue',edgecolor='black',linewidth=1) ax.plot([San_lon],[San_lat],marker='*',color='white',\ markersize=20,transform=ccrs.Geodetic(),markeredgecolor='black') gl=ax.gridlines(ylocs=np.arange(0,90,15.),xlocs=np.arange(-180.,180,15.),\ linewidth=2, linestyle="dotted") ax.coastlines(); # One problem is that **cartopy** does not support tick labels in the Lambert projection. I found this work around online (https://gist.github.com/ajdawson/dd536f786741e987ae4e) and let's see if it works. # In[17]: import shapely.geometry as sgeom from copy import copy def find_side(ls, side): """ Given a shapely LineString which is assumed to be rectangular, return the line corresponding to a given side of the rectangle. """ minx, miny, maxx, maxy = ls.bounds points = {'left': [(minx, miny), (minx, maxy)], 'right': [(maxx, miny), (maxx, maxy)], 'bottom': [(minx, miny), (maxx, miny)], 'top': [(minx, maxy), (maxx, maxy)],} return sgeom.LineString(points[side]) def lambert_xticks(ax, ticks): """Draw ticks on the bottom x-axis of a Lambert Conformal projection.""" te = lambda xy: xy[0] lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te) ax.xaxis.tick_bottom() ax.set_xticks(xticks) ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels]) def lambert_yticks(ax, ticks): """Draw ricks on the left y-axis of a Lamber Conformal projection.""" te = lambda xy: xy[1] lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te) ax.yaxis.tick_left() ax.set_yticks(yticks) ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels]) def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor): """Get the tick locations and labels for an axis of a Lambert Conformal projection.""" outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist()) axis = find_side(outline_patch, tick_location) n_steps = 30 extent = ax.get_extent(ccrs.PlateCarree()) _ticks = [] for t in ticks: xy = line_constructor(t, n_steps, extent) proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1]) xyt = proj_xyz[..., :2] ls = sgeom.LineString(xyt.tolist()) locs = axis.intersection(ls) if not locs: tick = [None] else: tick = tick_extractor(locs.xy) _ticks.append(tick[0]) # Remove ticks that aren't visible: ticklabels = copy(ticks) while True: try: index = _ticks.index(None) except ValueError: break _ticks.pop(index) ticklabels.pop(index) return _ticks, ticklabels # In[18]: proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33) fig = plt.figure(figsize=(6,6), frameon=True) # you need this frameon to be true ax = plt.axes(projection=proj) ax.set_extent([-130, -70, 20, 52], crs=ccrs.PlateCarree()) ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,facecolor='lightblue',edgecolor='black') ax.plot([San_lon],[San_lat],marker='*',color='white',\ markersize=20,transform=ccrs.Geodetic(),markeredgecolor='black') # *must* call draw in order to get the axis boundary used to add ticks: fig.canvas.draw() xticks=list(range(-180,-30,15)) yticks=list(range(0,90,15)) ax.gridlines(ylocs=yticks,xlocs=xticks,linewidth=2, linestyle="dotted") ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER) # you need this here ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)# you need this here, too lambert_xticks(ax, xticks) lambert_yticks(ax, yticks) ax.coastlines(); # Oh my, that is lovely! # But it sure would be nice to put on national and state boundaries. We can draw the national boundaries in a thick line like this: # In[19]: proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33) fig = plt.figure(figsize=(6,6), frameon=True) # you need this frameon to be true ax = plt.axes(projection=proj) ax.set_extent([-130, -70, 20, 52], crs=ccrs.PlateCarree()) ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,color='lightblue',linewidth=1) ax.plot([San_lon],[San_lat],marker='*',color='white',\ markersize=20,transform=ccrs.Geodetic(),markeredgecolor='black') # *must* call draw in order to get the axis boundary used to add ticks: fig.canvas.draw() xticks=list(range(-180,-30,15)) yticks=list(range(0,90,15)) ax.gridlines(ylocs=yticks,xlocs=xticks,linewidth=2, linestyle="dotted") ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER) # you need this here ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)# you need this here, too lambert_xticks(ax, xticks) lambert_yticks(ax, yticks) ax.coastlines(); ax.add_feature(BORDERS,linestyle='-',linewidth=2); # And finish off with the state boundaries like this: # In[20]: proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33) fig = plt.figure(figsize=(6,6), frameon=True) # you need this frameon to be true ax = plt.axes(projection=proj) ax.set_extent([-130, -70, 20, 52], crs=ccrs.PlateCarree()) ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,color='lightblue',linewidth=1) ax.plot([San_lon],[San_lat],marker='*',color='white',\ markersize=20,transform=ccrs.Geodetic(),markeredgecolor='black') # *must* call draw in order to get the axis boundary used to add ticks: fig.canvas.draw() xticks=list(range(-180,-30,15)) yticks=list(range(0,90,15)) ax.gridlines(ylocs=yticks,xlocs=xticks,linewidth=2, linestyle="dotted") ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER) # you need this here ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)# you need this here, too lambert_xticks(ax, xticks) lambert_yticks(ax, yticks) ax.coastlines(); ax.add_feature(BORDERS,linestyle='-',linewidth=2) states_provinces = cfeature.NaturalEarthFeature( category='cultural', name='admin_1_states_provinces_lines', scale='50m', edgecolor='purple', facecolor='none', linestyle='dotted') ax.add_feature(states_provinces); # You might also want to know that there are other options for background colors but **cartopy** has a very limited set: # the stock image, several Blue Marble version and a shaded relief map: # In[21]: proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33) ax = plt.axes(projection=proj) ax.set_extent([-130, -70, 20, 52], crs=ccrs.PlateCarree()) ax.stock_img(); # For more on background images, see this website: http://earthpy.org/cartopy_backgroung.html # And now for some science! # # # ### Earthquake locations - last 5 years # # Remember the data we've been using with latitude and longitude coordinates? Now we can place those on our maps! Here are all the earthquakes that have occured over a five year period with magnitudes greater than or equal to 5. The data come from: http://earthquake.usgs.gov/earthquakes/search/ # # Let's take a look: # # In[25]: open('Datasets/EarthquakeLocations/last5Years.csv').readlines()[0:5] # So we have to skip the first row (by setting the keyword **skiprows** to 1). # In[26]: eq_data=pd.read_csv('Datasets/EarthquakeLocations/last5Years.csv',skiprows=1) eq_data.columns # There is a lot of information in this file, but the most important (for now) is the latitude and longitude. We can deal with the depth and magnitude later. # # For now, let's extract that location information into arrays. # # Remember, this is how you get stuff out of a Pandas **DataFrame** into **NumPy** arrays: # # In[27]: lats=eq_data.latitude.values #make an array of the values of a pandas Series lons=eq_data.longitude.values lons=lons print (lats) print (lons) # Now we can plot these on a Mollweide projection. Everything is just as we did before, except I'm putting on the earthquake locations as white dots. The size of the dots can be set with the keyword **markersize**. # In[28]: plt.figure(1,(8,8)) # make a big figure (8x8) ax = plt.axes(projection=ccrs.Mollweide(central_longitude=180)) gl=ax.gridlines(crs=ccrs.PlateCarree(),color='black',linewidth=1,linestyle='dotted') gl.xlabels_top = False gl.ylocator=mticker.FixedLocator(np.arange(-90,91,30)) gl.xlocator=mticker.FixedLocator(np.arange(0,400,30)); gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,color='lightblue',linewidth=1) ax.plot([lons],[lats],marker='o',color='white',\ markersize=4,transform=ccrs.Geodetic(),markeredgecolor='black') ax.set_global() ax.coastlines(); # And did you notice how most of the earthquakes fall in narrow bands? Those are the plate boundaries you have have heard so much about of course. # # You might have also heard that earthquakes at trenches (like around the Pacific ocean's 'ring of fire') get deeper in a systematic way and reveals the location of the downgoing slabs. # # Included with this dataset, is the depth of the earthquakes, so we could color code the dots by depth and try to see the pattern. # # First let's see the range and frequency of different depths. We can plot them in a histogram to take a quick peek. # In[29]: plt.hist(eq_data.depth.values,bins=30,density=True) plt.xlabel('Depth (km)') plt.ylabel('Frequency'); # The majority of earthquakes occur in the top 30-50 km but they continue all the way down to 650 km. So let's make some bins to group the data by depth. We'll then color code the earthquakes, like the visible light spectrum, with the deepest earthquakes plotted as red. # # In[30]: depths=[33,70,150,300,400,650] # a list of depth bins colors=['purple','blue','lightblue','green','yellow','orange','red'] # Putting it all together: # In[31]: # this is just like before except for the background color (now white) plt.figure(1,(8,8)) # make a big figure (8x8) ax = plt.axes(projection=ccrs.Mollweide(central_longitude=180)) gl=ax.gridlines(crs=ccrs.PlateCarree(),color='black',linewidth=1,linestyle='dotted') gl.xlabels_top = False gl.ylocator=mticker.FixedLocator(np.arange(-90,91,30)) gl.xlocator=mticker.FixedLocator(np.arange(0,400,30)); gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER ax.add_feature(OCEAN,color='lightblue') ax.add_feature(LAND,color='orange') ax.add_feature(LAKES,color='lightblue',linewidth=1) # left out some fiddly bits to keep the map simple # now we put on the earthquakes last_bin=0 # this is the lower bound of the first bin for depths # or the upper bound of the last bin (same thing) for d in depths: # step through the depths list # use Pandas filtering to fish out depths in this range depth=eq_data[(eq_data.depth=last_bin)] # use the filtering of Pandas # convert to map coordinates # put this batch of earthquakes on the map plt.plot([depth.longitude.values],[depth.latitude.values],\ transform=ccrs.Geodetic(),marker='o',color=colors[depths.index(d)],markersize=5) # replace the lower depth bound with the upper one last_bin=d # increment the lower bound for the bin ax.set_global() ax.coastlines(); # Well, at least we see that the ridges all have shallow earthquakes and the deepest earthquakes are the farthest from the trench - to really see the so-called 'Benioff zones' we should use 3D plots and zoom in on the area of interest. We will learn 3D plotting tricks soon, so be patient. :) # # # # Assignment #6 # - Rename the notebook with the format LastnameFirstInitial_HomeworkNumber. For example, **CychB_6** # # # ### 1. Polynomial fits # # Certain isotopes are unstable and decay through _radioactive decay_. The formula for radioactive decay is: # # $N= N_o \exp^{-\lambda T}$ # # $\lambda$ is the decay constant (the time for $N$ to decay to $1/\exp$ of the original value # # $T$ is time # # $N$ is the number of nuclei remaining after time $T$ # # $N_o$ is the original number or parent nuclei # # The _half-life_ ($t_{1/2}$) is the time for $N$ to decay to $N_o$/2. The formula is: # # $t_{1/2} = \frac {ln 2}{\lambda}$ # # - Write a lambda function to calculate radioactive decay. The decay constant, time, and initial parent should be supplied as parameters. # - The half-life of radiocarbon is 5,730 yrs. Calculate the decay constant of radiocarbon # - Assume that the the initial parent, $N_o$, is 1, and time ranges from 0 to 7 half-lives. Use your function to calculate radioactive decay for radiocarbon # - Plot a curve of $N$ versus $T$ # - Calculate the best-fit polynomial to your curve. Use both **np.polyfit( )** and the **sklearn** way. Try different degrees for you polynomial, 3 and 13. # - Draw a red vertical line at the half-life of radiocarbon and a red horizontal line at .5. # # # ### 2. Seaborn pairplots # # - Argon-argon dating is one of the primary ways that geologists determine the age of their rocks. One interesting study attempted to date the age of hominid footprints in Africa (Liutkus-Pierce et al., 2016; http://dx.doi.org/10.1016/j.palaeo.2016.09.019) around Lake Natron in northern Tanzania. By using different methods and isotopic systems they dated the footprints as being older than 5,760 kyr, but younger than 19.1 kyr. As part of the study, the authors dated biotites and hornblends from the underlying substrate in an attempt to provide the older bound. They used three different argon-argon methods on six different samples. The data are in the file 'Datasets/radioisotopicAgeTanzania.xlsx' and the results are in the columns labelled 'Integ. Age', 'Plat. Age' and 'Isochron Age'. The ages are in They concluded that the ages were "too old", but the data set is non-the-less interesting for our purposes. # - read the data into a Pandas DataFrame. # - drop the rows that do not have all three methods. # - make a seaborn pair plot of the three different dating methods, coloring the points by sample name. # - What is the most reasonable guess of the age and why would this be "too old". # # ### 3. Maps # # People have been running observatories that measure the geomagnetic field since the mid-1800s and the US Geological Survey maintains a number of them throughout the US and US territories. The locations of the US array are available at this website: https://www.usgs.gov/natural-hazards/geomagnetism/science/observatories?qt-science_center_objects=0#qt-science_center_objects. The file that you download is a .kmz file suitable for importing into Google Earth, but we have translated the locations to a .csv file that can be loaded into a Pandas DataFrame. # - Read the data file "Datasets/geomagnetic_observatories.csv" into a pandas DataFrame # - Make a 10x10 figure object. # - Make an orthographic projection that shows all the locations of the observatories as red circles. # - Put on the gridlines as in the lecture. # - Label the points with the code of the observatory name. Make the text be on top and to the right of the dots. # # # # Your code must be fully commented. # # Turn in your notebook on the Canvas website # # You will recieve a "zero" if the assignment is late, copied from someone else, or not fully commented. If the degree of copying is serious enough, you WILL be reported, so just don't do it. # # If you have trouble - contact the TA or the instructor - we are here to help! # In[ ]: