News Column

A New Sea Surface Height-Based Code for Oceanic Mesoscale Eddy Tracking

May 1, 2014

McWilliams, James C


This paper presents a software tool that enables the identification and automated tracking of oceanic eddies observed with satellite altimetry in user-specified regions throughout the global ocean. As input, the code requires sequential maps of sea level anomalies such as those provided by Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO) data. Outputs take the form of (i) data files containing eddy properties, including position, radius, amplitude, and azimuthal (geostrophic) speed; and (ii) sequential image maps showing sea surface height maps with active eddy centers and tracks overlaid. The results given are from a demonstration in the Canary Basin region of the northeast Atlantic and are comparable with a published global eddy track database. Some discrepancies between the two datasets include eddy radius magnitude, and the distributions of eddy births and deaths. The discrepancies may be related to differences in the eddy identification methods, and also possibly to differences in the smoothing of the sea surface height maps. The code is written in Python and is made freely available under a GNU license ( emason/py-eddy-tracker/).

1. Introduction

Satellite altimetry has revealed the ubiquity of mesoscale eddies in the global ocean (e.g., Stammer 1997, 1998). Eddies range greatly in shape and size, are often asymmetric, and can have highly variable translational and rotational velocities (McWilliams 2008; Chelton et al. 2011b, hereafter CSS11; Early et al. 2011). Interest in mesoscale eddies arises from their role in the dynamics of the large-scale oceanic circulation; eddies are efficient carriers of mass and its physical, chemical, and biological properties, such that their presence modulates fluxes of heat and momentum and the dynamics of marine ecosystems (Chelton et al. 2011a; Gruber et al. 2011; Stramma et al. 2013).

Recent years have seen the emergence of several automated oceanic eddy tracking algorithms that contribute to knowledge of eddy properties and their variability. The techniques comprise threemainmethods: geometric (e.g., Chaigneau et al. 2008; Nencioli et al. 2010; CSS11); Okubo-Weiss (e.g., Isern-Fontanet et al. 2003; Morrow et al. 2004; Chelton et al. 2007; Ubelmann and Fu 2011); and wavelet (e.g., Doglioli et al. 2007; Rubio et al. 2009); and a comparative analysis of these approaches has been made by Souza et al. (2011). Novel techniques falling outside these methods are a hybrid geometric Okubo- Weiss approach (Halo et al. 2014), and an objective approach based on geodesic transport theory (Beron- Vera et al. 2013).

Two decades' worth of merged global sea surface height (SSH) from multiple satelliteborne altimeters are presently available from Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO) data, providing improved resolution of mesoscale features and their variability (Pascual et al. 2006). CSS11 applied an SSH anomaly-based eddy tracking algorithm to these data for the period 1992 to (at time of writing) 2012, and they make available a periodically updated database of the tracks and associated properties on their website.

In this article we evaluate and release a new SSHbased eddy tracking code that demonstrates comparable results to the CSS11 data. By providing code rather than a database, users may (i) track eddies using the most recent AVISO data, and (ii) adapt the code to specific regions and/or different data sources, such as oceanic numerical model outputs. The code is written in Python (e.g., Hunter 2007; Oliphant 2007); details of the algorithm are presented in section 2. The eddy tracker (pyeddy- tracker) is applied in section 3 to the region of the Canary Eddy Corridor (CEC; Sangra et al. 2009) offthe northwestern African coast. Comparisons with two other regions, the South Atlantic Ocean (SAO; Souza et al. 2011) and southeast Pacific (SEP; Chelton et al. 2011a), are summarized in tabular form. A discussion of some of the differences with CSS11 is provided in section 4, followed by brief conclusions in section 5.

2. Methodology

The py-eddy-tracker SSH-based approach is loosely based on procedures described by CSS11, Kurian et al. (2011), and Penven et al. (2005). Input sea level anomalies (SLA) are the same gridded weekly delayed-time global reference data from AVISO (e.g., AVISO 2013) as used by CSS11.

a. Eddy identification

Individual SLA fields are spatially high-pass filtered by removing a smoothed field obtained from a Gaussian filter with a zonal (meridional) major (minor) radius of 108 (58). An example filtered global SLA map on 28 August 1996 (not shown) makes good visual comparison with the high-pass filtered SLA field presented by CSS11 (their Fig. 1, bottom), although we are unable to make a quantitative comparison. The code works with interpolated contours of smoothed SLA; this is analogous to Kurian et al. (2011) but distinct from CSS11, who use raw SLA pixels. The impacts of any differences between ours and CSS11's results that may arise from the respective smoothing and contour handling procedures (henceforth inherent differences) are difficult to assess.

The eddy identification process requires specification of a tracking domain; the CEC tracking domain is shown in Fig. 1a. From a subset of the global SLA corresponding to the tracking domain, SLA contours are computed at 1-cm intervals for levels2100 to 100 cm. To identify cyclones (CC) [anticyclones (AC)], the contours are searched from 100 (2100)cm downward (upward).At each SLA interval, closed contours (CC) are sequentially identified and analyzed; for selection as a potential eddy, a CC must meet a series of criteria relating to its shape and interior characteristics:

(i) Pass a shape test with error%55%, where the error is defined as the ratio between the areal sum of CC deviations from its fitted circle and the area of that circle (see Fig. 1b; Kurian et al. 2011).

(ii) Contain a pixel count, I, satisfying Imin # I # Imax, where Imin 5 8 and Imax 5 1000.

(iii) Contain only pixels with SLA values above (below) the current SLA interval value for anticyclonic (cyclonic) eddies.

(iv) Contain no more than one local SLA maximum (minimum) for an anticyclone (cyclone). This constraint differs from CSS11, who permit multiple local maxima/minima; see section 4 for a comment on the potential implications of this difference.

(v) Have amplitude (A) that satisfies 1 # A # 150 cm.

The above criteria, barring the shape test and local minima/maxima threshold, are identical to CSS11 (see their section B3).

When a CC passes the above-mentioned tests, it is identified as an eddy and, henceforth, is referred to as the effective perimeter of the eddy (Ceff). Following CSS11, an associated effective radius (Leff) is defined as the radius of a circle with the same area as the region enclosed by Ceff. These features are illustrated in red in Fig. 1b; the centroid of Ceffis denoted by Peffand the same-area circle by 1Peff. Next, a speed-based eddy radius (Lspd) is found, defined as the radius of the circle with the same area as the region within the CC of SLA with maximum average geostrophic speed (U, the rotational speed of the eddy). Eddy radius Lspd is estimated by iterating from Ceffinward over all CCs whose pixel count satisfies Imin # I # Imax; Lspd and its associated contour, Cspd, are shown in green in Fig. 1b. The tracking centroid (Ptrk) in black in Fig. 1b corresponds to the contour of the last iteration (Ctrk, where I/Imin).

Finally, SLA pixels corresponding to the eddy are masked, making the region unavailable for further eddy identification. Figure 1a shows the positions of eddies identified following the above-mentioned procedures within the CEC tracking domain on 28 August 1996.

b. Eddy tracking

Eddy tracking is accomplished using positions Ptrk, in contrast to CSS11, who use Peff(Fig. 1b). We choose Ptrk rather than Peff(or Pspd) because, while in most instances the three coordinates are closely located, sometimes they are not. In these cases, tracking with Peffor Pspd occasionally resulted in spuriously dropped tracks; this behavior was improved by using Ptrk.

Cyclones and anticyclones are treated separately. Separation distances between all identified Ptrk positions at time steps k and k 1 1 are computed. Assignation of k 1 1 candidates to k active eddies is decided using the ellipse method of CSS11. If multiple k 1 1 eddies fall within the ellipse, then the eddy is assigned according to the minimum of a set of dimensionless similarity parameters, S (Penven et al. 2005), that are computed for each k 1 1 candidate:

(ProQuest: ... denotes formula omitted.)

where Dd is the separation distance between eddies at times k and k 1 1, Da is the variation of eddy area (based on Lspd), and DA is the variation of amplitude. Characteristic values for eddy separation distance, area, and amplitude are given by, respectively, d0, a0, and A0. We use the same values, namely, d05 25km (based on the AVISO time scale of one week), a0 5 p602km2, and A0 5 2 cm, in all the experiments presented in section 3.

Unused k 1 1 eddies are considered to be new eddies. The coordinates Ptrk and parameters Leff, Lspd, A, and U of eddies with lifetimes greater than a threshold minimum number of days are stored; examples of 28-day minimum eddy tracks (the same threshold is applied to the CSS11 data in section 3) are shown in Fig. 1a. The eddy tracking process continues by iteration over the temporal tracking domain.

3. Applications

To demonstrate the performance of py-eddy-tracker, we apply the algorithm to three distinct tracking domains: the CEC, SAO, and SEP. The runs cover the same 14 October 1992-4 April 2012 time period. Annual climatologies of eddy density (all observations, and also birth and death locations) and polarity are computed, as well as histograms of eddy lifetimes, radii, and amplitudes. Identical figures prepared from the CSS11 dataset are included for comparison. Figures are shown only for the CEC analysis domain, outlined in gray in Fig. 1a. Statistics from all three analysis domains are shown in tabular form.

Figure 2 compares spatial patterns of CEC anticyclone and cyclone density (eddy numbers per 18 square per year) and polarity from py-eddy-tracker and CSS11. Overall, good qualitative and quantitative agreement is found. Both methods identify elevated mesoscale activity associated with the Canary Islands (CI; Barton et al. 2004). The 95% confidence intervals from Student's t tests along 268N are approximately two eddy counts wide (not shown), indicating that the higher counts near the CI are real. The polarity maps also display closely matching patterns. Anticyclone dominance is strong along the upwelling and is also significant along a wide zonal band at ;258N. Pearson correlation coefficients (r) between the respective datasets are positive and large for both eddy signs (r . 0.8), although slightly smaller for polarity (Table 1).

Figure 3 presents eddy statistics from the two datasets. The annual mean numbers of eddies with lifetimes between 4 and ;25 weeks are similar (Fig. 3a); py-eddytracker eddies are slightly more numerous at 20 weeks and less. More than;30 weeks, the eddy count is mostly fewer than 10 eddies per year. The py-eddy-tracker records one or two cyclones with lifetimes exceeding 60 weeks; CSS11, however, find cyclones with ages greater than 100 weeks. Cyclones are generally dominant in both datasets (Fig. 3b).

Histograms of eddy radii (Lspd) from CSS11 and pyeddy- tracker in Fig. 3c show some differences. CSS11 find a wider distribution of radii, with a peak of ;40 eddies per year at Lspd ' 75 km; the py-eddy-tracker data have a tendency for more and smaller eddies, with a peak of ;55 eddies per year at Lspd ' 65 km. Explanations for these differences are explored in section 4. Both datasets show anticyclonic (cyclonic) dominance at smaller (larger) scales (Fig. 3d). Figures 3e,f reveal close agreement in the distribution and polarity of eddy amplitudes between the two datasets. There is also good agreement in eddy nonlinearity1 although py-eddytracker tends to predict slightly more nonlinear eddies (Figs. 3g,h).

Scatterplots of eddy radius versus amplitude for cyclones and anticyclones from the two datasets are shown in Figs. 3i,j. There is good general agreement between the datasets in the distributions for both eddy signs; however, a marked difference in the magnitude of the radii is visible at smaller amplitudes.

Annual mean CEC eddy birth and death location densities from py-eddy-tracker and CSS11 are compared in Fig. 4. The density distributions are generally comparable, although py-eddy-tracker consistently records more eddy births and deaths in the open ocean than CSS11. Both methods confirm the role of the CI in eddy generation, with births of both signs occurring in the lee of the archipelago. Pearson's r in Table 1 for pyeddy- tracker and CSS11 eddy births and deaths are positive but somewhat smaller than for eddy density, a reflection of the open ocean differences between the methods.

