physoce package

Submodules

physoce.graph module

physoce.graph.TS_contours(SP_range, T_range, sigma_levels, **kwargs)

Plot contours of density anomaly (sigma) on a T-S plot. Uses EOS-80 equation of state. If the T_range input is in-situ T, then the sigma-t values are contoured. If the temperature input is potential temperature (theta), then the sigma-theta values are contoured (see Stewart 2005 for definitions).

INPUTS: SP_range: Practical salinity, minimum and maximum values T_range: In-situ or potential temperature [C], minimum and maximum values sigma_levels: density anomaly values to contour **kwargs: these will be passed to the contour function, see http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.contour

RETURNS: cs: the matplotlib.contour.QuadContourSet object returned by contour (can

be used as input to pyplot.clabel function, for example).

REQUIRES: Seawater toolbox: https://pypi.python.org/pypi/seawater

Reference: Stewart, R. H. (2005) Introduction to Physical Oceanography. http://oceanworld.tamu.edu/resources/ocng_textbook/chapter06/chapter06_05.htm

physoce.graph.plts(date, y)

Plot a multi-year time series y as a function of time of year (rather than absolute time). The time series from each separate year is plotted in a different color, with months on the x-axis. Useful for visualizing seasonal patterns.

Inputs: date - a list of datetime objects y - a numpy array

physoce.io module

physoce.io.loadmat(filename)

This function is an alternative to directly calling scipy.io.loadmat and cures the problem of not properly recovering python dictionaries from mat files. It calls the function check_keys() to cure all entries which are still mat-objects.

Useful for loading comlicated structures with multiple nested levels.

Source: Stack Overflow http://stackoverflow.com/questions/7008608/scipy-io-loadmat-nested-structures-i-e-dictionaries Author: mergen http://stackoverflow.com/users/887597/mergen License: Creative Commons Attribution-ShareAlike 3.0 http://creativecommons.org/licenses/by-sa/3.0/

physoce.io.print_mat_nested(d, indent=0, nkeys=0)

Print nested structures in a dictionary, d, created with loadmat()

Source: PyHOGS http://pyhogs.github.io/reading-mat-files.html Author: JP Rinehimer

Inspired by: Stack Overflow http://stackoverflow.com/questions/3229419/pretty-printing-nested-dictionaries-in-python Author: sth http://stackoverflow.com/users/56338/sth

