Maps of useful dimensional and non-dimensional parameters for the coastal world ocean


In [167]:
%matplotlib inline
%load_ext secnum
%secnum

import numpy as np
import scipy as sp
import h5py
import netCDF4 

import sys
sys.path.append('/media/data/Work/tools/python/dc/')
from dcutils import *

import numexpr as ne
from scipy.io import loadmat
import scipy.interpolate

from scipy import ndimage
from skimage import filter
import skimage.morphology as morph
import skimage.segmentation as seg

from datetime import datetime
import time
from IPython import display
import pylab

import seawater as sw

# matplotlib section
import matplotlib as mpl
import matplotlib.pyplot as plt
# big figures
mpl.rcParams['savefig.dpi'] = 2 * mpl.rcParams['savefig.dpi']
# better color cycling
import brewer2mpl
colors =  brewer2mpl.get_map('Set2', 'qualitative', 8).mpl_colors
mpl.rcParams['axes.color_cycle'] = colors


The secnum extension is already loaded. To reload it, use:
  %reload_ext secnum

Reading data

Read in data files - etopo for bathymetry, WOA to calculate stratification


In [168]:
# load ETOPO2 bathymetry - convert to ETOPO1 eventually
bathyfile = 'data/etopo2_extract.mat'
bathy = loadmat(bathyfile, struct_as_record=False, squeeze_me=True)
topo = bathy['topo'] # pull out the topo structure

# pull out MAB for testing
ix = [find_nearest(topo.x, -75), find_nearest(topo.x, -60)]
iy = [find_nearest(topo.y, 35), find_nearest(topo.y, 45)]

#xx = topo.x[ix[0]:ix[1]]
#yy = topo.y[iy[0]:iy[1]]
#zz = topo.z[ix[0]:ix[1], iy[0]:iy[1]]

xx = topo.x
yy = topo.y
zz = topo.z

# land mask
landmask = zz > 0
# deep water mask
deepmask = zz < -4000

In [169]:
# dx,dy need to be converted from dlon,dlat to metres
dx = np.diff(xx).mean()*110000;
dy = np.diff(yy).mean()*110000;

# calculate gradients assuming dx=dy=1
hy, hx = np.gradient(zz)
hxy, hxx = np.gradient(hx)
hyy, _ = np.gradient(hy)

# get actual gradient
hx /= dx
hy /= dy
hxy /= (dx * dy)
hxx /= (dx * dx)
hyy /= (dy * dy)

# apply land mask
slope = np.ma.masked_array(np.sqrt(hx**2 + hy**2), 
                           landmask)

# get Gaussian curvature
curv = np.abs(np.ma.masked_array((hxx * hyy - (hxy **2)) / 
                          (1 + (hx **2) + (hy**2)) ** 2, 
                          landmask))

In [170]:
# plot
fig = plt.figure()
ax = fig.gca()
handle = ax.pcolorfast(xx, yy, zz.T, cmap='OrRd_r')
fig.colorbar(handle)
handle.set_clim((-5000,0))
fig.hold(True)
ax.contour(xx, yy, zz.T, levels=[0,-200], colors='k')
ax.set_title('Topo')

fig = plt.figure()
ax = fig.gca()
handle = ax.pcolorfast(xx, yy, slope.T, cmap='OrRd')
fig.colorbar(handle)
fig.hold(True)
ax.contour(xx, yy, zz.T, levels=[0,-200], colors='k')
ax.set_title('Magn. of bathymetric slope')

fig = plt.figure()
ax = fig.gca()
handle = ax.pcolorfast(xx, yy, 
                       np.ma.masked_array(filter.prewitt(zz).T, landmask.T), 
                       cmap='OrRd')
fig.colorbar(handle)
fig.hold(True)
ax.contour(xx, yy, zz.T, levels=[0,-200], colors='k')
ax.set_title('Magn. of bathymetric slope')


Out[170]:
<matplotlib.text.Text at 0x2556cfd0>