Figure 5 shows comparative maps of mean anticyclonic and cyclonic eddy radius and amplitude. The radius distributions from py-eddy-tracker and CSS11 are similar, with the largest scales in the open ocean to the south and west (Figs. 5a,b,e,f). However, the CSS11 radii are characteristically larger, as previously seen in Figs. 3c,i,j. Close agreement is evident between the datasets in the distribution and magnitude of mean eddy amplitudes (Figs. 5c,d,g,h).

The good comparative performance between pyeddy- tracker and CSS11 in the CEC analysis domain is replicated over the SAO and SEP regions, as evidenced by the respective Pearson's r in Table 1.

4. Discussion

The main discrepancies between the CEC py-eddytracker results and the eddy track data from CSS11 in section 3 concern Lspd eddy radius magnitude and the offshore distributions of eddy births and deaths. The pyeddy- tracker estimates generally smaller radii and more frequent instances of eddy generation and termination in the open ocean. Similar differences are observed for the SAO and SEP analysis domains.

The py-eddy-tracker eddy identification criteria are arguably stricter than CSS11's owing to the use of a shape test and the limit of a single local minima/maxima within Ceff(section 2a). As a result, the probability of identifying large irregularly shaped feature as eddies is greater in the CSS11 method (cf. the region Lspd * 125 km, A & 4.5 cm in Figs. 3i,j). When evaluating such features, py-eddytracker iterates farther toward the center before (if symmetry increases and local minima/maxima do not exceed one) eddies can be identified; those py-eddytracker eddies will clearly have smaller Leffthan the equivalent CSS11 eddies, and hence smaller Lspd. (Albeit an approximate relationship, CSS11 report Lspd ' 0.7 Leff; this value is confirmed for the three analysis domains herein.) A further factor that may influence the Lspd bias are the aforementioned inherent differences in SLA treatment (section 2a).

Concerning eddy birth and death densities, we find that in open ocean regions py-eddy-tracker predicts greater frequencies of both events than CSS11 (Fig. 4; also Fig. 3a). The eddy identification disparitymentioned above may again provide an explanation for the differences: having a stricter eddy identification threshold, py-eddy-tracker is more likely to terminate a track than CSS11 (hence more deaths); more terminations mean in turn a lower likelihood that genuine new eddies are spuriously linked to existing eddies (hence more births). This delicate balance highlights the important role of the eddy identification procedure. Unfortunately, not having access to CSS11 individual SLA fields, we are unable to directly compare the eddy identification performance of the two algorithms.

5. Concluding remarks

A new oceanic eddy tracking code written in the open source Python language has been introduced. Results from applying the py-eddy-tracker code to SLA data compare favorably with the eddy property database of CSS11. Both methods produce broadly consistent results in terms of eddy locations, scales, and amplitudes within three distinct analysis domains. Discrepancies between the methods relating to eddy scales, and distributions of eddy births and deaths are identified and discussed.

The py-eddy-tracker code works on global SLA over user-defined regions up to basin scale. At global scales the present code slows down significantly. Improving the iteration speed and the addition of oceanic numerical model eddy tracking capability (surface and subsurface) are important next steps.

There is arguably no uniquely correct eddy detection method, but several methods, including that of CSS11 and the code presented herein, produce useful enough results for most purposes of eddy counting. As eddy tracking techniques evolve and data resolution and coverage improve, comprehensive intercomparisons between different oceanic eddy trackers will reveal the extent to which results from the various codes are robust and where there may be large method-related uncertainties (e.g., Souza et al. 2011; Neu et al. 2013).

