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
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]:
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]:
Gaussian curvature looks noisy and there's a question of extracting the contour from it.
In [172]:
thresh = filter.threshold_otsu(zz)
print(thresh)
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]: