The python-gsw module, used in some of my previous posts, is an incomplete and unstable version for the Matlab gibbs toolbox. Because it is in constant development the best way to obtain it is to go straight to the source code.

The TEOS-10 primer is an excellent place to start before using this toolbox. It supersedes the now obsolete EOS-80 seawater module.

To obtain the latest version of this module you must install mercurial first, then you'll be able to clone the repository:

hg clone https://code.google.com/p/python-gsw/

Then just add the python-gsw/gsw directory into your PYTHONPATH so you can import it from anywhere. To update just go into that directory and type:

hg pull
hg update

My work-flow is to download all modules that are in a development phase into a "pymodules" directory and then create a symlink inside a local folder that is added to the PYTHONPATH in my .bashrc.

Here is a small example of some functions present in the module:

In [2]:

import gsw
import numpy as np
from oceans.sw_extras import gamma_GP_from_SP_pt
from ctd import DataFrame, Series, movingaverage, rosette_summary, lp_filter

def o2sat(s, pt):
    from seawater.library import T68conv
    t = T68conv(pt) + 273.16
    a = (-177.7888, 255.5907, 146.4813, -22.2040)
    b = (-0.037362, 0.016504, -0.0020564)
    lnC = (a[0] + a[1] * (100. / t) + a[2] * np.log(t / 100.) + a[3] *
           (t / 100.) +
           s * (b[0] + b[1] * (t / 100.) + b[2] * (t / 100.) ** 2))
    return np.exp(lnC) * 1000. / 22.392

In [3]:

cast = DataFrame.from_cnv('./data/CTD_001.cnv.gz', compression='gzip')

hours = cast['timeS'].max() / 60. / 60.  # Read metadata.
lon, lat = cast.longitude.mean(), cast.latitude.mean()

# Apply 'bad pump' and 'flag'.
cast = cast[cast['pumps']]  # True for good values.
cast = cast[~cast['flag']]  # True for bad values blame SBE!

keep = set(['t090C', 'c0S/m', 'sbeox0Mm/Kg', 'dz/dtM'])
cast = cast.drop(keep.symmetric_difference(cast.columns), axis=1)

cast['dz/dtM'] = movingaverage(cast['dz/dtM'], window_size=48)
cast['sbeox0Mm/Kg'] = movingaverage(cast['sbeox0Mm/Kg'], window_size=48)

kw = dict(sample_rate=24.0, time_constant=0.15)
cast.index = lp_filter(cast.index, **kw)

downcast, upcast = cast.split()

downcast = downcast.press_check()  # Remove pressure reversals.
downcast = downcast[downcast['dz/dtM'] >= 0.25]  # Threshold velocity.

kw = dict(n1=2, n2=20, block=150)
downcast = downcast.apply(Series.despike, **kw)

downcast = downcast.apply(Series.bindata, **dict(delta=1.))
downcast = downcast.apply(Series.interpolate)

pmax = max(cast.index)
if pmax >= 500:
    window_len = 21
elif pmax >= 100:
    window_len = 11
else:
    window_len = 5
kw = dict(window_len=window_len, window='hanning')
downcast = downcast.apply(Series.smooth, **kw)

In [4]:

p = downcast.index.values.astype(float)
downcast['depth'] = -gsw.z_from_p(p, lat)
# ctm [mS/cm] = ctm [S/m] * 10
downcast['SP'] = gsw.SP_from_C(downcast['c0S/m'].values * 10, downcast['t090C'].values, p)
downcast['SA'] = gsw.SA_from_SP(downcast['SP'].values, p, lon, lat)
downcast['CT'] = gsw.CT_from_t(downcast['SA'].values, downcast['t090C'].values, p)
downcast['sigma0_CT'] = gsw.sigma0_CT_exact(downcast['SA'].values, downcast['CT'].values)
downcast['sound'] = gsw.sound_speed_CT_exact(downcast['SA'].values, downcast['CT'].values, p)
downcast['gamma'] = gamma_GP_from_SP_pt(downcast['SA'].values, downcast['CT'].values, p, lon, lat)
# Does't fit the DataFrame because the index in at the middle (p_mid) point of the original data.
N2, p_mid = gsw.Nsquared(downcast['SA'].values, downcast['CT'].values, p, lat=lat)
downcast['aou'] = o2sat(downcast['SA'].values, downcast['CT'].values) - downcast['sbeox0Mm/Kg']

In [5]:

HTML(html)

Out[5]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

python-gsw

Comments