diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index 793723a..20dbed1 100644 --- a/OceanLab/__init.py__ +++ b/OceanLab/__init.py__ @@ -1,6 +1,6 @@ from .dyn import zeta, vmode_amp, psi2uv, eqmodes, vmodes from .eof import my_eof_interp, eoft from .oa import scaloa, vectoa -from utils import argdistnear +from .utils import argdistnear, meaneddy -__version__ = '0.0.7' +_version_ = '0.1.0' diff --git a/OceanLab/dyn.py b/OceanLab/dyn.py index 0b95795..e817cde 100644 --- a/OceanLab/dyn.py +++ b/OceanLab/dyn.py @@ -2,21 +2,36 @@ import numpy as np import seawater as sw +# =========================================================== +# DYNAMICAL MODES AMPLITUDE +# =========================================================== def vmode_amp(A,vi): ''' This function makes the projection of every vertical mode to - timeseries or section matrix to obtain its amplitude. It will be used to - data reconstruction. + timeseries or section matrix to obtain its amplitude. + It will be used to data reconstruction. Operation: ((A'*A)^-1)*(A'*vi) + ============================================================= + + INPUT: + A = + vi = + OUTPUT: + AMP = + + ============================================================== ''' - return np.dot(np.linalg.inv(np.dot(A.T,A)),np.dot(A.T,vi)) + return np.dot(np.linalg.inv(np.dot(A.T,A)),np.dot(A.T,vi)) +# =========================================================== +# RELATIVE VORTICITY +# =========================================================== def zeta(x,y,U,V): ''' - ZETA k-component of rotational by velocity field + ZETA k-component of rotational by velocity field Returns the scalar field of vorticity from U and V velocity components. X and Y are longitude and latitude matrices, respectively, as U and V. All matrices have @@ -39,25 +54,24 @@ def zeta(x,y,U,V): (JM = 0) o ---------- o ---------- o ---------- o --- ... (IM = 0) (IM = 1) (IM = 2) (IM = 3) - ----------------------------------------------------------------- + ================================================================ + USAGE: ZETA = zeta(x,y,u,v) - USAGE: u,v = psi2uv(x,y,psi) + INPUT: + x = decimal degrees (+ve E, -ve W) [-180..+180] + y = decimal degrees (+ve N, -ve S) [- 90.. +90] + u = velocity zonal component [m s^-1] + v = velocity meridional component [m s^-1] - INPUT: - x = decimal degrees (+ve E, -ve W) [-180..+180] - y = decimal degrees (+ve N, -ve S) [- 90.. +90] - u = velocity zonal component [m s^-1] - v = velocity meridional component [m s^-1] + OUTPUT: + ZETA = Relative vorticity field [s^-1] - OUTPUT: - ZETA = Relative vorticity field [s^-1] - - AUTHOR: + AUTHOR: Iury T. Simoes-Sousa and Wandrey Watanabe - 29 Jun 2016 Laboratório de Dinâmica Oceânica - IOUSP - ====================================================================== - + ================================================================ ''' + #função para cálculo de ângulos angcalc = lambda dy,dx: np.math.atan2(dy,dx) @@ -116,16 +130,21 @@ def zeta(x,y,U,V): return ZETA +# =========================================================== +# U AND V FROM STREAMFUNCTION +# =========================================================== def psi2uv(x,y,psi): ''' - PSI2UV Velocity components from streamfunction - Returns the velocity components U and V - from streamfunction PSI. X and Y are longitude and latitude matrices, - respectively, as PSI. All matrices have same dimension. + PSI2UV Velocity components from streamfunction - --> !!! IMPORTANT !!! <----------------------------------------- - The grid indexes IM and JM must have origin on the left-inferior - corner, increasing to right and to up, respectively. + Returns the velocity components U and V + from streamfunction PSI. X and Y are longitude and latitude + matrices, respectively, as PSI. + All matrices have same dimension. + + --> !!! IMPORTANT !!! <----------------------------------- + The grid indexes IM and JM must have origin on the + lower-left corner, increasing to right and to up. GRID : : : : @@ -139,23 +158,24 @@ def psi2uv(x,y,psi): (JM = 0) o ---------- o ---------- o ---------- o --- ... (IM = 0) (IM = 1) (IM = 2) (IM = 3) - ----------------------------------------------------------------- + ============================================================== - USAGE: u,v = psi2uv(x,y,psi) + USAGE: u,v = psi2uv(x,y,psi) - INPUT: - x = decimal degrees (+ve E, -ve W) [-180..+180] - y = decimal degrees (+ve N, -ve S) [- 90.. +90] - psi = streamfunction [m^2 s^-1] + INPUT: + x = decimal degrees (+ve E, -ve W) [-180..+180] + y = decimal degrees (+ve N, -ve S) [- 90.. +90] + psi = streamfunction [m^2 s^-1] - OUTPUT: - u = velocity zonal component [m s^-1] - v = velocity meridional component [m s^-1] + OUTPUT: + u = velocity zonal component [m s^-1] + v = velocity meridional component [m s^-1] - AUTHOR: - Iury T. Simoes-Sousa and Wandrey Watanabe - 12 May 2016 - Laboratório de Dinâmica Oceânica - IOUSP - ====================================================================== + AUTHOR: + Iury T. Simoes-Sousa and Wandrey Watanabe + May 2016 - Laboratório de Dinâmica Oceânica - IOUSP + Universidade de São Paulo + ========================================================== ''' #função para cálculo de ângulos @@ -210,45 +230,38 @@ def psi2uv(x,y,psi): return U,V - +# =========================================================== +# EQUATORIAL MODES +# =========================================================== def eqmodes(N2,z,nm,lat,pmodes=False): ''' This function computes the equatorial velocity modes - ======================================================== - - Input: - - N2 - Brunt-Vaisala frequency data array - - z - Depth data array (equaly spaced) - - lat - Latitude scalar - - nm - Number of modes to be computed - - pmodes - If the return of pressure modes is required. - Default is False - - ======================================================== - - Output: - - Si - Equatorial modes matrix with MxN dimension being: - M = z array size - N = nm - - Rdi - Deformation Radii array - - Fi - Pressure modes matrix with MxN dimension being: - M = z array size - N = nm - - Returned only if input pmodes=True - - made by Iury T. Simoes-Sousa, Hélio Almeida and Wandrey Watanabe - Laboratório de Dinâmica Oceânica - Universidade de São Paulo - 2016 + ============================================================ + + INPUT: + N2 = Brunt-Vaisala frequency data array + z = Depth data array (equaly spaced) + lat = Latitude scalar + nm = Number of modes to be computed + pmodes = If the return of pressure modes is required. + (Default is False) + + OUTPUT: + Si = Equatorial modes matrix with MxN dimension being: + M = z array size + N = nm + Rdi = Deformation Radii array + Fi = Pressure modes matrix with MxN dimension being: + M = z array size + N = nm + (Returned only if input pmodes=True) + + AUTHOR: + Iury T. Simoes-Sousa, Hélio Almeida and Wandrey Watanabe + 2016 Laboratório de Dinâmica Oceânica + Universidade de São Paulo + ============================================================ ''' #nm will be the number of baroclinic modes @@ -352,8 +365,22 @@ def eqmodes(N2,z,nm,lat,pmodes=False): return Si,radii - +# =========================================================== +# DYNAMICAL MODES +# =========================================================== def vmodes(N2,z,nm,lat,ubdy='N',lbdy='N'): + + """ + + ======================================================== + INPUT: + + OUTPUT: + + AUTHOR: + ======================================================== + """ + f0 = sw.f(lat) dz = np.abs(z[1] - z[0]) @@ -417,4 +444,4 @@ def vmodes(N2,z,nm,lat,ubdy='N',lbdy='N'): radii[0] = np.sqrt(9.81*H)/np.abs(f0)/1000 radii[1:] = 1/np.sqrt(ei[1:])/1000 - return fi, radii \ No newline at end of file + return fi, radii diff --git a/OceanLab/utils.py b/OceanLab/utils.py index 2044ce3..8e793a7 100644 --- a/OceanLab/utils.py +++ b/OceanLab/utils.py @@ -1,23 +1,162 @@ import numpy as np +import scipy.signal as sg +import xarray as xr -##### Function +import dask +from dask.distributed import Client, progress + +##### User functions #============================================================================= # NEAREST DISTANCE #============================================================================= def argdistnear(x,y,xi,yi): ''' - This function finds the index to nearest points in (xi,yi) from (x,y). - - usage: + This function finds the index to nearest points in (xi,yi) from (x,y) + ====================================================================== + + USAGE: x,y = [5,1,10],[2,6,3] xi,yi = np.linspace(0,19,20),np.linspace(-5,30,20) ind = argdistnear(x,y,xi,yi) INPUT: - --> (x,y): points [list] - --> (xi,yi): series to search nearest point [list] + (x,y) = points [list] + (xi,yi) = series to search nearest point [list] + + OUTPUT: + ind = index of the nearest points + ====================================================================== ''' + idxs = [np.argmin(np.sqrt((xi-xx)**2 + (yi-yy)**2)) for xx,yy in zip(x,y)] idxs = np.array(idxs) return idxs #============================================================================= + +#============================================================================= +# LOW PASS FILTER +#============================================================================= +def meaneddy(prop,days=60,ndim=1,DataArray=False,timedim=None): + + """ + Apply a low-pass filter (scipy.signal.butter) to 'prop' and obtain the + mean and eddy components. + ========================================================================== + + USAGE [1]: + # np.array() one year, 17 lat x 13 lon domain + Velocity = np.random.randn(365,17,13) + Filtered, Residual = meaneddy(Velocity, days=10, ndim=3, + DataArray=False,timedim=None) + + USAGE [2]: + # xr.DataArray() one year, 17 lat x 13 lon domain + Velocity = xr.DataArray(data=np.random.randn(365,17,13), + dims=["time","lat","lon"], + coords=dict(time=(["time"],range(0,365)), + lat=(["lat"],np.arange(-4,4.5,0.5)), + lon=(["lon"],np.arange(1,7.5,0.5)))) + Filtered, Residual = meaneddy(Velocity, days=10, + DataArray=True,timedim=["time"]) + + INPUT: + prop = 1, 2 or 3D array to be filtered + days = number of days to set up the filter + ndim = number of dimensions of the data + [only used for DataArray=False, max:3] + DataArray = True if prop is in xr.DataArray format + dim = name of time dimension to filter + [only used for DataArray=True] + + OUTPUT: + m_prop = mean (filtered) part of the property + p_prop = prime part of the property, essentially prop - m_prop + ========================================================================== + """ + + # creating filter + def timefilter(prop,filtdays=60): + filt_b,filt_a = sg.butter(4,1./filtdays) + return sg.filtfilt(filt_b,filt_a,prop) + + if DataArray: + m_prop = xr.apply_ufunc( + timefilter, # first the function + prop,# now arguments in the order expected by 'butter_filt' + input_core_dims=[timedim], # list with one entry per arg + output_core_dims=[timedim], # returned data + kwargs={'filtdays':days}, + vectorize=True, # loop over non-core dims + dask='vectorized') + + p_prop = prop - m_prop + + elif ndim ==3: + ti,lt,ln = prop.shape + prop = prop.reshape(ti,lt*ln) + + m_prop = [] + + for tot in prop.T: + # filtered series (mean) + m_prop.append(timefilter(tot,filtdays=days)) + + m_prop = np.array(m_prop) + m_prop = m_prop.T + + p_prop = prop - m_prop + + m_prop = m_prop.reshape(ti,lt,ln) + p_prop = p_prop.reshape(ti,lt,ln) + + elif ndim ==2: + + ti,lt = prop.shape + + m_prop = [] + + for tot in prop.T: + # filtered series (mean) + m_prop.append(timefilter(tot,filtdays=days)) + + m_prop = np.array(m_prop) + m_prop = m_prop.T + + p_prop = prop - m_prop + + m_prop = m_prop.reshape(ti,lt) + p_prop = p_prop.reshape(ti,lt) + + elif ndim ==1: + m_prop = timefilter(prop,filtdays=days) + p_prop = prop - m_prop + + return m_prop,p_prop +#============================================================================= + +##### Functions for relative imports +# ============================================================================= +# KERNEL FOR PARALLEL COMPUTING +# ============================================================================= +def _parallel_client(cpu_params=dict(tpw=2,nw=4,ml=7.5)): + """ + Create client kernel for parallel computing + ==================================================== + INPUT: + -> cpu_params: dict containing floats with keys + -> tpw: threads_per_worker + -> nw: n_workers + -> ml: memory_limit per worker [GB] + OUTPUT: + -> client: configuration of parallel computing + ==================================================== + """ + + client = Client(threads_per_worker=cpu_params['tpw'], + n_workers=cpu_params['nw'], + memory_limit=str(cpu_params['ml'])+'GB') + return client +#============================================================================= + + + diff --git a/requirements.txt b/requirements.txt index 9d87074..948d657 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,6 @@ numpy >= 1.8.2 -seawater >= 3.3.1 \ No newline at end of file +seawater >= 3.3.1 +scipy >= 1.6.3 +xarray >= 0.18.2 +dask >= 2021.06.0 +dask[distributed] >= 2021.06.0 diff --git a/setup.py b/setup.py index 5fae7c5..8a31867 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ def read(file): setup( name='OceanLab', - version='0.0.7', + version='0.1.0', packages=['OceanLab'], include_package_data=True, description='Python functions for Physical Oceanography', @@ -30,5 +30,9 @@ def read(file): install_requires = [ 'seawater ~= 3.3', 'numpy ~= 1.18', + 'scipy ~= 1.6', + 'xarray ~= 0.18', + 'dask ~= 2021.06', + 'dask[distributed] ~= 2021.06' ], -) \ No newline at end of file +)