In [1]:

    

url = ("http://omgsrv1.meas.ncsu.edu:8080/thredds/dodsC/fmrc/sabgom/"
       "SABGOM_Forecast_Model_Run_Collection_best.ncd")


    

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)
     units: 
     long_name: S-coordinate at RHO-points
     valid_min: -1.0
     valid_max: 0.0
     positive: up
     standard_name: ocean_s_coordinate_g1
     formula_terms: s: s_rho C: Cs_r eta: zeta depth: h depth_c: hc
     field: s_rho, scalar
     _CoordinateTransformType: Vertical
     _CoordinateAxisType: GeoZ
     _CoordinateZisPositive: up
     _CoordinateAxes: s_rho
 unlimited dimensions: 
 current shape = (36,)
 filling off, <type 'netCDF4._netCDF4.Variable'>
 float64 s_w(s_w)
     units: 
     long_name: S-coordinate at W-points
     valid_min: -1.0
     valid_max: 0.0
     positive: up
     standard_name: ocean_s_coordinate_g1
     formula_terms: s: s_w C: Cs_w eta: zeta depth: h depth_c: hc
     field: s_w, scalar
     _CoordinateTransformType: Vertical
     _CoordinateAxisType: GeoZ
     _CoordinateZisPositive: up
     _CoordinateAxes: s_w
 unlimited dimensions: 
 current shape = (37,)
 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'C', u'Cs_r'),
             (u'eta', u'zeta'),
             (u'depth', u'h'),
             (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'C', (u's_rho',)),
             (u'eta', (u'time', u'eta_rho', u'xi_rho')),
             (u'depth', (u'eta_rho', u'xi_rho')),
             (u'depth_c', ())])

In [5]:

    

from odvc import z_shape

new_shape = z_shape(nc, dims)

new_shape


    

    
Out[5]:





(496, 36, 320, 440)

In [6]:

    

from odvc import prepare_arrays


arrays = prepare_arrays(nc, formula_terms, new_shape)

arrays


    

    
Out[6]:





{u'C': dask.array<getitem..., shape=(1, 36, 1, 1), dtype=float64, chunksize=(1, 36, 1, 1)>,
 u'depth': dask.array<getitem..., shape=(1, 1, 320, 440), dtype=float64, chunksize=(1, 1, 320, 440)>,
 u'depth_c': array(5.0),
 u'eta': dask.array<getitem..., shape=(496, 1, 320, 440), dtype=float32, chunksize=(1, 1, 320, 440)>,
 u's': dask.array<getitem..., shape=(1, 36, 1, 1), dtype=float64, chunksize=(1, 36, 1, 1)>}

In [7]:

    

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


    

    
Out[7]:





[(u'depth', (1, 1, 320, 440)),
 (u's', (1, 36, 1, 1)),
 (u'eta', (496, 1, 320, 440)),
 (u'C', (1, 36, 1, 1)),
 (u'depth_c', ())]

In [8]:

    

from odvc import ocean_s_coordinate_g1

z = ocean_s_coordinate_g1(arrays['s'],
                          arrays['C'],
                          arrays['eta'],
                          arrays['depth'],
                          arrays['depth_c'])

print(z.shape)


    

    




(496, 36, 320, 440)

In [9]:

    

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, :, 5, 364]
t = temp[-1, :, 5, 364]
p = z[-1, :, 5, 364]


    

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