Acknowledgments. Evan Mason is supported by a Spanish government JAE Doc grant (CSIC), cofinanced by FSE. This work has been partially funded by the project MyOcean-2 EU FP7. We are grateful to Jaison Kurian for sharing his MATLAB functions for the shape test. We thank two anonymous reviewers, whose contributions helped to improve the manuscript; and also Dudley Chelton and Michael Schlax for their supportive comments on early versions of this work. The altimeter products were produced by SSALTO/ DUACS and distributed by AVISO, with support from CNES (

1U/c. 1 , here c is the translational speed of the eddy (CSS11).


AVISO, 2013: SSALTO/DUACS user handbook: (M)SLA and (M)ADT near-real time and delayed time products. SALP-MUP- EA-21065-CLS, 2nd ed. AVISO, 70 pp. [Available online at aviso-user-handbooks.html.].

Barton, E. D., J. Aristegui, P. Tett, and E. Navarro-Perez, 2004: Variability in the Canary Islands area of filament-eddy exchanges. Prog. Oceanogr., 62, 71-94, doi:10.1016/ j.pocean.2004.07.003.

Beron-Vera, F. J., Y. Wang, M. J. Olascoaga, G. J. Goni, and G. Haller, 2013: Objective detection of oceanic eddies and the Agulhas leakage. J. Phys. Oceanogr., 43, 1426-1438, doi:10.1175/JPO-D-12-0171.1.

Chaigneau, A., A. Gizolme, and C. Grados, 2008: Mesoscale eddies offPeru in altimeter records: Identification algorithms and eddy spatio-temporal patterns. Prog. Oceanogr., 79, 106-119, doi:10.1016/j.pocean.2008.10.013.

Chelton, D. B., M. G. Schlax, R. M. Samelson, and R. A. de Szoeke, 2007: Global observations of large oceanic eddies. Geophys. Res. Lett., 34, L15606, doi:10.1029/2007GL030812.

-, P. Gaube, M. A. Schlax, J. A. Early, and R. M. Samelson, 2011a: The influence of nonlinear mesoscale eddies on nearsurface oceanic chlorophyll. Science, 334, 6054, 328-332, doi:10.1126/science.1208897.

-, M. A. Schlax, and R. M. Samelson, 2011b: Global observations of nonlinear mesoscale eddies. Prog. Oceanogr., 91, 167- 216, doi:10.1016/j.pocean.2011.01.002.

Doglioli, A. M., B. Blanke, S. Speich, and G. Lapeyre, 2007: Tracking coherent structures in a regional ocean model with wavelet analysis: Application to cape basin eddies. J. Geophys. Res., 112, C05043, doi:10.1029/2006JC003952.

Early, J. A., R. M. Samelson, and D. B. Chelton, 2011: The evolution and propagation of quasigeostrophic ocean eddies. J. Phys. Oceanogr., 41, 1535-1555, doi:10.1175/2011JPO4601.1.

Gruber, N., Z. Lachkar, H. Frenzel, P. Marchesiello, M. Meurounnich, J. C. McWilliams, T. Nagai, and G.-K. Plattner, 2011: Eddyinduced reduction of biological production in eastern boundary upwelling systems. Nat.Geosci., 4, 787-792, doi:10.1038/ngeo1273.

Halo, I., B. Backeburg, P. Penven, I. Ansorge, C. Reason, and J. E. Ullgren, 2014: Eddy properties in the Mozambique Channel: A comparison between observations and two numerical ocean circulation models. Deep-Sea Res. II, 100, 38-53, doi:10.1016/ j.dsr2.2013.10.015.

Hunter, J. D., 2007: Matplotlib: A 2D graphics environment. Comput. Sci. Eng., 9, 90-95, doi:10.1109/MCSE.2007.55.

Isern-Fontanet, J., E. Garcia-Ladona, and J. Font, 2003: Identification ofmarine eddies fromaltimetricmaps. J. Atmos. Oceanic Technol., 20, 772-778, doi:10.1175/1520-0426(2003)20,772: IOMEFA.2.0.CO;2.

Kurian, J., F. Colas, X. J. Capet, J. C. McWilliams, and D. B. Chelton, 2011: Eddy properties in the California Current System. J. Geophys. Res., 116, C08027, doi:10.1029/ 2010JC006895.

McWilliams, J. C., 2008: The nature and consequences of oceanic eddies. Ocean Modeling in an Eddying Regime, Geophys. Monogr., Vol. 177, Amer. Geophys. Union, 5-15, doi:10.1029/ 177GM03.

Morrow, R., F. Birol, D. Griffin, and J. Sudre, 2004: Divergent pathways of cyclonic and anti-cyclonic ocean eddies. Geophys. Res. Lett., 31, L24311, doi:10.1029/2004GL020974.

Nencioli, F., C. Dong, T. D. Dickey, L. Washburn, and J. C. McWilliams, 2010: A vector geometry-based eddy detection algorithm and its application to a high-resolution numerical model product and high-frequency radar surface velocities in the Southern California Bight. J. Atmos. Oceanic Technol., 27, 564-579, doi:10.1175/2009JTECHO725.1.

Neu, R., and Coauthors, 2013: IMILAST: A community effort to intercompare extratropical cyclone detection and tracking algorithms. Bull. Amer. Meteor. Soc., 94, 529-547, doi:10.1175/ BAMS-D-11-00154.1.

Oliphant, T. E., 2007: Python for scientific computing. Comput. Sci. Eng., 9, 10, doi:10.1109/MCSE.2007.58.

Pascual, A., Y. Faugere, G. Larnicol, and P.-Y. Le Traon, 2006: Improved description of the ocean mesoscale variability by combining four satellite altimeters. Geophys. Res. Lett., 33, L02611, doi:10.1029/2005GL024633.

Penven, P., V. Echevin, J. Pasapera, F. Colas, and J. Tam, 2005: Average circulation, seasonal cycle, and mesoscale dynamics of the Peru Current System: A modeling approach. J. Geophys. Res., 110, C10021, doi:10.1029/2005JC002945.

Rubio, A., B. Blanke, S. Speich, N. Grima, and C. Roy, 2009: Mesoscale eddy activity in the southern Benguela upwelling system from satellite altimetry and model data. Prog. Oceanogr., 83, 288-295, doi:10.1016/j.pocean.2009.07.029.

Sangra, P., and Coauthors, 2009: The Canary Eddy Corridor:Amajor pathway for long-lived eddies in the subtropical North Atlantic. Deep-Sea Res. I, 56, 2100-2114, doi:10.1016/j.dsr.2009.08.008.

Souza, J. M. A. C., C. de Boyer Montegut, and P.-Y. Le Traon, 2011: Comparison between three implementations of automatic identification algorithms for the quantification and characterization of mesoscale eddies in the South Atlantic Ocean. Ocean Sci., 7, 317-334, doi:10.5194/os-7-317-2011.

Stammer, D., 1997: Global characteristics of ocean variability estimated from regional TOPEX/POSEIDON altimeter measurements. J. Phys. Oceanogr., 27, 1743-1769, doi:10.1175/ 1520-0485(1997)027,1743:GCOOVE.2.0.CO;2.

-,1998: On eddy characteristics, eddy transports, and mean flow properties. J. Phys. Oceanogr., 28, 727-739, doi:10.1175/ 1520-0485(1998)028,0727:OECETA.2.0.CO;2.

Stramma, L., H. W. Bange, R. Czeschel, A. Lorenzo, and M. Frank, 2013: On the role of mesoscale eddies for the biological productivity and biogeochemistry in the eastern tropical Pacific Ocean offPeru. Biogeosci. Discuss., 10, 9179-9211, doi:10.5194/ bgd-10-9179-2013.

Ubelmann, C., and L.-L. Fu, 2011: Vorticity structures in the tropical Pacific from a numerical simulation. J. Phys. Oceanogr., 41, 1455-1464, doi:10.1175/2011JPO4507.1.


Instituto Mediterraneo de Estudios Avanzados, Consejo Superior de Investigaciones CientÝ'ficas, University of the

Balearic Islands, Esporles, Illes Balears, Spain


Institute of Geophysics and Planetary Physics, University of California, Los Angeles, Los Angeles, California

(Manuscript received 24 January 2014, in final form 7 February 2014)

Corresponding author address: E. Mason, IMEDEA, CSIC- UIB, C./Miquel Marques 21, Esporles 07190, Islas Baleares, Spain.


For more stories covering the world of technology, please see HispanicBusiness' Tech Channel

Source: Journal of Atmospheric and Oceanic Technology

Story Tools Facebook Linkedin Twitter RSS Feed Email Alerts & Newsletters