In [1]:

    

url = ("http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/"
       "Forecasts/NECOFS_GOM3_FORECAST.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'>
 float32 siglay(siglay, node)
     long_name: Sigma Layers
     standard_name: ocean_sigma_coordinate
     positive: up
     valid_min: -1.0
     valid_max: 0.0
     formula_terms: sigma: siglay eta: zeta depth: h
 unlimited dimensions: 
 current shape = (40, 53087)
 filling off]

In [3]:

    

from odvc import get_formula_terms

formula_terms = get_formula_terms(var[0])

formula_terms


    

    
Out[3]:





OrderedDict([(u'sigma', u'siglay'), (u'eta', u'zeta'), (u'depth', u'h')])

In [4]:

    

from odvc import get_formula_terms_dims


dims = get_formula_terms_dims(nc, formula_terms)

dims


    

    
Out[4]:





OrderedDict([(u'sigma', (u'siglay', u'node')),
             (u'eta', (u'time', u'node')),
             (u'depth', (u'node',))])

In [5]:

    

from odvc import z_shape

new_shape = z_shape(nc, dims)

new_shape


    

    
Out[5]:





(145, 40, 53087)

In [6]:

    

from odvc import prepare_arrays

arrays = prepare_arrays(nc, formula_terms, new_shape)

arrays


    

    
Out[6]:





{u'depth': dask.array<getitem..., shape=(1, 1, 53087), dtype=float32, chunksize=(1, 1, 53087)>,
 u'eta': dask.array<getitem..., shape=(145, 1, 53087), dtype=float32, chunksize=(1, 1, 53087)>,
 u'sigma': dask.array<getitem..., shape=(1, 40, 53087), dtype=float32, chunksize=(1, 40, 53087)>}

In [7]:

    

from odvc import ocean_sigma_coordinate


z = ocean_sigma_coordinate(arrays['sigma'],
                           arrays['eta'],
                           arrays['depth'])
print(z.shape)


    

    




(145, 40, 53087)

In [8]:

    

salt = nc.get_variables_by_attributes(standard_name='sea_water_salinity')[0]
temp = nc.get_variables_by_attributes(standard_name='sea_water_potential_temperature')[0]

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


    

In [9]:

    

%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())