From 45210ca833f17c5ed708bfb1db5899dca030f2fb Mon Sep 17 00:00:00 2001 From: iury simoes-sousa Date: Tue, 8 Sep 2020 18:00:10 -0400 Subject: [PATCH 01/40] Manifest created --- MANIFEST.in | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 MANIFEST.in diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..b88fbb8 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,3 @@ +include *.txt +recursive-include examples *.ipynb +recursive-include src *.py From a9f006e0dcf19260d639f4c76da84e7e5e24ef67 Mon Sep 17 00:00:00 2001 From: iury simoes-sousa Date: Tue, 8 Sep 2020 19:59:48 -0400 Subject: [PATCH 02/40] v0.0.5 --- MANIFEST.in | 3 --- README.md | 14 ++------------ setup.py | 4 ++-- 3 files changed, 4 insertions(+), 17 deletions(-) delete mode 100644 MANIFEST.in diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index b88fbb8..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,3 +0,0 @@ -include *.txt -recursive-include examples *.ipynb -recursive-include src *.py diff --git a/README.md b/README.md index 0d68788..6090204 100644 --- a/README.md +++ b/README.md @@ -2,25 +2,15 @@ Package of Python scripts for Oceanography (Python +3.6) -## Synopsis - -This lib was made by **Iury T. Simoes-Sousa** (). - -Now this package have this subpackages below: - -- **OA** -- **DYN** -- **EOF** - ## Code Example -Check `examples` folder in our [github repository](github.com/iuryt/OceanLab). +Check `examples` folder in our [github repository](github.com/iuryt/OceanLab/examples). ## Installation `pip install OceanLab` -## Documentation +## Modules - **OA** - *vectoa()*: Objective analysis for vectorial fields; diff --git a/setup.py b/setup.py index de260cc..5472d5d 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,7 @@ setup( name='OceanLab', - version='0.0.1', + version='0.0.5', description='Python functions for Physical Oceanography', long_description=long_description, long_description_content_type='text/markdown', @@ -15,7 +15,7 @@ author_email='simoesiury@gmail.com', license='MIT', py_modules=['OA','EOF','DYN'], - package_dir={'': 'src'}, + package_dir={'': 'src/'}, classifiers=[ "Programming Language :: Python :: 3", "Programming Language :: Python :: 3.6", From da71fe4fdf206179b8d1d884f9d6977c8e4cf7be Mon Sep 17 00:00:00 2001 From: iury simoes-sousa Date: Tue, 8 Sep 2020 20:00:33 -0400 Subject: [PATCH 03/40] Add Manifest.in --- MANIFEST.in | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 MANIFEST.in diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..b88fbb8 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,3 @@ +include *.txt +recursive-include examples *.ipynb +recursive-include src *.py From b2fa21a532fbc4ee074d7a6ff500ad1387c1ed88 Mon Sep 17 00:00:00 2001 From: "Helio Almeida (CTW)" Date: Wed, 9 Sep 2020 20:00:21 +0100 Subject: [PATCH 04/40] Add init --- {src/OceanLab => OceanLab}/DYN.py | 0 {src/OceanLab => OceanLab}/EOF.py | 0 {src/OceanLab => OceanLab}/OA.py | 0 OceanLab/__init.py__ | 3 +++ setup.py | 1 + 5 files changed, 4 insertions(+) rename {src/OceanLab => OceanLab}/DYN.py (100%) rename {src/OceanLab => OceanLab}/EOF.py (100%) rename {src/OceanLab => OceanLab}/OA.py (100%) create mode 100644 OceanLab/__init.py__ diff --git a/src/OceanLab/DYN.py b/OceanLab/DYN.py similarity index 100% rename from src/OceanLab/DYN.py rename to OceanLab/DYN.py diff --git a/src/OceanLab/EOF.py b/OceanLab/EOF.py similarity index 100% rename from src/OceanLab/EOF.py rename to OceanLab/EOF.py diff --git a/src/OceanLab/OA.py b/OceanLab/OA.py similarity index 100% rename from src/OceanLab/OA.py rename to OceanLab/OA.py diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ new file mode 100644 index 0000000..8e1d83a --- /dev/null +++ b/OceanLab/__init.py__ @@ -0,0 +1,3 @@ +from DYN import * +from EOF import * +from OA import * \ No newline at end of file diff --git a/setup.py b/setup.py index 5472d5d..da46085 100644 --- a/setup.py +++ b/setup.py @@ -7,6 +7,7 @@ setup( name='OceanLab', version='0.0.5', + packages = ['OceanLab'], description='Python functions for Physical Oceanography', long_description=long_description, long_description_content_type='text/markdown', From fd9f19c2180661d2bac2d1b65c7a52ed5dbf1911 Mon Sep 17 00:00:00 2001 From: "Helio Almeida (CTW)" Date: Wed, 9 Sep 2020 20:37:46 +0100 Subject: [PATCH 05/40] Fix init --- OceanLab/__init.py__ | 8 +++++--- OceanLab/{DYN.py => dyn.py} | 0 OceanLab/{EOF.py => eof.py} | 0 OceanLab/{OA.py => oa.py} | 0 setup.py | 19 +++++++++++-------- 5 files changed, 16 insertions(+), 11 deletions(-) rename OceanLab/{DYN.py => dyn.py} (100%) rename OceanLab/{EOF.py => eof.py} (100%) rename OceanLab/{OA.py => oa.py} (100%) diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index 8e1d83a..43aaeed 100644 --- a/OceanLab/__init.py__ +++ b/OceanLab/__init.py__ @@ -1,3 +1,5 @@ -from DYN import * -from EOF import * -from OA import * \ No newline at end of file +from .dyn import zeta, vmode_amp, psi2uv, eqmodes, vmodes +from .eof import my_eof_interp, eoft +from .oa import scaloa, vectoa + +__version__ = '0.0.7' \ No newline at end of file diff --git a/OceanLab/DYN.py b/OceanLab/dyn.py similarity index 100% rename from OceanLab/DYN.py rename to OceanLab/dyn.py diff --git a/OceanLab/EOF.py b/OceanLab/eof.py similarity index 100% rename from OceanLab/EOF.py rename to OceanLab/eof.py diff --git a/OceanLab/OA.py b/OceanLab/oa.py similarity index 100% rename from OceanLab/OA.py rename to OceanLab/oa.py diff --git a/setup.py b/setup.py index da46085..5fae7c5 100644 --- a/setup.py +++ b/setup.py @@ -1,22 +1,25 @@ from setuptools import setup +def read(file): + return open(file, 'r').read() -with open('README.md','r') as fh: - long_description = fh.read() +LONG_DESCRIPTION = read('README.md') +LICENSE = read('LICENSE.txt') setup( name='OceanLab', - version='0.0.5', - packages = ['OceanLab'], + version='0.0.7', + packages=['OceanLab'], + include_package_data=True, description='Python functions for Physical Oceanography', - long_description=long_description, + long_description=LONG_DESCRIPTION, long_description_content_type='text/markdown', + download_url = 'https://pypi.python.org/pypi/OceanLab', url='https://github.com/iuryt/OceanLab', author='Iury T. Simoes-Sousa', author_email='simoesiury@gmail.com', - license='MIT', + license=LICENSE, py_modules=['OA','EOF','DYN'], - package_dir={'': 'src/'}, classifiers=[ "Programming Language :: Python :: 3", "Programming Language :: Python :: 3.6", @@ -28,4 +31,4 @@ 'seawater ~= 3.3', 'numpy ~= 1.18', ], -) +) \ No newline at end of file From 7052f5c8eff37334cb68735315156274beba4aeb Mon Sep 17 00:00:00 2001 From: dantecn Date: Mon, 10 May 2021 19:08:36 +0200 Subject: [PATCH 06/40] add function description --- OceanLab/eof.py | 54 ++++++++++++++++++++++++------------------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index 60a0244..5d8338c 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -1,5 +1,9 @@ import numpy as np +# functions +#========================================= +# COMPUTE EOFS +#========================================= def eoft(trmat,nm=None): ''' evals_perc,evecs_norm,amp = eoft(trmat) @@ -42,35 +46,31 @@ def eoft(trmat,nm=None): return evals_perc,evecs_norm,amp +#========================================= +# FILL GAPS WITH EOFS +#========================================= +def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): + """ + vi = my_eof_interp(M,nmodes) -def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): - #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vi = my_eof_interp(M,nmodes) - #% This function try to fill gappy data using its EOFs - #% - #% M is the data Matrix (generaly the diferent stations - #% are on the columns); - #% - #% nmodes is the number of statiscally significant EOF modes (to be used for the series reconstruction) - #% - #%% The method -- based on Beckers & Rixen (2003, JAOTech) - #% - #% The algorithm proposed by the authors is as follows. (see Eqs. (13) and (14) from Beckers & Rixen) - #% - #% (I) Put zero in the GAPPY positions ("unbiased INITIAL GUESS, by assuming that the removed mean was unbiased). - #% - #% (II) Then you can have a first estimative of the EOFs. - #% - #% (III) You can reconstruct the series using the statiscally significant modes (e.g.: determined using Monte Carlo - #% simulations as in Preisendorfer [1988] -- OBS: the authors propose another method using cross validation - #% during the iterative process) - #% - #% (IV) Substitute the reconstructed series ONLY in the GAPPY positions in your original matrix. - #% - #% (V) Run steps (II), (III) and (IV) iteractively till its convergence (i.e. compare the reconstructed series at time - #% t+1 with that at time t). - #% - #% OBS: It doesn't work for filling gaps simultaneously in all depths. + This function try to fill gappy data using its EOFs + + INPUT: + -> M: the data Matrix (generaly the diferent stations are on the columns); + -> nmodes: the number of statiscally significant EOF modes (to be used for the series reconstruction) + + The method -- based on Beckers & Rixen (2003, JAOTech) + + The algorithm proposed by the authors is as follows. (see Eqs. (13) and (14) from Beckers & Rixen) + (I) Put zero in the GAPPY positions ("unbiased INITIAL GUESS, by assuming that the removed mean was unbiased). + (II) Then you can have a first estimative of the EOFs. + (III) You can reconstruct the series using the statiscally significant modes (e.g.: determined using Monte Carlo simulations as in Preisendorfer [1988] -- OBS: the authors propose another method using cross validation during the iterative process) + (IV) Substitute the reconstructed series ONLY in the GAPPY positions in your original matrix. + (V) Run steps (II), (III) and (IV) iteractively till its convergence (i.e. compare the reconstructed series at time t+1 with that at time t). + + %%%%% OBS: my_eof_interp doesn't work for filling gaps simultaneously in all depths. + """ Mmean = np.nanmean(M,axis=1) From 8b739a2d2445bfa6851fcf824aaccbd59e19dc5d Mon Sep 17 00:00:00 2001 From: dantecn Date: Mon, 10 May 2021 19:09:33 +0200 Subject: [PATCH 07/40] add utils module with argdistnear function --- OceanLab/__init.py__ | 3 ++- OceanLab/utils.py | 23 +++++++++++++++++++++++ 2 files changed, 25 insertions(+), 1 deletion(-) create mode 100644 OceanLab/utils.py diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index 43aaeed..793723a 100644 --- a/OceanLab/__init.py__ +++ b/OceanLab/__init.py__ @@ -1,5 +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 -__version__ = '0.0.7' \ No newline at end of file +__version__ = '0.0.7' diff --git a/OceanLab/utils.py b/OceanLab/utils.py new file mode 100644 index 0000000..2044ce3 --- /dev/null +++ b/OceanLab/utils.py @@ -0,0 +1,23 @@ +import numpy as np + +##### Function +#============================================================================= +# NEAREST DISTANCE +#============================================================================= +def argdistnear(x,y,xi,yi): + ''' + 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] + ''' + idxs = [np.argmin(np.sqrt((xi-xx)**2 + (yi-yy)**2)) for xx,yy in zip(x,y)] + idxs = np.array(idxs) + return idxs +#============================================================================= From 7d41d38f9d1fc088d588fd513430a1d0a999cebd Mon Sep 17 00:00:00 2001 From: dantecn Date: Wed, 2 Jun 2021 16:22:54 +0200 Subject: [PATCH 08/40] added mean-eddy decomposition function to utils --- OceanLab/__init.py__ | 2 +- OceanLab/utils.py | 104 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 105 insertions(+), 1 deletion(-) diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index 793723a..0eb8a07 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' diff --git a/OceanLab/utils.py b/OceanLab/utils.py index 2044ce3..836cee5 100644 --- a/OceanLab/utils.py +++ b/OceanLab/utils.py @@ -1,4 +1,6 @@ import numpy as np +import scipy.signal as sg +import xarray as xr ##### Function #============================================================================= @@ -16,8 +18,110 @@ def argdistnear(x,y,xi,yi): INPUT: --> (x,y): points [list] --> (xi,yi): series to search nearest point [list] + + Iury T.Simões-Sousa + (IO-USP/ UMass-Dartmouth) ''' + 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]: + Velocity = np.random.randn(365,17,13) # one year, 17 lat x 13 lon domain + Filtered, Residual = meaneddy(Velocity, days=10, ndim=3, DataArray=False,timedim=None) + + usage [2]: + 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)))) # one year, 17 lat x 13 lon domain + Filtered, Residual = meaneddy(Velocity, days=10, DataArray=True,timedim=["time"]) + + + INPUT: + -> prop: 1, 2 or 3D array to filter + -> 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 + + v1 (February 2018) + Cesar B. Rocha + Dante C. Napolitano (dante.napolitano@legos.obs-mip.fr) + + v2 (December 2020) + Dante C. Napolitano (dante.napolitano@legos.obs-mip.fr) + Mariana M. Lage (mariana.lage@hereon.de) + """ + + # 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 +#============================================================================= From 1e2030524926f4ec61ede324d113aca81d5d0780 Mon Sep 17 00:00:00 2001 From: Dante Campagnoli Napolitano Date: Thu, 3 Jun 2021 14:26:32 +0200 Subject: [PATCH 09/40] Update README.md --- README.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 6090204..663bc46 100644 --- a/README.md +++ b/README.md @@ -24,7 +24,11 @@ Check `examples` folder in our [github repository](github.com/iuryt/OceanLab/exa - **EOF** - *eoft()*: Calculates the Empirical Orthogonal Functions; - *my_eof_interp()*: Fillgaps on matrix based on EOFs (translated from Cesar Rocha Matlab version); +- **UTILS** + - *argdistnear()*: Searchs the position of the closest points in an array to a reference point; + - *meaneddy()*: Performs an eddy-mean decomposing with a low-pass filter; + ## Contributors -Everyone can contribute to this code. Some of functions was based on Filipe Fernandes or Cesar Rocha functions and some of them was created with help of Dante C. Napolitano, Hélio M. R. Almeida and Wandrey Watanabe at Ocean Dynamics Lab of University of São Paulo (USP). \ No newline at end of file +Everyone can contribute to this code. Some of functions were based on Filipe Fernandes or Cesar Rocha functions and some of them were created with help of Dante C. Napolitano, Hélio M. R. Almeida and Wandrey Watanabe at Ocean Dynamics Lab of University of São Paulo (USP). From 4e72774dbad9cf56ebcfacb60c587881a0af397d Mon Sep 17 00:00:00 2001 From: dantecn Date: Sat, 5 Jun 2021 15:10:00 +0200 Subject: [PATCH 10/40] added parallel_client (dask) function to utils --- OceanLab/__init.py__ | 2 +- OceanLab/utils.py | 27 ++++++++++++++++++++++++++- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index 0eb8a07..fe3497a 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, meaneddy +from utils import parallel_client, argdistnear, meaneddy __version__ = '0.0.7' diff --git a/OceanLab/utils.py b/OceanLab/utils.py index 836cee5..bfd1e2d 100644 --- a/OceanLab/utils.py +++ b/OceanLab/utils.py @@ -2,7 +2,32 @@ import scipy.signal as sg import xarray as xr -##### Function +import dask +from dask.distributed import Client, progress + +##### Functions +# ============================================================================= +# 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 +#============================================================================= + #============================================================================= # NEAREST DISTANCE #============================================================================= From 69948351f0b318218c2bd8bb81c2ea0a7e22a1c9 Mon Sep 17 00:00:00 2001 From: dantecn Date: Sat, 5 Jun 2021 16:10:30 +0200 Subject: [PATCH 11/40] updated help information on some of dyn.py functions --- OceanLab/dyn.py | 175 ++++++++++++++++++++++++++++-------------------- 1 file changed, 101 insertions(+), 74 deletions(-) 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 From 2ebeac5cb92cdc12bb0b4ff05c3a192548c28adb Mon Sep 17 00:00:00 2001 From: dantecn Date: Mon, 7 Jun 2021 13:20:55 +0200 Subject: [PATCH 12/40] remove parallel_client from __init__. Can be only imported between modules --- OceanLab/__init.py__ | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index fe3497a..0eb8a07 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 parallel_client, argdistnear, meaneddy +from utils import argdistnear, meaneddy __version__ = '0.0.7' From 1bf69ae6befb31ecaabfe73c5973fdc92f59cb76 Mon Sep 17 00:00:00 2001 From: dantecn Date: Tue, 8 Jun 2021 09:25:30 +0200 Subject: [PATCH 13/40] added version 0.1.0 utils scripts --- OceanLab/__init.py__ | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index 0eb8a07..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, meaneddy +from .utils import argdistnear, meaneddy -__version__ = '0.0.7' +_version_ = '0.1.0' From dcfc78065b1940af0cb99fc18877a7422214284d Mon Sep 17 00:00:00 2001 From: dantecn Date: Tue, 8 Jun 2021 09:33:49 +0200 Subject: [PATCH 14/40] updated functions documentation --- OceanLab/utils.py | 74 +++++++++++++++++++++++++---------------------- 1 file changed, 40 insertions(+), 34 deletions(-) diff --git a/OceanLab/utils.py b/OceanLab/utils.py index bfd1e2d..adb81f0 100644 --- a/OceanLab/utils.py +++ b/OceanLab/utils.py @@ -12,7 +12,7 @@ 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 @@ -20,6 +20,7 @@ def parallel_client(cpu_params=dict(tpw=2,nw=4,ml=7.5)): -> ml: memory_limit per worker [GB] OUTPUT: -> client: configuration of parallel computing + ==================================================== """ client = Client(threads_per_worker=cpu_params['tpw'], @@ -33,19 +34,21 @@ def parallel_client(cpu_params=dict(tpw=2,nw=4,ml=7.5)): #============================================================================= 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] - - Iury T.Simões-Sousa - (IO-USP/ UMass-Dartmouth) + (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)] @@ -59,36 +62,39 @@ def argdistnear(x,y,xi,yi): 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]: - Velocity = np.random.randn(365,17,13) # one year, 17 lat x 13 lon domain - Filtered, Residual = meaneddy(Velocity, days=10, ndim=3, 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]: - 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)))) # one year, 17 lat x 13 lon domain - Filtered, Residual = meaneddy(Velocity, days=10, DataArray=True,timedim=["time"]) - + 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 filter - -> 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) + 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 - - v1 (February 2018) - Cesar B. Rocha - Dante C. Napolitano (dante.napolitano@legos.obs-mip.fr) - - v2 (December 2020) - Dante C. Napolitano (dante.napolitano@legos.obs-mip.fr) - Mariana M. Lage (mariana.lage@hereon.de) + m_prop = mean (filtered) part of the property + p_prop = prime part of the property, essentially prop - m_prop + ========================================================================== """ # creating filter From 0d0de7114c9a86b4ddbf1ab7f2e4e66d38dfec89 Mon Sep 17 00:00:00 2001 From: dantecn Date: Tue, 8 Jun 2021 09:34:51 +0200 Subject: [PATCH 15/40] updated requirements for current version --- requirements.txt | 6 +++++- setup.py | 45 +++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 48 insertions(+), 3 deletions(-) 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..8ae11e2 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,46 @@ 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 +) + +# OLDER VERSION (Before MeanEddy branch Jun 8 2021) + +# from setuptools import setup + +# def read(file): +# return open(file, 'r').read() + +# LONG_DESCRIPTION = read('README.md') +# LICENSE = read('LICENSE.txt') + +# setup( +# name='OceanLab', +# version='0.0.7', +# packages=['OceanLab'], +# include_package_data=True, +# description='Python functions for Physical Oceanography', +# long_description=LONG_DESCRIPTION, +# long_description_content_type='text/markdown', +# download_url = 'https://pypi.python.org/pypi/OceanLab', +# url='https://github.com/iuryt/OceanLab', +# author='Iury T. Simoes-Sousa', +# author_email='simoesiury@gmail.com', +# license=LICENSE, +# py_modules=['OA','EOF','DYN'], +# classifiers=[ +# "Programming Language :: Python :: 3", +# "Programming Language :: Python :: 3.6", +# "Programming Language :: Python :: 3.7", +# "License :: OSI Approved :: MIT License", +# "Operating System :: OS Independent", +# ], +# install_requires = [ +# 'seawater ~= 3.3', +# 'numpy ~= 1.18', +# ], +# ) \ No newline at end of file From b4a537c5205078cfced07efba532f471a503bed3 Mon Sep 17 00:00:00 2001 From: dantecn Date: Wed, 9 Jun 2021 12:32:26 +0200 Subject: [PATCH 16/40] removed coomented old version --- setup.py | 37 ------------------------------------- 1 file changed, 37 deletions(-) diff --git a/setup.py b/setup.py index 8ae11e2..8a31867 100644 --- a/setup.py +++ b/setup.py @@ -36,40 +36,3 @@ def read(file): 'dask[distributed] ~= 2021.06' ], ) - -# OLDER VERSION (Before MeanEddy branch Jun 8 2021) - -# from setuptools import setup - -# def read(file): -# return open(file, 'r').read() - -# LONG_DESCRIPTION = read('README.md') -# LICENSE = read('LICENSE.txt') - -# setup( -# name='OceanLab', -# version='0.0.7', -# packages=['OceanLab'], -# include_package_data=True, -# description='Python functions for Physical Oceanography', -# long_description=LONG_DESCRIPTION, -# long_description_content_type='text/markdown', -# download_url = 'https://pypi.python.org/pypi/OceanLab', -# url='https://github.com/iuryt/OceanLab', -# author='Iury T. Simoes-Sousa', -# author_email='simoesiury@gmail.com', -# license=LICENSE, -# py_modules=['OA','EOF','DYN'], -# classifiers=[ -# "Programming Language :: Python :: 3", -# "Programming Language :: Python :: 3.6", -# "Programming Language :: Python :: 3.7", -# "License :: OSI Approved :: MIT License", -# "Operating System :: OS Independent", -# ], -# install_requires = [ -# 'seawater ~= 3.3', -# 'numpy ~= 1.18', -# ], -# ) \ No newline at end of file From 5cc9edecb86945af0d716968c62df8b3e2353761 Mon Sep 17 00:00:00 2001 From: dantecn Date: Wed, 9 Jun 2021 12:33:12 +0200 Subject: [PATCH 17/40] added underscore to function definition _parallel_client --- OceanLab/utils.py | 52 +++++++++++++++++++++++++---------------------- 1 file changed, 28 insertions(+), 24 deletions(-) diff --git a/OceanLab/utils.py b/OceanLab/utils.py index adb81f0..8e793a7 100644 --- a/OceanLab/utils.py +++ b/OceanLab/utils.py @@ -5,30 +5,7 @@ import dask from dask.distributed import Client, progress -##### Functions -# ============================================================================= -# 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 -#============================================================================= - +##### User functions #============================================================================= # NEAREST DISTANCE #============================================================================= @@ -156,3 +133,30 @@ def timefilter(prop,filtdays=60): 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 +#============================================================================= + + + From 33d8c01baa89886a76c14ee22ad6a9bcc357313f Mon Sep 17 00:00:00 2001 From: iury simoes-sousa Date: Fri, 11 Jun 2021 08:30:59 -0400 Subject: [PATCH 18/40] Update oa.py --- OceanLab/oa.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/OceanLab/oa.py b/OceanLab/oa.py index 3a47dc3..9853dcf 100644 --- a/OceanLab/oa.py +++ b/OceanLab/oa.py @@ -34,12 +34,7 @@ def vectoa(Xg,Yg,X,Y,U,V,corrlenx,corrleny,err,b=0): PSI - Streamfuction field matrix with MxN dimension. The dimension of the output is always the same of XC & YC - ====================================================================== - - PYTHON VERSION by: - Iury Sousa and Hélio Almeida - 30 May 2016 - Laboratório de Dinâmica Oceânica - IOUSP - ======================================================================''' + ''' # making sure that the input variables aren't changed xc,yc,x,y,u,v=Xg.copy(),Yg.copy(),X.copy(),Y.copy(),U.copy(),V.copy() From 322e34a97f127164efd6c915b23102d02a1e2fcc Mon Sep 17 00:00:00 2001 From: iury simoes-sousa Date: Fri, 11 Jun 2021 08:32:42 -0400 Subject: [PATCH 19/40] Update dyn.py --- OceanLab/dyn.py | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/OceanLab/dyn.py b/OceanLab/dyn.py index e817cde..cf923c4 100644 --- a/OceanLab/dyn.py +++ b/OceanLab/dyn.py @@ -65,10 +65,7 @@ def zeta(x,y,U,V): OUTPUT: ZETA = Relative vorticity field [s^-1] - - AUTHOR: - Iury T. Simoes-Sousa and Wandrey Watanabe - 29 Jun 2016 - Laboratório de Dinâmica Oceânica - IOUSP + ================================================================ ''' @@ -171,10 +168,6 @@ def psi2uv(x,y,psi): u = velocity zonal component [m s^-1] v = velocity meridional component [m s^-1] - AUTHOR: - Iury T. Simoes-Sousa and Wandrey Watanabe - May 2016 - Laboratório de Dinâmica Oceânica - IOUSP - Universidade de São Paulo ========================================================== ''' @@ -257,10 +250,6 @@ def eqmodes(N2,z,nm,lat,pmodes=False): 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 ============================================================ ''' From e8036bd19eb043cafa8152a70c0fec2dce7c7483 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Tue, 15 Jun 2021 13:15:41 -0300 Subject: [PATCH 20/40] teste ceof --- OceanLab/eof.py | 119 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index 5d8338c..33ea50e 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -121,3 +121,122 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): vi = (M.T+Mmean).T return vi + +#========================================= +# PERFORM COMPLEX EOF +#========================================= +from scipy.signal import hilbert +import scipy.linalg as la +import numpy as np +from dask import delayed +from scipy.signal import hilbert +import scipy.linalg as la +import numpy as np + +def ceof(data, nkp = 10): + ''' Complex (Hilbert) EOF + First written in MATLAB and found in Prof. Daniel J. Vimont webpage + (https://www.aos.wisc.edu/~dvimont/matlab/Stat_Tools/complex_eof.html) + + Translated to Python by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Mar/2021 + ============================================================================== + Inputs: + -> data: original data set [time, space] + -> nkp: number of modes to output (default = 10) + + ============================================================================== + Outputs: + -> lamda: eigenvalues + -> load: First nkp Complex Loadings or eigenvectors [space, nkp] + -> pcs: First nkp Complex Principal Components or amplitudes [time, nkp] + -> per: percent variance explained (real eigenvalues) + ============================================================================== + Version 2.0.0 by Felipe Vilela da Silva on 25/May/2021. + Now, it is possible to input data with NaN values. + The quantity of lambda is related the amount of not Nan values in data + ''' + load_real = np.zeros([data.shape[1], nkp])*np.nan + load_imag = np.zeros([data.shape[1], nkp])*np.nan + # It is necessary to remove the nan values of the matrix to solve the eigenvalue problem + nan_values = np.isnan(data[0,:]) # We can just look at each coordinate along a single time + data = data[:,~nan_values] # Then, we remove all these coordinates in all of the occurences + + ntim, npt = data.shape + + # Hilbert transform: input sequence x and returns a complex result of the same length + print('1: Performing Hilbert transform') + data_hilbert = hilbert(data) + # Compute the covariance matrix in the Hilbert transform + print('2: Computing covariance matrix') + c = delayed(np.dot)(data_hilbert.conjugate().T, data_hilbert).compute()/ntim + print('3: Solving the eigenvalue problem') + lamda, loadings = delayed(la.eig)(c).compute() # lamda: eigenvalue, loadings: eigenvectors + + l = lamda.conjugate().T; k = np.argsort(l) + lamda, loadings = np.flip(l[k]), np.fliplr(loadings[:,k]) + loadings = loadings[:,:nkp] + + per = lamda.real*100/np.sum(lamda.real) + pcs = np.dot(data_hilbert,loadings) + + load_real[~nan_values,:] = loadings.real.copy() + load_imag[~nan_values,:] = loadings.imag.copy() + load = load_real + 1j*load_imag + print('Done! \U0001F600') + + return lamda, load, pcs, per + +def amplitude_phase(evecs, amp): + ''' Complex (Hilbert) EOF + First written in MATLAB and found in the webpage below + (https://www.jsg.utexas.edu/fu/files/GEO391-W11-CEOF.pdf) + + Translated to Python by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Mar/2021 + =========================================================================== + Inputs: + -> evecs: First nkp Complex Loadings or eigenvectors [space, nkp] + -> amp: First nkp Complex Principal Components or amplitudes [time, nkp] + + =========================================================================== + Outputs: + -> sp_amp: Spatial amplitude [space, nkp] + -> sp_phase: Spatial phase [space, nkp] + -> t_amp: Temporal amplitude [time, nkp] + -> t_phase: Temporal phase [time, nkp] + =========================================================================== + ''' + # Spatial amplitude + sp_amp = pow(np.multiply(evecs, np.conj(evecs)),0.5) + theta = np.arctan2(evecs.imag, evecs.real) + # Spatial phase + sp_phase = np.divide(np.multiply(theta, 180), np.pi) + + # Temporal amplitude + t_amp = pow(np.multiply(amp, np.conj(amp)), 0.5) + # Temporal phase + phit = np.arctan2(amp.imag, amp.real) + t_phase = np.divide(np.multiply(phit, 180), np.pi) + + return sp_amp, sp_phase, t_amp, t_phase + +def reconstruct_ceof(amp, modes, n, day): + ''' Reconstrucion of CEOF modes individually following Majumder et al. (2019) + + Written by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Apr/2021 + =========================================================================== + Inputs: + -> amp: amplitude or coefficient of expansion [time, n] + -> modes: eigenvector or loading [lat, lon, n] + -> n: mode to be reconstructed. 0 is the first + -> day: day to be reconstructed. 0 is the first + =========================================================================== + Output: + -> Rec_ceof: Reconstruction of the CEOF field [time, lat, lon] + =========================================================================== + ''' + + # Majumder et al (2019) compute the reconstructed CEOF field as the real part of the multiplication between + # the coefficient of expansion (i.e., amplitude) and the complex conjugate of the loading (i.e., mode) + Rec_ceof = amp[day,n]*np.conj(modes[:,:,n]) + return Rec_ceof.real + From 92df601a844b3650bd222cb68765f29e574e431a Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Fri, 25 Jun 2021 14:52:29 -0300 Subject: [PATCH 21/40] Now, it is possible to input data with NaN values. The function organizes the data as time vs space and subtract the mean field in each coordinate. It returns all the variables related to CEOF in a DataArray. --- OceanLab/eof.py | 88 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 64 insertions(+), 24 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index 33ea50e..b821ecd 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -127,13 +127,11 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): #========================================= from scipy.signal import hilbert import scipy.linalg as la -import numpy as np from dask import delayed from scipy.signal import hilbert -import scipy.linalg as la -import numpy as np +import xarray as xr -def ceof(data, nkp = 10): +def ceof(lon, lat, data, nkp = 10): ''' Complex (Hilbert) EOF First written in MATLAB and found in Prof. Daniel J. Vimont webpage (https://www.aos.wisc.edu/~dvimont/matlab/Stat_Tools/complex_eof.html) @@ -141,31 +139,49 @@ def ceof(data, nkp = 10): Translated to Python by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Mar/2021 ============================================================================== Inputs: - -> data: original data set [time, space] + -> lon: longitudes (array) + -> lat: latitude (array) + -> data: original data set [time, lat, lon] + -> Note: the mean field in each coordinate is subtracted within the function + -> Thus, do not remove the time-mean before inputing in the fuction. -> nkp: number of modes to output (default = 10) ============================================================================== - Outputs: - -> lamda: eigenvalues - -> load: First nkp Complex Loadings or eigenvectors [space, nkp] - -> pcs: First nkp Complex Principal Components or amplitudes [time, nkp] + Outputs in a DataArray: -> per: percent variance explained (real eigenvalues) + -> modes: First nkp Complex Loadings or eigenvectors [lat, lon, nkp] + -> sp_amp: Spatial amplitude [lat, lon, nkp] + -> sp_phase: Spatial phase [lat, lon, nkp] + -> pcs: First nkp Complex Principal Components or amplitudes [time, nkp] + -> t_amp: Temporal amplitude [time, nkp] + -> t_phase: Temporal phase [time, nkp] ============================================================================== Version 2.0.0 by Felipe Vilela da Silva on 25/May/2021. Now, it is possible to input data with NaN values. - The quantity of lambda is related the amount of not Nan values in data - ''' - load_real = np.zeros([data.shape[1], nkp])*np.nan - load_imag = np.zeros([data.shape[1], nkp])*np.nan + ============================================================================== + Version 3.0.0 by Felipe Vilela da Silva on 17/Jun/2021. + New modifications: + (i) The function organizes the data as time vs space and subtract the mean field in each coordinate + (ii) It returns all the variables related to CEOF in a DataArray + ''' + # Organizing the data as time vs space + data_ceof = org_data_ceof(lon, lat, data) + # We need to remove the mean field (i.e., the trend) in each coordinate to + # evaluate the variability + data_ceof = data_ceof - np.nanmean(data_ceof,axis=0) + + # The variables below are useful later + load_real = np.zeros([data_ceof.shape[1], nkp])*np.nan + load_imag = np.zeros([data_ceof.shape[1], nkp])*np.nan # It is necessary to remove the nan values of the matrix to solve the eigenvalue problem - nan_values = np.isnan(data[0,:]) # We can just look at each coordinate along a single time - data = data[:,~nan_values] # Then, we remove all these coordinates in all of the occurences + nan_values = np.isnan(data_ceof[0,:]) # We can just look at each coordinate along a single time + data_ceof = data_ceof[:,~nan_values] # Then, we remove all these coordinates in all of the occurences - ntim, npt = data.shape + ntim, npt = data_ceof.shape # Hilbert transform: input sequence x and returns a complex result of the same length print('1: Performing Hilbert transform') - data_hilbert = hilbert(data) + data_hilbert = hilbert(data_ceof) # Compute the covariance matrix in the Hilbert transform print('2: Computing covariance matrix') c = delayed(np.dot)(data_hilbert.conjugate().T, data_hilbert).compute()/ntim @@ -175,16 +191,39 @@ def ceof(data, nkp = 10): l = lamda.conjugate().T; k = np.argsort(l) lamda, loadings = np.flip(l[k]), np.fliplr(loadings[:,k]) loadings = loadings[:,:nkp] + # In case there were nan values in the orginal data, we need to perform the approach below: + load_real[~nan_values,:] = loadings.real.copy() + load_imag[~nan_values,:] = loadings.imag.copy() + load = load_real + 1j*load_imag + modes = load.reshape((len(lat),len(lon), nkp)) per = lamda.real*100/np.sum(lamda.real) + per = per[:nkp].copy() pcs = np.dot(data_hilbert,loadings) - load_real[~nan_values,:] = loadings.real.copy() - load_imag[~nan_values,:] = loadings.imag.copy() - load = load_real + 1j*load_imag + sp_amp, sp_phase, t_amp, t_phase = amplitude_phase(load, pcs) + sp_amp = sp_amp.reshape((len(lat),len(lon), nkp)) + sp_phase = sp_phase.reshape((len(lat),len(lon), nkp)) + print('Done! \U0001F600') - return lamda, load, pcs, per + dims = ["lat", "lon", "nkp", "time"] + ds = xr.Dataset({"per":(dims[2], per),"modes":(dims[:-1], modes),"sp_amp":(dims[:-1], sp_amp), + "sp_phase":(dims[:-1], sp_phase),"pcs":(dims[-2:][::-1], pcs),"t_amp":(dims[-2:][::-1], t_amp), + "t_phase":(dims[-2:][::-1], t_phase)}, + coords={"lat":(dims[0], lat), "lon":(dims[1], lon), "nkp":(dims[2], np.arange(nkp)), + "time":(dims[3], np.arange(len(data_ceof)))}) + + return ds + +def org_data_ceof(lon, lat, data): + data_ceof = np.zeros((len(data), len(lon)*len(lat)))*np.nan + k = 0 + for i in range(len(lat)): + for j in range(len(lon)): + data_ceof[:, k] = data[:, i, j] + k += 1 + return data_ceof def amplitude_phase(evecs, amp): ''' Complex (Hilbert) EOF @@ -219,12 +258,13 @@ def amplitude_phase(evecs, amp): return sp_amp, sp_phase, t_amp, t_phase -def reconstruct_ceof(amp, modes, n, day): +def reconstruct_ceof(data_mean, amp, modes, n, day): ''' Reconstrucion of CEOF modes individually following Majumder et al. (2019) Written by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Apr/2021 =========================================================================== Inputs: + -> data_mean: time-mean of the original data [lat, lon] (e.g., np.nanmean(data,axis=0)) -> amp: amplitude or coefficient of expansion [time, n] -> modes: eigenvector or loading [lat, lon, n] -> n: mode to be reconstructed. 0 is the first @@ -238,5 +278,5 @@ def reconstruct_ceof(amp, modes, n, day): # Majumder et al (2019) compute the reconstructed CEOF field as the real part of the multiplication between # the coefficient of expansion (i.e., amplitude) and the complex conjugate of the loading (i.e., mode) Rec_ceof = amp[day,n]*np.conj(modes[:,:,n]) - return Rec_ceof.real - + Rec_return = Rec_ceof.real + data_mean + return Rec_return From 8db085edd3abeee196b970a80ee96900a320d711 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Fri, 25 Jun 2021 15:04:20 -0300 Subject: [PATCH 22/40] Update README.md I added the description for the ceof and reconstruct_ceof modules. --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 663bc46..aae0180 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,8 @@ Check `examples` folder in our [github repository](github.com/iuryt/OceanLab/exa - **EOF** - *eoft()*: Calculates the Empirical Orthogonal Functions; - *my_eof_interp()*: Fillgaps on matrix based on EOFs (translated from Cesar Rocha Matlab version); + - *ceof()*: Performs the Complex (or Hilbert) Empirical Orthogonal Functions decomposition; + - *reconstruct_ceof()*: Reconstructs the CEOF modes individually; - **UTILS** - *argdistnear()*: Searchs the position of the closest points in an array to a reference point; - *meaneddy()*: Performs an eddy-mean decomposing with a low-pass filter; From 780d05a14ef4d005112d251517d8507ca76e8cf7 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Fri, 25 Jun 2021 15:31:48 -0300 Subject: [PATCH 23/40] Update __init.py__ I added the ceof, reconstruct_ceof modules in the init file --- OceanLab/__init.py__ | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/OceanLab/__init.py__ b/OceanLab/__init.py__ index 20dbed1..a1a2eae 100644 --- a/OceanLab/__init.py__ +++ b/OceanLab/__init.py__ @@ -1,5 +1,5 @@ from .dyn import zeta, vmode_amp, psi2uv, eqmodes, vmodes -from .eof import my_eof_interp, eoft +from .eof import my_eof_interp, eoft, ceof, reconstruct_ceof from .oa import scaloa, vectoa from .utils import argdistnear, meaneddy From 927f9bbb8c3af8a0b4f107501a6db22a9e9a885a Mon Sep 17 00:00:00 2001 From: Dante Campagnoli Napolitano Date: Mon, 12 Jul 2021 16:41:38 +0200 Subject: [PATCH 24/40] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 663bc46..4a86840 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,7 @@ Check `examples` folder in our [github repository](github.com/iuryt/OceanLab/exa - *my_eof_interp()*: Fillgaps on matrix based on EOFs (translated from Cesar Rocha Matlab version); - **UTILS** - *argdistnear()*: Searchs the position of the closest points in an array to a reference point; - - *meaneddy()*: Performs an eddy-mean decomposing with a low-pass filter; + - *meaneddy()*: Performs an eddy-mean decomposition with a low-pass filter; ## Contributors From cfc9febbfed83890a3ad425e23f9bfa032f701b7 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Mon, 12 Jul 2021 13:00:09 -0300 Subject: [PATCH 25/40] Update eof.py Following Ryukamusa request about "Imports should be placed in the top of each file" --- OceanLab/eof.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index b821ecd..a65982c 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -1,4 +1,8 @@ import numpy as np +import scipy.linalg as la +from dask import delayed +from scipy.signal import hilbert +import xarray as xr # functions #========================================= @@ -125,12 +129,6 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): #========================================= # PERFORM COMPLEX EOF #========================================= -from scipy.signal import hilbert -import scipy.linalg as la -from dask import delayed -from scipy.signal import hilbert -import xarray as xr - def ceof(lon, lat, data, nkp = 10): ''' Complex (Hilbert) EOF First written in MATLAB and found in Prof. Daniel J. Vimont webpage From ab76e7d50b4f07e65b7a6ad37383bd8f89763d43 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Mon, 12 Jul 2021 13:08:48 -0300 Subject: [PATCH 26/40] Removed authorship comments Following dantecn request --- OceanLab/eof.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index a65982c..ae2e101 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -133,8 +133,7 @@ def ceof(lon, lat, data, nkp = 10): ''' Complex (Hilbert) EOF First written in MATLAB and found in Prof. Daniel J. Vimont webpage (https://www.aos.wisc.edu/~dvimont/matlab/Stat_Tools/complex_eof.html) - - Translated to Python by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Mar/2021 + Version 1.0.0 on 15/Mar/2021 ============================================================================== Inputs: -> lon: longitudes (array) @@ -154,10 +153,10 @@ def ceof(lon, lat, data, nkp = 10): -> t_amp: Temporal amplitude [time, nkp] -> t_phase: Temporal phase [time, nkp] ============================================================================== - Version 2.0.0 by Felipe Vilela da Silva on 25/May/2021. + Version 2.0.0 on 25/May/2021. Now, it is possible to input data with NaN values. ============================================================================== - Version 3.0.0 by Felipe Vilela da Silva on 17/Jun/2021. + Version 3.0.0 on 17/Jun/2021. New modifications: (i) The function organizes the data as time vs space and subtract the mean field in each coordinate (ii) It returns all the variables related to CEOF in a DataArray @@ -228,7 +227,6 @@ def amplitude_phase(evecs, amp): First written in MATLAB and found in the webpage below (https://www.jsg.utexas.edu/fu/files/GEO391-W11-CEOF.pdf) - Translated to Python by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Mar/2021 =========================================================================== Inputs: -> evecs: First nkp Complex Loadings or eigenvectors [space, nkp] @@ -259,7 +257,6 @@ def amplitude_phase(evecs, amp): def reconstruct_ceof(data_mean, amp, modes, n, day): ''' Reconstrucion of CEOF modes individually following Majumder et al. (2019) - Written by Felipe Vilela da Silva @ IMAS, UTas, AU on 15/Apr/2021 =========================================================================== Inputs: -> data_mean: time-mean of the original data [lat, lon] (e.g., np.nanmean(data,axis=0)) From aac365f4df1d63309c9c4740adb50b3c88ca8da5 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Mon, 12 Jul 2021 13:18:55 -0300 Subject: [PATCH 27/40] Update eof.org_data_ceof() Now, org_data_ceof works faster and returns a DataArray --- OceanLab/eof.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index ae2e101..3c8abe5 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -165,7 +165,7 @@ def ceof(lon, lat, data, nkp = 10): data_ceof = org_data_ceof(lon, lat, data) # We need to remove the mean field (i.e., the trend) in each coordinate to # evaluate the variability - data_ceof = data_ceof - np.nanmean(data_ceof,axis=0) + data_ceof = data_ceof - data_ceof.mean('time') # The variables below are useful later load_real = np.zeros([data_ceof.shape[1], nkp])*np.nan @@ -214,12 +214,11 @@ def ceof(lon, lat, data, nkp = 10): return ds def org_data_ceof(lon, lat, data): - data_ceof = np.zeros((len(data), len(lon)*len(lat)))*np.nan - k = 0 - for i in range(len(lat)): - for j in range(len(lon)): - data_ceof[:, k] = data[:, i, j] - k += 1 + ''' Version 2.0.0 on 12/Jul/2021: Now, org_data_ceof works faster and returns a DataArray ''' + dims = ["time", "lat", "lon"] + datxarray = a = xr.Dataset({"data_latlon": (dims, data)}, + coords={'lat':(dims[1], lat), 'lon':(dims[2], lon)}) + data_ceof = datxarray.stack(lat_lon=("lat", "lon")).data_latlon return data_ceof def amplitude_phase(evecs, amp): From 3bc6cb78f64dca1bd366f1452190038be54fd2a2 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Mon, 12 Jul 2021 14:15:14 -0300 Subject: [PATCH 28/40] Update inputs and outputs labels Followed @dantecn suggestion. --- OceanLab/eof.py | 102 ++++++++++++++++++++++++------------------------ 1 file changed, 52 insertions(+), 50 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index 3c8abe5..14959cc 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -131,27 +131,27 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): #========================================= def ceof(lon, lat, data, nkp = 10): ''' Complex (Hilbert) EOF + Note: the mean field in each coordinate is subtracted within the function. + Do not remove the time-mean before inputing in the fuction. + First written in MATLAB and found in Prof. Daniel J. Vimont webpage (https://www.aos.wisc.edu/~dvimont/matlab/Stat_Tools/complex_eof.html) Version 1.0.0 on 15/Mar/2021 ============================================================================== - Inputs: - -> lon: longitudes (array) - -> lat: latitude (array) - -> data: original data set [time, lat, lon] - -> Note: the mean field in each coordinate is subtracted within the function - -> Thus, do not remove the time-mean before inputing in the fuction. - -> nkp: number of modes to output (default = 10) - - ============================================================================== - Outputs in a DataArray: - -> per: percent variance explained (real eigenvalues) - -> modes: First nkp Complex Loadings or eigenvectors [lat, lon, nkp] - -> sp_amp: Spatial amplitude [lat, lon, nkp] - -> sp_phase: Spatial phase [lat, lon, nkp] - -> pcs: First nkp Complex Principal Components or amplitudes [time, nkp] - -> t_amp: Temporal amplitude [time, nkp] - -> t_phase: Temporal phase [time, nkp] + INPUT: + lon = longitudes (array) + lat = latitude (array) + data = original data set [time, lat, lon] + nkp = number of modes to return (default = 10) + + OUTPUT: + per = percent variance explained (real eigenvalues) + modes = first nkp complex loadings or eigenvectors [lat, lon, nkp] + SpAmp = spatial amplitude [lat, lon, nkp] + SpPhase = spatial phase [lat, lon, nkp] + pcs = first nkp complex principal components or amplitudes [time, nkp] + TAmp = temporal amplitude [time, nkp] + TPhase = temporal phase [time, nkp] ============================================================================== Version 2.0.0 on 25/May/2021. Now, it is possible to input data with NaN values. @@ -205,9 +205,9 @@ def ceof(lon, lat, data, nkp = 10): print('Done! \U0001F600') dims = ["lat", "lon", "nkp", "time"] - ds = xr.Dataset({"per":(dims[2], per),"modes":(dims[:-1], modes),"sp_amp":(dims[:-1], sp_amp), - "sp_phase":(dims[:-1], sp_phase),"pcs":(dims[-2:][::-1], pcs),"t_amp":(dims[-2:][::-1], t_amp), - "t_phase":(dims[-2:][::-1], t_phase)}, + ds = xr.Dataset({"per":(dims[2], per),"modes":(dims[:-1], modes),"SpAmp":(dims[:-1], sp_amp), + "SpPhase":(dims[:-1], sp_phase),"pcs":(dims[-2:][::-1], pcs),"TAmp":(dims[-2:][::-1], t_amp), + "TPhase":(dims[-2:][::-1], t_phase)}, coords={"lat":(dims[0], lat), "lon":(dims[1], lon), "nkp":(dims[2], np.arange(nkp)), "time":(dims[3], np.arange(len(data_ceof)))}) @@ -227,50 +227,52 @@ def amplitude_phase(evecs, amp): (https://www.jsg.utexas.edu/fu/files/GEO391-W11-CEOF.pdf) =========================================================================== - Inputs: - -> evecs: First nkp Complex Loadings or eigenvectors [space, nkp] - -> amp: First nkp Complex Principal Components or amplitudes [time, nkp] - - =========================================================================== - Outputs: - -> sp_amp: Spatial amplitude [space, nkp] - -> sp_phase: Spatial phase [space, nkp] - -> t_amp: Temporal amplitude [time, nkp] - -> t_phase: Temporal phase [time, nkp] + INPUT: + evecs = first nkp complex loadings or eigenvectors [lat, lon, nkp] + amp = first nkp complex principal components or amplitudes [time, nkp] + + OUTPUT: + SpAmp = spatial amplitude [lat, lon, nkp] + SpPhase = spatial phase [lat, lon, nkp] + TAmp = temporal amplitude [time, nkp] + TPhase = temporal phase [time, nkp] =========================================================================== ''' # Spatial amplitude - sp_amp = pow(np.multiply(evecs, np.conj(evecs)),0.5) + SpAmp = pow(np.multiply(evecs, np.conj(evecs)),0.5) theta = np.arctan2(evecs.imag, evecs.real) # Spatial phase - sp_phase = np.divide(np.multiply(theta, 180), np.pi) + SpPhase = np.divide(np.multiply(theta, 180), np.pi) # Temporal amplitude - t_amp = pow(np.multiply(amp, np.conj(amp)), 0.5) + TAmp = pow(np.multiply(amp, np.conj(amp)), 0.5) # Temporal phase phit = np.arctan2(amp.imag, amp.real) - t_phase = np.divide(np.multiply(phit, 180), np.pi) + TPhase = np.divide(np.multiply(phit, 180), np.pi) - return sp_amp, sp_phase, t_amp, t_phase + return SpAmp, SpPhase, TAmp, TPhase -def reconstruct_ceof(data_mean, amp, modes, n, day): - ''' Reconstrucion of CEOF modes individually following Majumder et al. (2019) +def reconstruct_ceof(DataMean, amp, modes, n, day): + ''' Reconstrucion of daily CEOF modes individually similar to Majumder et al. (2019). + Here, the mean field in each coordinate is added within the function. + Besides, each mode is reconstructed individually, instead of computing the sum of the + reconstruction of different modes. - =========================================================================== - Inputs: - -> data_mean: time-mean of the original data [lat, lon] (e.g., np.nanmean(data,axis=0)) - -> amp: amplitude or coefficient of expansion [time, n] - -> modes: eigenvector or loading [lat, lon, n] - -> n: mode to be reconstructed. 0 is the first - -> day: day to be reconstructed. 0 is the first - =========================================================================== - Output: - -> Rec_ceof: Reconstruction of the CEOF field [time, lat, lon] - =========================================================================== + ========================================================================================= + INPUT: + DataMean = time-mean of the original data [lat, lon] (e.g., np.nanmean(data,axis=0)) + amp = principal components or amplitudes [time, nkp] + mode = eigenvectors or loadings [lat, lon, nkp] + n = mode of variability to be reconstructed + day = day to be reconstructed. + + OUTPUT: + RecCEOF = reconstruction of a CEOF mode on a chosen day [lat, lon] + ========================================================================================= ''' # Majumder et al (2019) compute the reconstructed CEOF field as the real part of the multiplication between # the coefficient of expansion (i.e., amplitude) and the complex conjugate of the loading (i.e., mode) Rec_ceof = amp[day,n]*np.conj(modes[:,:,n]) - Rec_return = Rec_ceof.real + data_mean - return Rec_return + RecCEOF = Rec_ceof.real + DataMean + return RecCEOF From 38bd15ba4134e3a24b5b25e09684355bb9e6ba52 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Tue, 13 Jul 2021 11:24:39 -0300 Subject: [PATCH 29/40] Update eof.org_data_ceof() I just realized a typo. --- OceanLab/eof.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index 14959cc..2dfd772 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -216,8 +216,8 @@ def ceof(lon, lat, data, nkp = 10): def org_data_ceof(lon, lat, data): ''' Version 2.0.0 on 12/Jul/2021: Now, org_data_ceof works faster and returns a DataArray ''' dims = ["time", "lat", "lon"] - datxarray = a = xr.Dataset({"data_latlon": (dims, data)}, - coords={'lat':(dims[1], lat), 'lon':(dims[2], lon)}) + datxarray = xr.Dataset({"data_latlon": (dims, data)}, + coords={'lat':(dims[1], lat), 'lon':(dims[2], lon)}) data_ceof = datxarray.stack(lat_lon=("lat", "lon")).data_latlon return data_ceof From 317b7d23e6273aa99acace710ce568e12e77173c Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Wed, 14 Jul 2021 14:00:16 -0300 Subject: [PATCH 30/40] Removal of versions in the ceof description I just followed dantecn suggestion --- OceanLab/eof.py | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index 2dfd772..eb5b8a6 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -132,11 +132,12 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): def ceof(lon, lat, data, nkp = 10): ''' Complex (Hilbert) EOF Note: the mean field in each coordinate is subtracted within the function. - Do not remove the time-mean before inputing in the fuction. + Do not subtract the time-mean field before inputing. + NaN values are removed in the algorithm. + The user can input the data as it is. First written in MATLAB and found in Prof. Daniel J. Vimont webpage (https://www.aos.wisc.edu/~dvimont/matlab/Stat_Tools/complex_eof.html) - Version 1.0.0 on 15/Mar/2021 ============================================================================== INPUT: lon = longitudes (array) @@ -145,21 +146,15 @@ def ceof(lon, lat, data, nkp = 10): nkp = number of modes to return (default = 10) OUTPUT: + The variables below return inside a DataArray. per = percent variance explained (real eigenvalues) modes = first nkp complex loadings or eigenvectors [lat, lon, nkp] SpAmp = spatial amplitude [lat, lon, nkp] SpPhase = spatial phase [lat, lon, nkp] pcs = first nkp complex principal components or amplitudes [time, nkp] TAmp = temporal amplitude [time, nkp] - TPhase = temporal phase [time, nkp] - ============================================================================== - Version 2.0.0 on 25/May/2021. - Now, it is possible to input data with NaN values. + TPhase = temporal phase [time, nkp] ============================================================================== - Version 3.0.0 on 17/Jun/2021. - New modifications: - (i) The function organizes the data as time vs space and subtract the mean field in each coordinate - (ii) It returns all the variables related to CEOF in a DataArray ''' # Organizing the data as time vs space data_ceof = org_data_ceof(lon, lat, data) @@ -172,7 +167,7 @@ def ceof(lon, lat, data, nkp = 10): load_imag = np.zeros([data_ceof.shape[1], nkp])*np.nan # It is necessary to remove the nan values of the matrix to solve the eigenvalue problem nan_values = np.isnan(data_ceof[0,:]) # We can just look at each coordinate along a single time - data_ceof = data_ceof[:,~nan_values] # Then, we remove all these coordinates in all of the occurences + data_ceof = data_ceof[:,~nan_values] # Then, we remove all these coordinates in all of the occurences ntim, npt = data_ceof.shape @@ -214,7 +209,6 @@ def ceof(lon, lat, data, nkp = 10): return ds def org_data_ceof(lon, lat, data): - ''' Version 2.0.0 on 12/Jul/2021: Now, org_data_ceof works faster and returns a DataArray ''' dims = ["time", "lat", "lon"] datxarray = xr.Dataset({"data_latlon": (dims, data)}, coords={'lat':(dims[1], lat), 'lon':(dims[2], lon)}) From a32944621a19c4b0ce4cc860031bf82f5c85b4f3 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Mon, 2 Aug 2021 15:38:19 -0300 Subject: [PATCH 31/40] Update auxiliary functions "_" at the beginning of auxiliary functions. --- OceanLab/eof.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index eb5b8a6..f33fcc2 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -130,7 +130,7 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): # PERFORM COMPLEX EOF #========================================= def ceof(lon, lat, data, nkp = 10): - ''' Complex (Hilbert) EOF + ''' Complex (Hilbert) EOF to detect propagating features: waves, meanders, etc. Note: the mean field in each coordinate is subtracted within the function. Do not subtract the time-mean field before inputing. NaN values are removed in the algorithm. @@ -157,7 +157,7 @@ def ceof(lon, lat, data, nkp = 10): ============================================================================== ''' # Organizing the data as time vs space - data_ceof = org_data_ceof(lon, lat, data) + data_ceof = _org_data_ceof(lon, lat, data) # We need to remove the mean field (i.e., the trend) in each coordinate to # evaluate the variability data_ceof = data_ceof - data_ceof.mean('time') @@ -193,7 +193,7 @@ def ceof(lon, lat, data, nkp = 10): per = per[:nkp].copy() pcs = np.dot(data_hilbert,loadings) - sp_amp, sp_phase, t_amp, t_phase = amplitude_phase(load, pcs) + sp_amp, sp_phase, t_amp, t_phase = _amplitude_phase(load, pcs) sp_amp = sp_amp.reshape((len(lat),len(lon), nkp)) sp_phase = sp_phase.reshape((len(lat),len(lon), nkp)) @@ -208,14 +208,14 @@ def ceof(lon, lat, data, nkp = 10): return ds -def org_data_ceof(lon, lat, data): +def _org_data_ceof(lon, lat, data): dims = ["time", "lat", "lon"] datxarray = xr.Dataset({"data_latlon": (dims, data)}, coords={'lat':(dims[1], lat), 'lon':(dims[2], lon)}) data_ceof = datxarray.stack(lat_lon=("lat", "lon")).data_latlon return data_ceof -def amplitude_phase(evecs, amp): +def _amplitude_phase(evecs, amp): ''' Complex (Hilbert) EOF First written in MATLAB and found in the webpage below (https://www.jsg.utexas.edu/fu/files/GEO391-W11-CEOF.pdf) From 2c3c5898473c08c6a85376ee0394868877875a06 Mon Sep 17 00:00:00 2001 From: vsilvafelipe Date: Wed, 4 Aug 2021 15:42:13 -0300 Subject: [PATCH 32/40] Create client within ceof I added the parallel option to build a standard client within the ceof function. In case the user doesn't want a standard client, he/she can turn the option to False and build the client with other settings. --- OceanLab/eof.py | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/OceanLab/eof.py b/OceanLab/eof.py index f33fcc2..3461cde 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -3,6 +3,7 @@ from dask import delayed from scipy.signal import hilbert import xarray as xr +from dask.distributed import Client, LocalCluster # functions #========================================= @@ -129,7 +130,7 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): #========================================= # PERFORM COMPLEX EOF #========================================= -def ceof(lon, lat, data, nkp = 10): +def ceof(lon, lat, data, nkp = 10, parallel = True): ''' Complex (Hilbert) EOF to detect propagating features: waves, meanders, etc. Note: the mean field in each coordinate is subtracted within the function. Do not subtract the time-mean field before inputing. @@ -140,22 +141,29 @@ def ceof(lon, lat, data, nkp = 10): (https://www.aos.wisc.edu/~dvimont/matlab/Stat_Tools/complex_eof.html) ============================================================================== INPUT: - lon = longitudes (array) - lat = latitude (array) - data = original data set [time, lat, lon] - nkp = number of modes to return (default = 10) + lon = longitudes (array) + lat = latitude (array) + data = original data set [time, lat, lon] + nkp = number of modes to return (default = 10) + parallel = create a standard client kernel for parallel computing + [switch parallel to False, in case you created your own client] OUTPUT: The variables below return inside a DataArray. - per = percent variance explained (real eigenvalues) - modes = first nkp complex loadings or eigenvectors [lat, lon, nkp] - SpAmp = spatial amplitude [lat, lon, nkp] - SpPhase = spatial phase [lat, lon, nkp] - pcs = first nkp complex principal components or amplitudes [time, nkp] - TAmp = temporal amplitude [time, nkp] - TPhase = temporal phase [time, nkp] + per = percent variance explained (real eigenvalues) + modes = first nkp complex loadings or eigenvectors [lat, lon, nkp] + SpAmp = spatial amplitude [lat, lon, nkp] + SpPhase = spatial phase [lat, lon, nkp] + pcs = first nkp complex principal components or amplitudes [time, nkp] + TAmp = temporal amplitude [time, nkp] + TPhase = temporal phase [time, nkp] ============================================================================== ''' + # Configure client for parallel computing + if parallel: + cluster = LocalCluster() + client = Client(cluster) + # Organizing the data as time vs space data_ceof = _org_data_ceof(lon, lat, data) # We need to remove the mean field (i.e., the trend) in each coordinate to From ea4c132c6fa161cf8e9ecccc1a9c499aacb0694d Mon Sep 17 00:00:00 2001 From: Dante Campagnoli Napolitano Date: Tue, 19 Oct 2021 10:31:36 +0200 Subject: [PATCH 33/40] Update README.md fix link to examples (master) directory --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 6dd1f34..31e5670 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ Package of Python scripts for Oceanography (Python +3.6) ## Code Example -Check `examples` folder in our [github repository](github.com/iuryt/OceanLab/examples). +Check `examples` folder in our [github repository](github.com/OceanLabPy/OceanLab/tree/master/examples). ## Installation From ca4f2e1c5781a93f78d7a242bf2e2edbae241265 Mon Sep 17 00:00:00 2001 From: Dante Campagnoli Napolitano Date: Tue, 19 Oct 2021 10:32:58 +0200 Subject: [PATCH 34/40] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 31e5670..9b57986 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ Package of Python scripts for Oceanography (Python +3.6) ## Code Example -Check `examples` folder in our [github repository](github.com/OceanLabPy/OceanLab/tree/master/examples). +Check `examples` folder in our [github repository](tree/master/examples). ## Installation From 57cd9456ee7c96c4ab96c87b3b1f2f5a141c562d Mon Sep 17 00:00:00 2001 From: Dante Campagnoli Napolitano Date: Tue, 19 Oct 2021 10:34:04 +0200 Subject: [PATCH 35/40] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 9b57986..31e5670 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ Package of Python scripts for Oceanography (Python +3.6) ## Code Example -Check `examples` folder in our [github repository](tree/master/examples). +Check `examples` folder in our [github repository](github.com/OceanLabPy/OceanLab/tree/master/examples). ## Installation From 7f3f6c9da3b7b5ad2c6de9d46b36451f590197e9 Mon Sep 17 00:00:00 2001 From: Dante Campagnoli Napolitano Date: Tue, 19 Oct 2021 10:36:10 +0200 Subject: [PATCH 36/40] Update README.md fixed link to examples --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 31e5670..8153edd 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ Package of Python scripts for Oceanography (Python +3.6) ## Code Example -Check `examples` folder in our [github repository](github.com/OceanLabPy/OceanLab/tree/master/examples). +Check `examples` folder in our [github repository](../../tree/master/examples). ## Installation From 4cc0b8dc346e8c7d57659e7a99456e88c7adecc7 Mon Sep 17 00:00:00 2001 From: Ou Wang Date: Thu, 31 Aug 2023 15:45:08 -0700 Subject: [PATCH 37/40] Fix a typo --- examples/DYN_EOF.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/DYN_EOF.ipynb b/examples/DYN_EOF.ipynb index 7feb8d1..6f04e4f 100644 --- a/examples/DYN_EOF.ipynb +++ b/examples/DYN_EOF.ipynb @@ -284,7 +284,7 @@ } ], "source": [ - "n = 100\n", + "q = 100\n", "xi,yi = np.meshgrid(np.linspace(0,15,q),np.linspace(0,15,q))\n", "\n", "vorti = DYN.zeta(xi,yi,lambdify((x,y),u)(xi,yi),lambdify((x,y),v)(xi,yi))\n", From 7b98af8fc654e5cfe163485548b0a1ed162773b4 Mon Sep 17 00:00:00 2001 From: Felipe Vilela da Silva Date: Fri, 13 Dec 2024 14:27:46 +1100 Subject: [PATCH 38/40] Update dyn.py I added the functions of curvature and surface ageostrophic velocity --- OceanLab/dyn.py | 67 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/OceanLab/dyn.py b/OceanLab/dyn.py index cf923c4..9b9a469 100644 --- a/OceanLab/dyn.py +++ b/OceanLab/dyn.py @@ -2,6 +2,73 @@ import numpy as np import seawater as sw +# =========================================================== +# SURFACE AGEOSTROPHIC VELOCITY +# =========================================================== +def curvature(ssh_grad, dxt, dyt): + ''' Curvature of the sea surface height contours + ============================================================================== + INPUT: + SSH = sea surface height gradient [Y,X] + dxt = zonal distance between the grid cells in meter [Y,X] + dyt = meridional distance between the grid cells in meter [Y,X] + + OUTPUT: + k = curvature of the flow + ============================================================================== + ''' + + ssh_grad_y, ssh_grad_x = ssh_grad[0], ssh_grad[1] + + # Terms in the norminator + nominator_1st = (np.gradient(ssh_grad_y/dyt)[0]/dyt)*(ssh_grad_x/dxt)**2 + nominator_2nd = (np.gradient(ssh_grad_x/dxt)[1]/dxt)*(ssh_grad_y/dyt)**2 + nominator_3rd = 2*(np.gradient(ssh_grad_y/dyt)[1]/dxt)*(ssh_grad_x/dxt)*(ssh_grad_y/dyt) + + # Terms in the denominator + denominator_1st = (ssh_grad_x/dxt)**2 + denominator_2nd = (ssh_grad_y/dyt)**2 + denominator = pow(denominator_1st + denominator_2nd, 3/2) + + k = (-1*nominator_1st - 1*nominator_2nd + 2*nominator_3rd)/denominator + + return k + +def ageostrophic_vel(LON, LAT, SSH): + ''' Surface ageostrophic velocity from gradient wind balance through the + combination of the surface geostrophic velocity and curvature of the flow + ============================================================================== + INPUT: + LON = longitudes [Y,X] + LAT = latitude [Y,X] + SSH = sea surface height in meter [Y,X] + + OUTPUT: + u_ag_vel = zonal component of the surface ageostrophic velocity + v_ag_vel = meridional component of the surface ageostrophic velocity + ============================================================================== + ''' + g = 9.8 + f = sw.f(LAT) + [trash,dlonx] = np.gradient(LON) + dx = dlonx*np.cos(LAT*np.pi/180)*111195. + [dlaty,trash] = np.gradient(LAT) + dy = dlaty*111195. + SSH_gradient = np.gradient(SSH) + k = curvature(SSH_gradient, dx, dy) + + ## Estimating the surface geostrophic velocity from SSH + u_geo_vel = (-g/f)*(SSH_gradient[0]/dy) + v_geo_vel = (g/f)*(SSH_gradient[1]/dx) + ## Reconstructing the surface velocity from the combinantion of the geostrophic velocity and curvature of the flow + u_total = (2*u_geo_vel)/(1 + np.sqrt(1 + 4*k*np.sqrt(u_geo_vel**2 + v_geo_vel**2)/f)) + v_total = (2*v_geo_vel)/(1 + np.sqrt(1 + 4*k*np.sqrt(u_geo_vel**2 + v_geo_vel**2)/f)) + ## Computing the surface ageostrophic velocity as the difference between the total and geostrophic fields + u_ag_vel = u_total - u_geo_vel + v_ag_vel = v_total - v_geo_vel + + return u_ag_vel, v_ag_vel + # =========================================================== # DYNAMICAL MODES AMPLITUDE # =========================================================== From b40b7bc3e3c8daeb243ab929bac43c29afe98e4e Mon Sep 17 00:00:00 2001 From: Felipe Vilela-Silva Date: Wed, 5 Mar 2025 16:19:20 +0100 Subject: [PATCH 39/40] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 8153edd..4d95112 100644 --- a/README.md +++ b/README.md @@ -33,4 +33,4 @@ Check `examples` folder in our [github repository](../../tree/master/examples). ## Contributors -Everyone can contribute to this code. Some of functions were based on Filipe Fernandes or Cesar Rocha functions and some of them were created with help of Dante C. Napolitano, Hélio M. R. Almeida and Wandrey Watanabe at Ocean Dynamics Lab of University of São Paulo (USP). +Everyone can contribute to this code. Some of functions were based on Filipe Fernandes or Cesar Rocha functions and some of them were created with help of Dante C. Napolitano, Hélio M. R. Almeida, Wandrey Watanabe, and Felipe Vilela-Silva at Ocean Dynamics Lab of University of São Paulo (USP). From 4953af7d80d4376fcf06a36a1eaa3708111fa6a2 Mon Sep 17 00:00:00 2001 From: DavidVadnais Date: Thu, 22 May 2025 00:32:00 +0000 Subject: [PATCH 40/40] Remove unnecessary semicolons at end of Python statements --- OceanLab/dyn.py | 2 +- OceanLab/eof.py | 7 ++++--- OceanLab/oa.py | 10 +++++----- OceanLab/utils.py | 3 +-- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/OceanLab/dyn.py b/OceanLab/dyn.py index 9b9a469..761b994 100644 --- a/OceanLab/dyn.py +++ b/OceanLab/dyn.py @@ -396,7 +396,7 @@ def eqmodes(N2,z,nm,lat,pmodes=False): # F_j(z)=-g.he_j d/dzS_j # #### - if pmodes==True: + if pmodes: Fi=np.zeros(Si.shape) for i in np.arange(nm): Fi[1:-1,i] =\ diff --git a/OceanLab/eof.py b/OceanLab/eof.py index 3461cde..2f1e922 100644 --- a/OceanLab/eof.py +++ b/OceanLab/eof.py @@ -44,7 +44,7 @@ def eoft(trmat,nm=None): amp = np.dot(evecs_norm.T,trmat) #if was chosen a number of eigenvectors - if nm!=None: + if nm is not None: evecs_norm = evecs_norm[:,:nm] evals_perc = evals_perc[:nm] amp = amp[:nm,:] @@ -117,7 +117,7 @@ def my_eof_interp(M,nmodes,errmin=1e-15,repmax=None): print(err) if err < errmin: break - if repmax!=None: + if repmax is not None: if rep>=repmax: break rep += 1 @@ -188,7 +188,8 @@ def ceof(lon, lat, data, nkp = 10, parallel = True): print('3: Solving the eigenvalue problem') lamda, loadings = delayed(la.eig)(c).compute() # lamda: eigenvalue, loadings: eigenvectors - l = lamda.conjugate().T; k = np.argsort(l) + l = lamda.conjugate().T + k = np.argsort(l) lamda, loadings = np.flip(l[k]), np.fliplr(loadings[:,k]) loadings = loadings[:,:nkp] # In case there were nan values in the orginal data, we need to perform the approach below: diff --git a/OceanLab/oa.py b/OceanLab/oa.py index 9853dcf..2360eb6 100644 --- a/OceanLab/oa.py +++ b/OceanLab/oa.py @@ -94,12 +94,12 @@ def vectoa(Xg,Yg,X,Y,U,V,corrlenx,corrleny,err,b=0): tc = np.array(tc) tc.shape = ppc[0].shape d2=((np.tile(xc,(n,1)).T-np.tile(x,(nv,1)))**2+(np.tile(yc,(n,1)).T-np.tile(y,(nv,1)))**2) - R=np.exp(-lambd*d2)+bmo; + R=np.exp(-lambd*d2)+bmo P=np.zeros((nv,2*n)) # streamfunction-velocity covariance - P[:,0:n]=np.sin(tc)*np.sqrt(d2)*R; - P[:,n:2*n]=-np.cos(tc)*np.sqrt(d2)*R; + P[:,0:n]=np.sin(tc)*np.sqrt(d2)*R + P[:,n:2*n]=-np.cos(tc)*np.sqrt(d2)*R PSI=np.dot(P,np.linalg.solve(A,uv)) # solvi the linear system PSI=PSI.reshape(nv2,nv1).T @@ -163,7 +163,7 @@ def scaloa(xc, yc, x, y, t=[], corrlenx=None,corrleny=None, err=None, zc=None): # Gauss-Markov to get the weights that minimize the variance (OI). tp = None ep = 1 - np.sum(C.T * np.linalg.solve(A, C.T), axis=0) / (1 - err) - if any(t)==True: ##### was t!=None: + if any(t): ##### was t!=None: t = np.reshape(t, (n, 1)) tp = np.dot(C, np.linalg.solve(A, t)) #if 0: # NOTE: `scaloa2.m` @@ -174,7 +174,7 @@ def scaloa(xc, yc, x, y, t=[], corrlenx=None,corrleny=None, err=None, zc=None): # tp = tp + mD * np.ones(tp.shape) return tp, ep - if any(t)==False: ##### was t==None: + if not any(t): ##### was t==None: print("Computing just the interpolation errors.") #Normalized mean error. Taking the squared root you can get the #interpolation error in percentage. diff --git a/OceanLab/utils.py b/OceanLab/utils.py index 8e793a7..6b8c35b 100644 --- a/OceanLab/utils.py +++ b/OceanLab/utils.py @@ -2,8 +2,7 @@ import scipy.signal as sg import xarray as xr -import dask -from dask.distributed import Client, progress +from dask.distributed import Client ##### User functions #=============================================================================