Modified by T. Connolly for Python 3 compatibility, as suggested by “Drunken Master” (http://stackoverflow.com/users/4592067/drunken-master) in response to sth’s SO post above

License: Creative Commons Attribution-ShareAlike 3.0 http://creativecommons.org/licenses/by-sa/3.0/

physoce.oceans module

physoce.oceans.ubstream(Hsig, wavper, h, formula='LH', kN=0.03)

Uo = ubstream(Hsig,wavper,h,kN,formula=’LH’)

Calculate bottom streaming velocity (wave-averaged Eulerian velocity at the top of the wave boundary layer).

Inputs: Hsig: significant wave height [m] per: wave period [s] h: bottom depth [m] formula: ‘LH’ - Longuet-Higgins (1953) [default]

‘K’ - Kranenburg (2012)

kN: Nikuradse roughness length (kN = 30*zo)

[default = 0.03m, not needed for Longuet-Higgins formula]

Output: Uo = velocity at top of wave boundary layer [m/s]

References: Longuet-Higgins, M. S. (1953) Mass transport in water waves, Philos. Trans. R.

Soc. London, Ser. A, 245(903), 535-581, doi:10.1098/rsta.1953.0006.

Kranenburg, W. M. et al. (2012) Net currents in the wave bottom boundary layer:

On waveshape streaming and progressive wave streaming, 117, F03005, doi:10.1029/2011JF002070.

physoce.oceans.ubwave(Hsig, wavper, h)

Ub = ubwave(Hsig,wavper,h)

Calculate bottom wave orbital velocity

Inputs: Hsig: significant wave height [m] per: wave period [s] h: bottom depth [m]

Output: Ub: wave orbital velocity near the bottom [m/s]

physoce.oceans.ustokes(Hsig, wavper, h, zst=())

ust,zst = ustokes(Hsig,wavper,h,zst=()): Stokes drift magnitude ——————— Inputs: Hsig - significant wave height [m], single value or 1D array

wavper - wave period [s], single value or 1D array h - bottom depth [m], single value zst - depths of calculation [m] (optional, every 1m if not given)

Outputs: ust - Stokes drift [m/s], 1D or 2D array

zst - depths of calculation [m]

physoce.oceans.ustokesda(Hsig, wavper, h)

Inputs: Hsig - significant wave height [m]

wavper - wave period [s] h - bottom depth [m]

Outputs: depth-averaged Stokes drift [m/s]

physoce.oceans.wavedisp(wavper, h)

Returns [omega,k,Cph,Cg]

Inputs (can use arrays):

wavper - wave period [s] h - water depth [m]

Outputs:

omega - angular wave frequency [radians/s] k - angular wave number [radians/m] Cph - phase speed [m/s] Cg - group velocity [m/s]

physoce.stats module

physoce.stats.maxcorr(x, y, **options)

(rmax,lag,ind) = maxcorr(x,y,**’maxlag’=int(len(x)/4)): Calculate the maximum lagged correlation between two 1D arrays Inputs: x,y are 1D arrays Options ‘maxlag’ the maximum number of lagged correlations to calculate (default: 1/4 of array length) Output: r is the correlation coefficient with the maximum absolute value lag is the lag of the maximum correlation (positive: y lags x)

physoce.stats.nancorr(x, y)

r = nancorr(x,y) Calculate correlation matrix, treating NaN values as missing data

physoce.stats.rcrit(nu, sig=0.05)

Critical r (correlation coefficient), given significance level and degrees of freedom.

INPUTS: nu - degrees of freedom (N-2) sig - significance level (default 0.05)

OUTPUT: rcrit - critical r value

Values for 0.05 and 0.01 correspond with Appendix E in Emery and Thomson (2004) Data Analysis Methods in Physical Oceanography

physoce.stats.rsig(r, nu)

p-value for correlation coefficient r, and degrees of freedom nu

INPUTS: r - correlation coefficient nu - degrees of freedom (N-2)

OUTPUT: p - significance level/p-value

significance levels of 0.05 and 0.01 correspond with Appendix E in Emery and Thomson (2004) Data Analysis Methods in Physical Oceanography

physoce.tseries module

physoce.tseries.depthavg(x, z, h, ssh=None, surface='mixed', bottom='zero')

Compute depth average of each row in 2D array x, with corresponding depths z.

Designed to accomodate upward looking ADCP data, with moving sea surface and blank bins with no data near surface. If no sea surface height is specified, it is assumed to be at z=0 for all times.

x: variable to be depth-averaged, 2D array with shape N rows, M columns z: measurement depths (z=0 is surface, z=-h is bottom), array of length M h: bottom depth (positive value) ssh: sea surface height (optional, set to zero if None or where value is undefined) surface: boundary condition for surface

‘mixed’ (default) or ‘extrap’

bottom: boundary condition for bottom

‘zero’ (default),’mixed’ or ‘extrap’

physoce.tseries.fillgapwithnan(x, date)

Fill in missing data with NaN values. This is intended for a regular time series that has gaps where no data are reported.

Although this function is intended for regularly-sampled time series (with gaps), it does allow for some slight irregularity (e.g. 13:00,14:00,15:00,15:59,16:59…)

Inputs: x - a numpy array of data (1 or 2 dimensions) date - datetime values that correspond to x

Returns: (newx,newdate) newx - new array of data newdate = new datetime values

physoce.tseries.hanning(x, N)

Filter a time series x with a Hanning window of length N. If x is 2D, the time series will be filtered along columns.

Inputs: x - a numpy array to be filtered N - width of window

Output: numpy array of filtered time series

physoce.tseries.lancz(x, dt=1, T=40)

Filter a time series x with cosine-Lanczos filter. If x is 2D, the time series will be filtered along columns.

The default half amplitude period of 40 hours corresponds to a frequency of 0.6 cpd. A half amplitude period of 34.29h corresponds to 0.7 cpd. The 40 hour half amplitude period is more effective at reducing diurnal-band variability but shifts periods of variability in low passed time series to >2 days.

Inputs: x - a numpy array to be filtered. dt - sample interval (hours), default = 1 T - half-amplitude period (hours), default = 40

Output: numpy array of filtered time series, same size as input with ends NaN values at start and end.

Reference: Emery and Thomson, 2004, Data Analysis Methods in Physical Oceanography. 2nd Ed., pp. 539-540. Section 5.10.7.4 - The Hanning window.

physoce.tseries.lombscargle(t, x, ofac=4, hifac=1, t0=None, return_onesided=True, return_zero=False, window='boxcar', scaling='classical')

Compute the discrete Fourier transform and periodogram for unevenly-spaced data using the Lomb-Scargle periodogram. Follows methods outlined in Scargle (1989).

INPUTS

t - array of numerical time values (length N) x - array of data values, may be complex (length N)

RETURNS

f - array of frequencies ftx - array of complex coefficients

discrete Fourier transform of x

px - periodogram of ftx, proportional to |ftx|**2

with the default “classical” scaling used in Scargle (1989), px = (1/N)*|ftx|**2

OPTIONAL PARAMETERS

ofac - oversampling parameter

ratio of number of frequencies used to number of samples in x (default 4)

hifac - high frequency parameter

ratio of highest frequency to pseudo-Nyquist frequency (default 1)

t0 - time origin

reference point for phase calculation (default - None, first value in t array is used)

return_onesided - boolean for returning a one-sided spectrum

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum for real data. Note that for complex data, a two-sided spectrum is always returned. (default - True)

return_zero - boolean for evaluating zero frequency

If True, include zero frequency. If False, do not include zero frequency. Uses expressions for the limit as frequency approaches zero, following Scargle (1989). (default - False)

window - String specifying desired window to use. See scipy.signal.get_window for

a list of windows and required parameters.

scaling - Selects between computing the classical periodogram used by Scargle

(‘classical’) or the power spectral density (‘density’). The classical periodogram has units of x**2. The power spectral density has units of x**2/f. The scaling determines how the periodogram px is calculated from the discrete Fourier transform ftx: ‘classical’: px = (1/N)*|ftx|**2, where N is the number of samples ‘density’: px = (deltat/N)*|ftx|**2, where deltat is the average time step

REFERENCE

Scargle, J.D. (1989) Studies in astronomical time series analysis III: Fourier transforms, autocorrelation functions, and cross-correlation functions of unevenly spaced data. The Astrophysical Journal, 343, 874-887

physoce.tseries.pl64(x, dt=1, T=33)

Filter a time series x with the PL64 filter. If x is 2D, the time series will be filtered along columns.

Inputs: x - a numpy array to be filtered. dt - sample interval (hours), default = 1 T - half-amplitude period (hours), default = 33

Output: numpy array of filtered time series, same size as input with ends NaN values at start and end.

Reference: CODE-2: Moored Array and Large-Scale Data Report, WHOI 85-35

physoce.tseries.pl66(x, dt=1, T=33)

Filter a time series x with the PL66 filter. If x is 2D, the time series will be filtered along columns.

Inputs: x - a numpy array to be filtered. dt - sample interval (hours), default = 1 T - half-amplitude period (hours), default = 33

Output: numpy array of filtered time series, same size as input with ends NaN values at start and end.

Reference: Rosenfeld (1983), WHOI Technical Report 85-35 Matlab code: http://woodshole.er.usgs.gov/operations/sea-mat/bobstuff-html/pl66tn.html

physoce.tseries.princax(u, v=None)

Principal axes of a vector time series.

Usage: theta,major,minor = princax(u,v) # if u and v are real-valued vector components

or

theta,major,minor = princax(w) # if w is a complex vector

Input: u,v - 1-D arrays of vector components (e.g. u = eastward velocity, v = northward velocity)

or

w - 1-D array of complex vectors (u + 1j*v)

Output: theta - angle of major axis (math notation, e.g. east = 0, north = 90) major - standard deviation along major axis minor - standard deviation along minor axis

Reference: Emery and Thomson, 2001, Data Analysis Methods in Physical Oceanography, 2nd ed., pp. 325-328. Matlab function: http://woodshole.er.usgs.gov/operations/sea-mat/RPSstuff-html/princax.html

physoce.tseries.rot(u, v, theta)

Rotate a vector counter-clockwise OR rotate the coordinate system clockwise.

Usage: ur,vr = rot(u,v,theta)

Input: u,v - vector components (e.g. u = eastward velocity, v = northward velocity) theta - rotation angle (degrees)

Output: ur,vr - rotated vector components

Example: rot(1,0,90) returns (0,1)

physoce.tseries.surface_flux(u, z, tr, ssh=None, surface='mixed')

Compute integrated surface-layer flux of a tracer tr, due to velocity component u with corresponding depths z.

Integration is performed from the sea surface to the first zero crossing in the u profile.

If u has units of m/s, and tr has units of C, result will have units of (C m^2)/s

Designed to accomodate upward looking ADCP data, with moving sea surface and blank bins with no data near surface. If no sea surface height is specified, it is assumed to be at z=0 for all times.

INPUT: u: velocity, 2D array with shape N rows, M columns z: measurement depths (z=0 is surface, negative below surface), array of length M tr: tracer values (N rows, M columns corresponding to depths z) ssh: sea surface height (optional, set to zero if None or where value is undefined) surface: boundary condition for surface

‘mixed’ (default) or ‘extrap’

OUTPUT: trflux - tracer flux integrated from surface to first zero crossing of velocity profile tr0 - tracer value at first zero crossing of velocity profile

physoce.tseries.surface_transport(u, z, ssh=None, surface='mixed')

Compute surface transport Us from velocity component u, with corresponding depths z.

Integration is performed from the sea surface to the first zero crossing in the u profile.

If u has units of m/s, Us will have units of m^2/s.

Designed to accomodate upward looking ADCP data, with moving sea surface and blank bins with no data near surface. If no sea surface height is specified, it is assumed to be at z=0 for all times.

INPUT: x: variable to be depth-averaged, 2D array with shape N rows, M columns z: measurement depths (z=0 is surface, negative below surface), array of length M ssh: sea surface height (optional, set to zero if None or where value is undefined) surface: boundary condition for surface

‘mixed’ (default) or ‘extrap’

OUTPUT: Us - transport integrated from surface to first zero crossing of velocity profiled zs - depth of first zero crossing

physoce.util module

General purpose functions.

physoce.util.compasstransform(theta)

Converts angles between compass direction (clockwise from true North) to direction in polar coordinates (counter-clockwise from x-axis, pointing East).

Note that regardless of which way the conversion is being done (compass -> polar, or polar -> compass), the output of this function will be the same for a given input.

INPUT: - theta: direction (degrees), numpy array

OUTPUT: - converted direction (degrees), numpy array of same size

(0 to 360 degrees)

physoce.util.haversine(lat, lon)

Calculate the great circle distances between a set of points on the earth (specified in decimal degrees).

INPUTS: lat,lon - 1D arrays of coordinate pairs (length N)

OUTPUT: great circle distances between consecutive points in km (length N-1)

physoce.util.inside(xpt, ypt, xpoly, ypoly)

Check whether a set of (x,y) points are within a polygon.

INPUTS: xpt,ypt - define a set of points (arrays of length N) xpoly,ypoly - define the vertices of an arbitrary polygon (arrays of length M)

OUTPUT: returns an array of Booleans, length N, True where (xpt,ypt) is inside the polygon

physoce.util.list2date(datestr_list, fmt='%a %b %d %H:%M:%S %Y')

Convert a list of date strings to datetime format.

INPUT: datestr_list: a list of strings that represent dates fmt: format of the date string, as would be input to strftime() or strptime()

see https://docs.python.org/library/datetime.html#strftime-and-strptime-behavior

OUTPUT: list of datetimes

physoce.util.matlab2datetime64(datenum, unit='s')

Convert Matlab serial date number to NumPy datetime64 format.

INPUTS: datenum - Matlab serial date number, can be array unit - time unit of datetime64 output (default ‘s’)

OUTPUT: array of datetime64 objects

Module contents