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>

Thresholding this image (limiting the max) returns a pretty sharp edge. I think I can use this and some existing edge detection algorithm.


In [171]:
f = plt.figure()
ax = f.gca();
handle = ax.pcolorfast(xx, yy, curv.T *~deepmask.T, cmap='OrRd')
f.colorbar(handle)
handle.set_clim((0,1e-11))
fig.hold(True)
ax.contour(xx, yy, zz.T, levels=[0,-200], colors='k')
plt.title('Gaussian topography curvature')


Out[171]:
<matplotlib.text.Text at 0x308c8d50>

Gaussian curvature looks noisy and there's a question of extracting the contour from it.


In [172]:
thresh = filter.threshold_otsu(zz)
print(thresh)


-2158

Shelfbreak depth / contour & Slope Width


In [195]:
# based on http://scikit-image.org/docs/dev/auto_examples/plot_canny.html
# create image
#im = (slope > 0.004)

gradient_thresh = 5e-3

gradient = filter.sobel(zz, mask=(~landmask))
im = gradient > gradient_thresh;
im_morph = morph.binary_closing(im, morph.square(8))
#im_morph = morph.dilation(slope*~landmask*~deepmask, morph.square(10))
edges = seg.find_boundaries(im_morph)

zzshallow = np.ma.masked_array(zz.T, ~deepmask)
filtered = np.ma.masked_array(zz.T, abs(im_morph.T-1))
edge_topo = np.ma.masked_array(zz.T, abs(edges.T))

f = plt.figure()
a = f.gca()
plt.imshow(gradient.T)
plt.colorbar()
a.invert_yaxis()
a.hold(True)
a.contour(zz.T, levels=[0,-200], colors='w')
a.set_title('Sobel gradient')

f = plt.figure()
a = f.gca()
plt.imshow(im.T)
plt.colorbar()
a.invert_yaxis()
a.hold(True)
a.contour(zz.T, levels=[0,-200], colors='w')
a.set_title('Thresholded gradient mask')

f = plt.figure()
a = f.gca()
plt.imshow(im_morph.T)
plt.colorbar()
a.invert_yaxis()
a.hold(True)
a.contour(zz.T, levels=[0,-200], colors='w')
a.set_title('Processed thresholded gradient mask')

fig = plt.figure()
ax = fig.gca()
hsb = ax.pcolorfast(xx, yy, filtered, cmap='RdYlBu')
ax.hold(True)
ax.contour(xx, yy, zz.T, levels=[0,-200], colors='k')
fig.colorbar(hsb)
#ax.set_xlim((-75, -60))
#ax.set_ylim((35, 45))
hsb.set_clim((-6000,0))
ax.set_title('topo * processed mask')

fig = plt.figure()
ax = fig.gca()
hsb = ax.pcolorfast(xx, yy, zz.T)
ax.hold(True)
ax.contour(xx, yy, edges.T, colors='k', alpha=0.5)
ax.contour(xx, yy, zz.T, levels=[0], colors='w')
fig.colorbar(hsb)
#ax.set_xlim((-75, -60))
#ax.set_ylim((35, 45))
hsb.set_clim((-6000,0))
ax.set_title('Topo + detected edge')


fig = plt.figure()
ax = fig.gca()
hsb = ax.pcolorfast(xx, yy, edge_topo)
ax.hold(True)
ax.contour(xx, yy, zz.T, levels=[0], colors='w')
fig.colorbar(hsb)
#ax.set_xlim((-75, -60))
#ax.set_ylim((35, 45))
ax.set_title('Topo * detected edge')


Out[195]:
<matplotlib.text.Text at 0x5dcee610>

Slope Burger number