Exercise 2

In the second exercise you need to plot ocean depth against ocean age. You have the topography / bathymetry data and the age data. You should do this for the "deep oceans" to avoid including the continental shelves and, for that matter, the continents. The age dataset is only valid where there are magnetic stripes. Other places in the grid are flagged as "Not a Number" or NaN and will not plot. You might have to strip out these points for plotting / curve fitting though.


In [2]:
%pylab inline

from osgeo import gdal

import cartopy.crs as ccrs
import matplotlib.pyplot as plt


Populating the interactive namespace from numpy and matplotlib

Hint

You will need to have the two datasets at the same resolution and valid data shape. The data are stored in an image so re-interpolation can be done with the image resizing fanciness available in ... (guess where !)


In [11]:
from scipy.misc import imresize

etopoH = gdal.Open("../../Data/Resources/ETOPO1_Ice_c_geotiff.tif") 
etopoH_img = etopoH.ReadAsArray()[::4,::4]
del(etopoH)

ages = np.load("../../Data/Resources/global_age_data.3.6.z.npz")["ageData"]

etopoH_1 = imresize(etopoH_img, ages.shape, interp='bilinear', mode="F")
etopoH_1[ np.isnan(ages) ] = np.nan

print etopoH_img.shape
print etopoH_1.shape
print ages.shape


(2700, 5400)
(1801, 3601)
(1801, 3601)

Task

Make an image plot of etopoH_1 and ages (check they are the same shape and same NaN mask)


In [ ]:
## Code here

In [ ]:
## Code here

Task

Make a scatter plot of the depth against age. You might need to downsample the number of points as there are a lot of them:

ages_decimated = ages[::10]
topo_decimated = etopoH[::10]

Do you see any systematic trend ?

Hint: don't worry if you do not


In [ ]:
## Code here

Smoothing

Now try interpolating the bathymetry to a smoothed / downsampled version of the ages.

ages_reduced = ages[::8,::8] # This is arbitrary, you should try some different ones

# OR

ages_reduced2 = imresize(ages, (226,451), interp='bilinear', mode="F")

# OR

ages_reduced2 = imresize(ages, (226,451), interp='bicubic', mode="F")

# OR

ages_reduced2 = imresize(ages, (226,451), interp='lanczos', mode="F")

Comment on which of these is the more effective.


In [34]:
ages_reduced  = ages[::8,::8] # This is arbitrary, you should try some different ones
plt.imshow(ages_reduced)
plt.show()

ages_reduced2 = imresize(ages, (226,451), interp='bicubic', mode="F")
plt.imshow(ages_reduced2)
plt.show()


## Have a go ...

Task

Plot the age/depth data and see if it looks smoother.

Make a scatterplot with this data too (note: you will have to downsample the bathymetry too)

Does smoothing help bring out a trend ?


In [ ]:
### Your code here

Task

Fit $ \textrm{depth} = A + B\sqrt{\textrm{age}} $ to this data and create a plot

You should be able to use

from scipy.optimize import curve_fit
help(curve_fit)

In [ ]:
### Your code here

Discussion

Obviously there are problems with the data - looking at every pixel in the image does not account for regions where there are sea mounts or other features on the ocean floor that we might consider anomalous. Actually, the other reason this doesn't work very well is that it fails to account for sediment accumulation and loading.