In [1]:

    

url = ("http://geoport.whoi.edu/thredds/dodsC/examples/bora_feb.nc")


    

In [2]:

    

from netCDF4 import Dataset

from odvc import get_formula_terms_variables


nc = Dataset(url)

var = get_formula_terms_variables(nc)

var


    

    
Out[2]:





[<type 'netCDF4._netCDF4.Variable'>
 float64 s_rho(s_rho)
     long_name: S-coordinate at RHO-points
     valid_min: -1.0
     valid_max: 0.0
     standard_name: ocean_s_coordinate
     formula_terms: s: s_rho eta: zeta depth: h a: theta_s b: theta_b depth_c: hc
     field: s_rho, scalar
 unlimited dimensions: 
 current shape = (20,)
 filling off, <type 'netCDF4._netCDF4.Variable'>
 float64 s_w(s_w)
     long_name: S-coordinate at W-points
     valid_min: -1.0
     valid_max: 0.0
     standard_name: ocean_s_coordinate
     formula_terms: s: s_w eta: zeta depth: h a: theta_s b: theta_b depth_c: hc
     field: s_w, scalar
 unlimited dimensions: 
 current shape = (21,)
 filling off]

In [3]:

    

from odvc import get_formula_terms

formula_terms = get_formula_terms(var[0])

formula_terms


    

    
Out[3]:





OrderedDict([(u's', u's_rho'),
             (u'eta', u'zeta'),
             (u'depth', u'h'),
             (u'a', u'theta_s'),
             (u'b', u'theta_b'),
             (u'depth_c', u'hc')])

In [4]:

    

from odvc import get_formula_terms_dims


dims = get_formula_terms_dims(nc, formula_terms)

dims


    

    
Out[4]:





OrderedDict([(u's', (u's_rho',)),
             (u'eta', (u'ocean_time', u'eta_rho', u'xi_rho')),
             (u'depth', (u'eta_rho', u'xi_rho')),
             (u'a', ()),
             (u'b', ()),
             (u'depth_c', ())])

In [5]:

    

from odvc import z_shape

new_shape = z_shape(nc, dims)

new_shape


    

    
Out[5]:





(8, 20, 60, 160)

In [6]:

    

from odvc import prepare_arrays


arrays = prepare_arrays(nc, formula_terms, new_shape)

arrays


    

    
Out[6]:





{u'a': array(4.0),
 u'b': array(0.9),
 u'depth': dask.array<getitem..., shape=(1, 1, 60, 160), dtype=float64, chunksize=(1, 1, 60, 160)>,
 u'depth_c': array(4.0),
 u'eta': dask.array<getitem..., shape=(8, 1, 60, 160), dtype=float32, chunksize=(1, 1, 60, 160)>,
 u's': dask.array<getitem..., shape=(1, 20, 1, 1), dtype=float64, chunksize=(1, 20, 1, 1)>}

In [7]:

    

[(var, arr.shape) for var, arr in arrays.items()]


    

    
Out[7]:





[(u'a', ()),
 (u's', (1, 20, 1, 1)),
 (u'eta', (8, 1, 60, 160)),
 (u'depth_c', ()),
 (u'b', ()),
 (u'depth', (1, 1, 60, 160))]

In [8]:

    

from odvc import ocean_s_coordinate


z = ocean_s_coordinate(arrays['s'],
                       arrays['eta'],
                       arrays['depth'],
                       arrays['a'],
                       arrays['b'],
                       arrays['depth_c'])

print(z.shape)


    

    




(8, 20, 60, 160)

In [9]:

    

salt = nc.get_variables_by_attributes(long_name='averaged salinity')[0]
temp = nc.get_variables_by_attributes(long_name='averaged potential temperature')[0]

s = salt[-1, :, 45, 137]
t = temp[-1, :, 45, 137]
p = z[-1, :, 45, 137]


    

In [10]:

    

%matplotlib inline

import matplotlib.pyplot as plt
import mpl_toolkits.axisartist as AA
from mpl_toolkits.axes_grid1 import host_subplot

import seaborn
seaborn.set(style='ticks')


def adjust_xlim(ax, offset):
    x1, x2 = ax.get_xlim()
    ax.set_xlim(x1 - offset, x2 + offset)


fig = plt.figure(figsize=(6, 9))
fig.subplots_adjust(bottom=0.1, top=0.85)
ax0 = host_subplot(111, axes_class=AA.Axes)

ax0.invert_yaxis()
ax1 = ax0.twiny()

new_axis0 = ax0.get_grid_helper().new_fixed_axis
new_axis1 = ax1.get_grid_helper().new_fixed_axis

ax0.axis["top"] = new_axis0(loc="top", axes=ax0, offset=(0, 0))
ax1.axis["top"] = new_axis1(loc="top", axes=ax1, offset=(0, 40))

ax0.axis["bottom"].toggle(all=False)
ax1.axis["bottom"].toggle(all=False)

ax0.set_ylabel("Depth (m)")
ax0.set_xlabel(u"Temperature (\xb0C)")
ax1.set_xlabel(r"Salinity (g kg$^{-1}$)")

kw = dict(linewidth=2.5)
l0, = ax0.plot(t, -p, **kw)
l1, = ax1.plot(s, -p, **kw)

adjust_xlim(ax0, offset=0.05)
adjust_xlim(ax1, offset=0.05)

ax0.axis["top"].label.set_color(l0.get_color())
ax1.axis["top"].label.set_color(l1.get_color())