Trackeddy is an algorithm which identifies and tracks eddies assuming that their outer most contours is well fitted by an ellipse (A. Fernandes and S. Nascimento, (2006)). Additionally, the area of the eddy contour should be smaller than the first baroclinic Rossby Radius of Deformation (Mesoscale Scaling) defined as: \begin{equation} 2\pi L_r \end{equation} (Klocker, A., & Abernathey, R. (2014)).
To avoid the detection of meanders (jets), and due to the coherience of eddies, the eccentricity of the fitted ellipse eddies should be smaller than \begin{equation} \frac{b}{2a} \end{equation} which qualitatively corresponds to a semi-major axis two times larger than the semi-minor axis.
Optionally, an additional criterion is whe the 2D Gaussian fitting is implemented, this criterion identifies eddies only if the fitted 2D Gaussian has a fitness over 90%.
To install trackeddy please follow the instructions at Trackeddy ReadTheDocs.
Please be aware that this package still under development, therfore any issue can be reported at the Github Page. or please send an email to josue.martinezmoreno@anu.edu.au.
This code will be transfer to classes as soon as an stable version is released.
In [1]:
# Importing all libraries.
from pylab import *
from netCDF4 import Dataset
import os
import cmocean as cm
from trackeddy.tracking import *
from trackeddy.datastruct import *
from trackeddy.geometryfunc import *
from trackeddy.init import *
from trackeddy.physics import *
from trackeddy.plotfunc import *
import matplotlib.pyplot as plt, mpld3
#Allow interactive plots
%matplotlib inline
In [2]:
# Load Data
filepath = '../input/dt_global_allsat_phy_l4_20160901.nc'
In [3]:
# Open netcdf Dataset.
ncfile = Dataset(filepath)
# Load data into memory
sla = ncfile.variables['sla'][:]
lon = ncfile.variables['longitude'][:]
lat = ncfile.variables['latitude'][:]
In [4]:
# Define area of study
areamap = array([[0,len(lon)],[0,len(lat)]]) # Global option
The next code section contain the parameters used by track eddy to define and identify eddies. To know more about it, please refer to the Trackeddy Documentation. The following parameters correspond to the default values, except for the spatial filter.
In [5]:
# Time and spatial filter
filters = {'time':{'type':None,'t':None,'t0':None,'value':None},
'spatial':{'type':'moving','window':50,'mode':'uniform'}}
# Mesoscale scaling
checkarea={'mesoscale':2*np.pi}
# Eddy definition criteria
preferences={'ellipse':0.85,'eccentricity':0.85,'gaussian':0.8}
In [6]:
# Levels to be analysed and to extract positive eddies from anomaly
levels = {'max':sla[0,:,:].max(),'min':0.01,'step':0.03}
Make sure the dimention of the input matrix is 3D, which where the first index corresponds to time and the other two correspond to spatial dimentions. For example:
shape(sla)
>>> (1, 720, 1440)where 1 corresponds to the time step of the 1st of September 2016, and 720 grid cells in latitude and 1440 grid cells in longitude.
In [7]:
positive_eddies=analyseddyzt(sla,lon,lat,0,1,1,levels,preferences=preferences
,areamap=areamap,areaparms=checkarea,filters=filters
,maskopt='contour',diagnostics=False,pprint=True)
The output of analyseddyzt is a dictionary which contain all the parameters that define an eddy.
To convert it back to a field, its necesary to use the following function:
In [8]:
positive_eddy_field=reconstruct_syntetic(shape(sla),lon,lat,positive_eddies)
In [9]:
basemap_mplot(lon,lat,positive_eddy_field[0,:,:],title="Positive Eddies",projection='mbtfpq',
lat_0=-90,lon_0=-180,resolution='c',vmin=-0.80,vmax=0.80,
cmap=cm.cm.balance,xan=1,yan=1,figsize=(5,2.5),fontsize=12)
colorbar()
Out[9]:
In [10]:
# Levels to be analysed and to extract negative eddies from anomaly
levels = {'max':sla[0,:,:].min(),'min':-0.01,'step':-0.03}
In [11]:
negative_eddies=analyseddyzt(sla,lon,lat,0,1,1,levels,preferences=preferences
,areamap=areamap,areaparms=checkarea,filters=filters
,maskopt='contour',diagnostics=False,pprint=True)
In [12]:
negative_eddy_field=reconstruct_syntetic(shape(sla),lon,lat,negative_eddies)
In [13]:
basemap_mplot(lon,lat,negative_eddy_field[0,:,:],title="Negative Eddies",projection='mbtfpq',
lat_0=-90,lon_0=-180,resolution='c',vmin=-0.80,vmax=0.80,
cmap=cm.cm.balance,xan=1,yan=1,figsize=(5,2.5),fontsize=12)
colorbar()
Out[13]:
In [15]:
basemap_mplot(lon,lat,[sla[0,:,:],+negative_eddy_field[0,:,:]+positive_eddy_field[0,:,:]],
title=["SLA","Reconstructed SSHa"],projection='mbtfpq',
lat_0=-90,lon_0=-180,resolution='c',vmin=-0.80,vmax=0.80,
cmap=cm.cm.balance,xan=2,yan=1,figsize=(6,2.5),fontsize=12)
Out[15]:
In [18]:
basemap_mplot(lon,lat,sla[0,:,:]-(negative_eddy_field[0,:,:]+positive_eddy_field[0,:,:]),
title="SLA - Reconstructed SSHa",projection='mbtfpq',
lat_0=-90,lon_0=-180,resolution='c',vmin=-0.80,vmax=0.80,
cmap=cm.cm.balance,xan=1,yan=1,figsize=(6,2.5),fontsize=12)
colorbar()
Out[18]:
In [ ]: