Scalable Feature Extraction and Tracking (SCAFET): a general framework for feature extraction from large climate data sets

This study describes a generalized computational mathematical framework, Scalable Feature Extraction and Tracking (SCAFET), that extracts and tracks features from large climate data sets. SCAFET utilizes novel shape-based metrics that can identify and compare features from different mean states, data sets, and between distinct regions. Features of interest such as atmospheric rivers, tropical and extratropical cyclones, and jet streams are extracted by segmenting the data based on a scale-independent bounded variable called the shape index (SI). The SI gives a quantitative measurement of the local geometric shape of the field with respect to its surroundings. Compared to other widely used frameworks in feature detection, SCAFET does not use a posteriori assumptions about the climate model or mean state to extract features of interest and levelize the comparison between different models and scenarios. To demonstrate the capabilities of the method, we illustrate the detection of atmospheric rivers, tropical and extratropical cyclones, sea surface temperature fronts, and jet streams. Cyclones and atmospheric rivers are extracted to show how the algorithm identifies and tracks both the nodes and areas from climate data sets. The extraction of sea surface temperature fronts exemplifies how SCAFET effectively handles curvilinear grids. Last, jet streams are extracted to demonstrate how the algorithm can also detect three-dimensional features. As a generalized framework, SCAFET can be implemented to extract and track many weather and climate features across scales, grids, and dimensions.

1 Introduction

The volume of climate data is growing exponentially, owing to rapid expansions in both observational capabilities and computational power, which are driven in particular by the precision and insights offered by higher-resolution models (Overpeck et al., 2011; Balaji et al., 2018). Frontier research like global cloud-resolving and large-ensemble simulations leads not only to increased volume but also to inflated velocity, variety, veracity, and value (5Vs) (Marr, 2015; Guo, 2017; van Genderen et al., 2019) of climate data. This makes the detection and comparative analysis of important atmospheric and oceanic features, such as atmospheric rivers (ARs), tropical and extratropical cyclones, sea surface temperature fronts (SSTFs), and jet streams, an onerous task. Although these climate phenomena influence regional and global weather and climate with immense societal, economic, and ecological impacts, the volume of data representing these events and features is a small percentage of the whole simulation. Feature extraction considerably reduces the volume of data that needs to be stored, thus improving computational efficiency in analyzing these features (Yang et al., 2016). Moreover, the mean, variability, and characteristics of features can be compared to observational data sets as a measure of bias within model simulations, thereby improving our understanding of the causal differences between observations and models (Sellars et al., 2013). Thus, efficient and reliable feature extraction is vital to climate data processing, analysis, and model development.

Despite the importance of feature extraction in climate data analysis and model development, there is little consensus on standard best practices for feature extraction. The simplest method for extracting a feature is to use a physical threshold or its derivative for some climate variable (SST, precipitation, wind speed, humidity, etc.), or a combination thereof, to identify ARs, fronts, jet streams, or tropical and extratropical cyclones (Bengtsson et al., 1982, 1995; Vitart et al., 1997; Hewson, 1998; Koch et al., 2006; Strong and Davis, 2007; Rutz et al., 2014; Guan and Waliser, 2015). The limitations and discrepancies in these methods arise from the somewhat arbitrary choice of physical thresholds in relation to the spatiotemporal distributions of the climate variables. In other words, many studies choose a physical threshold that is not theoretically defined but rather a function of the location, time span, and data set used. Validation can then unfortunately come down to the intelligent but subjective human eye or, in other words, tuning an absolute or relative threshold until it appears to have captured all the features of interest, while leaving out the background noise (Zarzycki and Ullrich, 2017; Vishnu et al., 2020).

Choosing an absolute threshold from climate variables for feature extraction that applies to different climate models and spans multiple mean states and model scenarios is not straightforward. Thresholds are often applied to climate variables or derivatives in which the features are most visible, such as relative vorticity (RV) and sea level pressure anomalies for tropical cyclones (e.g., Vitart et al., 1997), integrated water vapor transport (IVT) for ARs (e.g., Guan and Waliser, 2015), or the first derivative of sea surface temperature (SST) for SST fronts (Castelao et al., 2006). Thresholds are often either empirically derived from observational studies or calculated from a model-specific distribution; though, even within the same data set, a particular choice of threshold may be suitable for one region but not for another, given varying regional characteristics and topography. In the case where the feature extraction threshold is an a posteriori assumption of the data set used, one must pre-process large, representative data sets just to calculate reasonable thresholds. While some detection methods have done well to streamline their algorithms to reduce total runtime, the process of posterior threshold calculation for higher-resolution and large-ensemble data sets inherently becomes increasingly less efficient, highlighting the need to develop feature extraction methods that do not use posterior assumptions.

Aside from the sensitivity of feature detection to inter-model and inter-simulation differences, feature detection is further complicated when trying to detect and compare features between present and future climate change scenarios, as the underlying spatiotemporal climate variable distributions change under global warming. Feature detection must be reconsidered when applied to variables with significant and/or non-linear changes in their means and extremes in response to external forcings such as doubling or quadrupling carbon dioxide concentrations. It should be emphasized that applying different arbitrary thresholds can and does lead to contradictory conclusions regarding the response of these features to greenhouse gas warming (Horn et al., 2014; Zhao, 2020; O'Brien et al., 2022; Nellikkattil et al., 2023). To counter these uncertainties, methods based on topology, machine learning, ridge extraction, edge detection, and various other image-processing techniques have been proposed over the years (Dixon and Wiener, 1993; Post et al., 2003; Molnos et al., 2017; Biard and Kunkel, 2019; Xu et al., 2020). While these methods offer an alternative for the extraction of features in data sets spanning different mean states, many of these methods were developed for detecting specific rather than general features.

The need for a general framework for extracting and tracking features from large climate data sets has been raised in various climate science communities for the last several decades. In a pioneer study, Hodges (1994) developed a general framework for extracting and tracking features from meteorological data sets in the following three steps: segmentation, filtering, and tracking. In the segmentation step, the field is split into distinct regions by applying a threshold and defining each of the connected regions as an object. Segmented regions are then filtered based on the characteristics of each object, and feature nodes are defined for the remaining objects. Finally, the feature nodes are tracked over time to produce the final output for further analysis. This framework was further developed for cyclones, storm tracks, convective systems, ocean eddies, monsoon depressions, etc. (Hodges, 1995; Hogg et al., 2005; Hodges et al., 2011; Burston et al., 2014; Hurley and Boos, 2014; Pinheiro et al., 2016; Priestley et al., 2020; Torres-Alavez et al., 2021; Karmakar et al., 2021). However, it is limited to the detection of points of local maxima in two-dimensional scalar fields, which do not always fully characterize various features.

In 2012, a team from the Lawrence Berkeley National Laboratory developed the Toolkit for Extreme Climate Analysis (TECA), integrating pre-existing, physical threshold-dependent detection methods and algorithms into a comprehensive software package that was parallelized to make the algorithms more suitable for large data sets (Prabhat et al., 2012). In a more recent effort, a team led by Paul Ullrich at the University of California, Davis, created TempestExtremes (Ullrich and Zarzycki, 2017; Ullrich et al., 2021), which is another computationally efficient algorithm package that uses C$++$ and several core functions to detect a variety of features. These functions are being actively developed for extraction, characterization, and uncertainty quantification of weather extremes. Both TECA and TempestExtremes have been widely implemented by the climate community and have been monumental in advancing scientific understanding of meso- and synoptic-scale processes and their connections to long-term climate variability.

In this study, we present a novel method called Scalable Feature Extraction and Tracking (SCAFET), which serves as a versatile and general framework for detecting and tracking features of various shapes and intensities across scales, grid types, and dimensions. Simply put, SCAFET uses the curvature measurements of a given scalar field to identify distinct emergent shapes corresponding to features of interest. The local shape calculation is finite, bounded, and scale-independent, and it can be tuned, depending on the specified feature of interest. Unlike traditional methods that rely on physical thresholds often derived from data-specific, posterior conditions, this method relies on shape-based thresholds. As such, it separates the feature detection process from inter- and intra-model variation, making it less sensitive to these differences. Furthermore, this approach allows for the complete parallelization of feature extraction along the time dimension, since the detection operates independently of time. Time-independent feature extraction offers two key advantages. First, it has the potential to boost computational efficiency by enabling data pre-processing such as smoothing to occur in parallel, rather than requiring a single pre-processing step before feature extraction. Second, it holds the promise of being developed and implemented for real-time feature extraction during critical events like hurricanes and tornadoes. Importantly, the code for this framework is fully open-source and written in Python in an easy-to-use package so that even individuals with beginner-level Python skills can readily implement the algorithm (see https://github.com/nbarjun/SCAFET/blob/master/scafet_demo.ipynb, last access: 17 October 2023, for a simple working example).

The novelty of SCAFET compared to pre-existing methods lies in feature detection that does not use a posteriori assumptions and is based on the overall “shape” of a climate variable field, rather than arbitrary thresholding of that field or derivative. The core methodology for the detection of any feature is the same and can be tuned using just two variables: one for the spatial scale and the other for the shape of the features one is looking for. For example, between the two variables, one can tune the difference between a long filament-shaped atmospheric river and a shorter, round-shaped cyclone. The algorithm applies to both rectilinear and curvilinear grids and can also be extended to detect three-dimensional (3D) features. Even in the context of recent advancements in feature extraction such as TempestExtremes and TECA, SCAFET is a comprehensive, efficient, and easily implementable framework that aims to upgrade the feature extraction process with a novel shape-based approach that does not rely on iterative posterior conditions and could prove to be a robust method for detecting a diverse set of features under different mean climate states. Further discussion on the differences between SCAFET and other detection algorithms can be found in Sect. A.

The paper is organized as follows, Sect. 2 introduces the fundamentals of SCAFET and how it is implemented in a two-dimensional (2D) field. Section 3 presents three specific use cases of SCAFET, demonstrating its capabilities in detecting various climate features across different grid types. Extraction of 3D features using jet steams as an example will be discussed in Sect. 4. Though the application of SCAFET is not limited to the features described here, this study focuses on atmospheric rivers, cyclones, SST fronts, and jet streams, as these examples cover a broad range of phenomena, providing users with insights on how to adapt SCAFET to their specific use cases and requirements.

2 Description of Scalable Feature Extraction and Tracking

SCAFET adopts the same three-step approach as outlined by Hodges (1994): segmentation (yellow boxes in Fig. 1), filtering (orange boxes in Fig. 1), and tracking (green boxes in Fig. 1). However, before commencing these steps, SCAFET requires initialization with essential information describing the data sets and the specific feature to be extracted (as indicated by blue boxes in Fig. 1). The key inputs for this initialization include the following:

Figure 1 Overall schematic of SCAFET workflow and components. Inputs to the algorithm are depicted in blue, while the algorithm's outputs are shown in pink boxes. Processes related to the segmentation step are highlighted in yellow boxes, whereas the orange boxes represent the filtering processes. The tracking step is denoted by green boxes. Arrows on the periphery of the boxes illustrate the flow of the algorithm. Each section is elaborated upon in detail within the text.

• Primary field (ϕ_p). This is a gridded data set in which the target feature is most easily distinguishable. For instance, cyclones are readily identified using the RV field, ARs emerge from IVT, and SSTFs are distinguished using the SST gradient. Optionally, one or more secondary fields can be used to further constrain the detected features.

• Grid properties. Information on the primary field's grid, including grid cell area/volume, grid distance, and coastlines, is required for calculating derivatives of the basic field and identifying locations of landfall.

• Feature properties. The algorithm requires information on the properties of the target feature. This includes estimated spatial scale, shape, eccentricity (for 2D features only), minimum length, minimum area, minimum volume (for 3D feature only), minimum duration, and maximum distance per time step.

In the SCAFET scheme, segmentation, filtering, and tracking are developed and coded as separate Python libraries. This design allows users to substitute any of these components with their own methods, while still being able to execute the algorithm. Once all three steps have been executed, the algorithm yields two outputs: one provides information about the properties of the detected objects and the other produces a labeled mask highlighting the feature of interest on the input grid (pink boxes in Fig. 1).

2.1 Segmentation

The core operation for the feature extraction involves categorizing points within a scalar field into one of five shapes. This categorization is achieved using curvature measurements obtained from the eigenvalues of the Hessian of the basic field. These five selected shapes (see Fig. 2) are an abridged version of the shapes described in previous studies (Koenderink and van Doorn, 1992). Depending on the specific feature of interest, one or more shapes are extracted from the primary field. The segmentation process starts with scale–space selection of the field to remove smaller scales of variability that are background noise compared to the feature of interest. Last, the algorithm calculates SI to estimate the local geometric shape at each point.

2.1.1 Scale–space selection

Scale–space selection is a widely used technique in image processing, signal processing, and computer vision (Lindeberg, 2014). In our current study, scale–space selection involves applying a Gaussian smoothing kernel to suppress variability smaller than the chosen smoothing scale (σ) (see https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.smooth_gaussian.html, last access: 20 December 2023, for the implementation of Gaussian smoothing). Mathematically, scale–space selection is performed by convolving the primary field (ϕ_p) with a Gaussian function, which is expressed as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\varphi }}_{\mathrm{s}}(x,y,\mathrm{\dots })={\mathit{\varphi }}_{\mathrm{p}}(x,y,\mathrm{\dots })\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }\mathit{\sigma }}}{e}^{-(x^{\mathrm{2}}+y^{\mathrm{2}}+\mathrm{\dots })/\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}.\end{array}$$

In the context of the meso–synoptic-scale processes examined in this study, the scale–space selection filters out smaller microscale features to isolate features like cyclonic vortexes or atmospheric rivers. Notably, this function can be adjusted to the spatial scale of interest and could also be used to filter out synoptic-scale features in isolating micro- and mesoscale processes. In climate data sets, grid spacing is not always uniform. To account for that, we adapt the above equation to be “grid-aware”. The input for the smoothing scale is provided in kilometers, and based on this input, we calculate the value of σ, while considering the grid size. Notably, the value of σ remains constant when smoothing is applied along each longitude, but it varies along each circle of latitude. For future studies, researchers may explore other, more advanced scale–space selection methods to further refine their analyses.

2.1.2 Local shape extraction

The local geometric shape of the field, ϕ_s is calculated as a function of the eigenvalues (k₁ and k₂) of the Hessian of the magnitude of the field ($\left|{\mathit{\varphi }}_{\mathrm{s}}\right|$), where the Hessian is given by

$$\begin{array}{}\text{(2)}& \mathcal{H}\left(\left|{\mathit{\varphi }}_{\mathrm{s}}\right|\right)=\left[\begin{array}{cc}\frac{{\partial }^{\mathrm{2}}\left|{\mathit{\varphi }}_{\mathrm{s}}\right|}{\partial x^{\mathrm{2}}}& \frac{{\partial }^{\mathrm{2}}\left|{\mathit{\varphi }}_{\mathrm{s}}\right|}{\partial x\partial y}\\ \frac{{\partial }^{\mathrm{2}}\left|{\mathit{\varphi }}_{\mathrm{s}}\right|}{\partial y\partial x}& \frac{{\partial }^{\mathrm{2}}\left|{\mathit{\varphi }}_{\mathrm{s}}\right|}{\partial y^{\mathrm{2}}}.\end{array}\right]\end{array}$$

Figure 2 Selected shapes used in this study and the values of the shape index associated with each of them. The x and y axes are a set of general axes, while $Z(X,Y)=\sin \left(\mathrm{2}X\right)+\cos \left(\mathrm{2}Y\right)$. Regions within Z(X,Y) satisfying conditions for different shapes are isolated to show the geometry associated with them.

In the context of simple differential geometry, we can determine whether a point is a local maximum or a local minimum based on the eigenvalues k₁ and k₂. Specifically, if $k_{\mathrm{2}}\le k_{\mathrm{1}}<\mathrm{0}$, then the point under consideration is a local maximum, whereas if $k_{\mathrm{1}}\ge k_{\mathrm{2}}>\mathrm{0}$, then point is a local minimum. The criterion is primarily applicable to nodal features such as tropical cyclones or monsoon depressions. To expand our ability to identify a range of features, we use the shape index (SI) (Koenderink and van Doorn, 1992), a quantitative measure of the local shape of the field defined as

$$\begin{array}{}\text{(3)}& \text{SI}(k_{\mathrm{1}},k_{\mathrm{2}})={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}{\tan }^{-\mathrm{1}}\left[{\displaystyle \frac{k_{\mathrm{2}}+k_{\mathrm{1}}}{k_{\mathrm{2}}-k_{\mathrm{1}}}}\right],\end{array}$$

where k₁ and k₂ are the two eigenvalues, satisfying k₁≥k₂, for the Hessian matrix. It is important to clarify that in the original work by Koenderink and van Doorn (1992), the principal curvatures, not the eigenvalues of the field, are utilized to calculate the SI. However, the disparity between the SI calculated using principal curvatures and the SI derived from eigenvalues is exceedingly minimal in climate data analysis. The SI is used to categorize the primary field into distinct shapes (see Fig. 2). The choice of SI values is contingent upon the specific type of feature to be extracted. For example, we select caps and domes when extracting features such as atmospheric depressions or cyclones, whereas ridges, caps, and domes are chosen when targeting features like ARs and fronts.

SI is designed to be a bounded value (range −1 to 1) independent of the magnitude of the field (Fig. 3). In simple terms, SI provides a continuous and quantitative measurement of the geometric shape of the field with respect to its immediate background field. This concept is similar to how a climate scientist's trained eye identifies features based on differences in color or value contrast, though the SI is arguably a more objective and precise measure of geometric shape. These characteristics make the SI particularly well-suited for feature extraction from data sets with varying mean states, which is in contrast to traditional physical threshold-based methods. In addition to the two eigenvalues, the shape extraction provides us with corresponding eigenvectors. The eigenvector for k₁ points perpendicular to the local ridge direction, while that of k₂ is parallel to it. This allows us to impose further constraints, such as the coherence of transport or flow with respect to the local ridge when ϕ_s is a vector field. This capability is aptly demonstrated in the context of AR detection, as discussed in Sect. 3.1.

Figure 3 Sensitivity of the shape index (SI) to eigenvalues k₁ and k₂. The x and y axes represent values of the two eigenvalues used for calculating the shape index, while the color indicates the value of the shape index. Shapes corresponding to the SI regimes are labeled. Shapes corresponding to SI regimes are labeled.

2.2 Filtering

Once the target features are extracted, properties like area, location, mean, minimum, and maximum values of different properties are calculated for each of the objects. A series of filtering is carried out to remove objects that do not satisfy certain conditions regarding (a) grid properties like area, length, and region masks, (b) primary field properties like magnitude and direction, and (c) constraints from the secondary field(s). The primary aim of the filtering process is to remove small, weak, or ephemeral objects.

2.3 Tracking

The properties extracted for each object include key positional details, such as its centroid, and endpoints, as well as the locations of maximum and minimum intensity of input field within each object. To follow objects through time, one of these positional attributes is tracked. In the present study, we employ a straightforward tracking method. For each object at time step n, we identify the closest object to it at time n+1. If this identified object is closer than a pre-defined radius r, we consider it to be the same object in motion. The radius r is defined in kilometers based on the maximum translation speed of the object and the temporal frequency of the input data. At this stage, it is possible to filter out short-lived features as needed. While this uncomplicated tracking approach may not be suitable for microscale processes, it can be adapted to incorporate greater complexity if necessary.

3 Application to 2D features

In this section, we showcase how SCAFET is employed to detect cyclonic vortices, ARs, and SSTFs from various climate data sets. These examples serve to illustrate the versatility of SCAFET as applied to different types of features and grids, though all the examples in this study follow the same general process shown in Fig. 1. Each subsection has a table of parameters detailing the properties of the desired feature. The properties include the feature's typical spatial scales, shape index (SI) regime, minimum length, minimum area, object eccentricity, and minimum duration of its track. To determine the quantitative values for these properties, we refer to a consensus among previous studies, which are cited within each section. A detailed examination of the sensitivity of these parameters in relation to the detected features, using AR detection as an example, can be found in Sect. S1 in the Supplement. In addition to the results discussed in the following sections, videos are also included in the Supplement for each of the features. The primary objective of this work is to demonstrate SCAFET's capability to detect a variety of features. Consequently, we present results for the long-term climatology of each of the features, enabling a comparison with other published detection algorithms.

3.1 Atmospheric rivers

According to the American Meteorology Society's glossary of meteorology, ARs are “long, narrow, and transient corridors of strong horizontal water vapor transport that are typically associated with a low-level jet stream ahead of the cold front of an extratropical cyclone” (Ralph et al., 2018). A substantial portion of the precipitation and water vapor transport in midlatitude regions is concentrated within ARs (Guan and Waliser, 2015). These atmospheric phenomena play a significant role in midlatitude hydrology, contributing to more than 50 % of the extreme precipitation and wind events in the region (Waliser and Guan, 2017; Nash et al., 2018). The ability to accurately detect, forecast, and project future ARs is of utmost importance for both extreme weather preparedness, as well as for water resource management in basins worldwide.

The ambiguity in AR projections and AR detection tools (ARDTs) stems from the lack of a clear quantitative definition of ARs in strength, length, narrowness, and other such parameters used in detection. In comparison with other criteria, the choice of threshold for AR strength has a significant effect on the inferences drawn between the detection schemes (Zhao, 2020; O'Brien et al., 2022; Nellikkattil et al., 2023). Many ARDTs determine this threshold empirically from the data set itself, which renders them sensitive to spatiotemporal variations and changes in mean-state conditions (Shields et al., 2018). SCAFET defines ARs as long (length >2000 km), narrow (eccentricity >0.75) regions of strong water vapor transport (SI >0.375) and significant precipitation (minimum AR precipitation >1 mm d⁻¹) (see Table 1 for complete details). The sensitivity of these parameters in AR detection to the characteristics of detected ARs is discussed in Sect. S1. This approach reduces the sensitivity of AR characteristics to arbitrary strength thresholds, making it easier to compare ARs across different mean state conditions.

To illustrate how SCAFET identifies ARs, we utilized daily mean data from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis Version 5 (ERA5; Hersbach et al., 2020) for the period 2000 to 2019. The key fields of interest included the daily mean integrated water vapor transport (IVT) as the primary field and the daily mean total precipitation as the secondary variable. All the data sets employed share a spatial resolution of 0.25^∘ × 0.25^∘. The vector field, IVT, is calculated as

$$\begin{array}{}\text{(4)}& {\displaystyle }\text{IVT}x=-{\displaystyle \frac{\mathrm{1}}{g}}\underset{\mathrm{1000}\mathrm{hPa}}{\overset{\mathrm{300}\mathrm{hPa}}{\int }}qu\mathrm{d}p\text{(5)}& {\displaystyle }\text{IVT}y=-{\displaystyle \frac{\mathrm{1}}{g}}\underset{\mathrm{1000}\mathrm{hPa}}{\overset{\mathrm{300}\mathrm{hPa}}{\int }}qv\mathrm{d}p\text{(6)}& {\displaystyle }\left|\text{IVT}\right|=\sqrt{\text{IVT}{x}^{\mathrm{2}}+\text{IVT}{y}^{\mathrm{2}}}.\end{array}$$

where u, v, and q are the zonal wind speed, meridional wind speed, and specific humidity, respectively.

To detect AR-like structures, SCAFET employs a search for specific shapes, such as ridges, caps, and domes (see Fig. 2). Following the process outlined in Fig. 1, the SI is calculated after applying a grid-aware smoothing technique that suppresses variability smaller than 1000 km (Fig. 4a). Once the SI is calculated for |IVT| (Fig. 4b), regions where the SI >0.375 are passed on to the next stage for filtering. To maximally utilize the vector qualities of the primary field, we ensure that the local transport direction (arrows in Fig. 4a) and local ridge direction (arrows in Fig. 4b) do not deviate by more than 45^∘. The local ridge direction is identified as the eigenvector corresponding to the smallest eigenvalue (k₂). Filtering based on the grid properties removes candidates that are too small (length <2000 km and area $<\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ km²) or too wide (eccentricity ≤0.75). To eliminate AR-like objects with low strength (precipitation <1 mm d⁻¹), we constrain our results with the secondary field, total precipitation, to within the object's area. The use of precipitation as a strength indicator is relevant, given its significant socioeconomic impact. In line with other ARDTs, we impose a regional mask to filter out AR-like structures along the equatorial belt. All these steps can be applied in parallel along the time axis, and at each time step, AR-like structures similar to those shown in Fig. 4c are identified. Once all ARs are detected, the tracking algorithm is applied to the daily data to filter out ARs that last shorter than 1 d. Tracking is performed based on the centroid of each identified object. The closest objects within a distance of 4000 km between two consecutive time steps are considered to be the same object evolving over time (Fig. 4d). The annual mean frequency of the detected AR objects and their seasonality are shown in Fig. 4e, f, and g. SCAFET's identification of ARs is consistent with other ARDTs, both in terms of detecting single events and determining their mean climatology, as further detailed in the Sect. S2 (see also Lora et al., 2020).

Table 1 The table presents the values for various parameters used in the detection of ARs using SCAFET. The rows for each step, including segmentation, filtering, and tracking are grouped together and labeled.

Figure 4 Major steps in the detection and tracking of atmospheric rivers. (a) Smoothed primary field of vertically integrated water vapor transport (IVT). Smoothing removes variability smaller than 1000 km from the IVT. The arrows in panel (a) represent the direction of unsmoothed IVT. (b) Magnitude (shading) of the shape index (SI) and direction of the local ridge (arrows) direction calculated from smoothed IVT. (c) Labeled AR objects after filtering out weak, small, and ephemeral candidates. (d) Example of tracked AR centroids and marked time; the inlay shows the object's area mean IVT over time. (e) AR annual mean frequency for the period 2000 to 2019. (f, g) AR frequency anomaly relative to the annual mean for (f) November to March and (g) May to September.

3.2 Tropical and extratropical cyclones

In the scientific literature, cyclones are generally described as large weather systems ranging from 500–4000 km in size, characterized by strong cyclonic circulation, low pressure at their center, and exceptionally high winds around them (Emanuel, 2003; Schultz et al., 2019; Encyclopaedia, 2022). The dynamics and characteristics of cyclones can vary, depending on factors such as their genesis location and translation speeds. For instance, cyclones generated near the Equator, commonly referred to as tropical cyclones, are typically smaller in size compared to those formed in midlatitudes, known as extratropical cyclones. Regardless of their origin, cyclones have the potential to unleash intense rainfall and powerful winds along their path and can lead to flooding, landslides, and severe damage to coastal infrastructure when they make landfall (Knutson et al., 2010; Mendelsohn et al., 2012; Ranson et al., 2014). Moreover, the impact of cyclones is becoming a subject of heightened public concern, due to rising sea levels and the potential for increased cyclone intensity in response to global warming. Thus, the identification and future projection of cyclones are a subject of growing attention and importance for the climate community (Woodruff et al., 2013).

Table 2 Same as in Table 1 but for parameters and values relevant to detecting tropical and extratropical cyclones.

Once again, discrepancies among different detection algorithms can be attributed to varying choices of physical thresholds or constraints related to factors such as size, wind speeds, vorticity, or surface pressure anomalies. While most studies generally agree on the present and future characteristics of cyclones, resolving details such as the changes in genesis rate and durations is complicated by the uncertainties in the detection methods (Ulbrich et al., 2009; Neu et al., 2013; Horn et al., 2014; Walsh et al., 2015). In this study, SCAFET identifies cyclones as regions of strong local maxima of cyclonic circulation (SI >0.625) and maximum wind speeds exceeding 10 m s⁻¹. This definition enables the detection of robust cyclonic vorticities worldwide, including but not limited to tropical and extratropical cyclones. The primary field used for cyclone detection is the absolute value of cyclonic relative vorticity (ζ) defined as

$$\begin{array}{}\text{(7)}& \mathit{\zeta }=\mathrm{\nabla }\times \mathit{U},\end{array}$$

where U is the 6 h wind speeds at 10 m above the surface obtained from the ERA5 reanalysis data set, with a spatial resolution of 0.25^∘ × 0.25^∘ (Hersbach et al., 2020). The magnitude of the wind speed at 10 m is utilized as the secondary field to constrain detection. Additional cyclone-related variables such as surface pressure anomaly and potential temperature can also serve as secondary fields for the identification and classification of cyclones.

In contrast with ARs, the detection of cyclones relies on a scalar field, specifically in this case the cyclonic relative vorticity $\left|\mathit{\zeta }\right|$. First, the data are pre-processed with grid-aware Gaussian smoothing to suppress spatial variability smaller than 750 km (Fig. 5a). The chosen smoothing scale allows us to identify both tropical and extratropical cyclones. Caps and dome shapes (SI >0.625) are then identified within the smoothed $\left|\mathit{\zeta }\right|$ field as potential cyclones (Fig. 5b). Subsequently, objects with an area smaller than 10⁵ km² and a diameter shorter than 20 km are filtered out. Once these spatial criteria are met, we can further refine our selection by excluding weak cyclonic vorticities $\left|\mathit{\zeta }\right|<{\mathrm{10}}^{-\mathrm{6}}$ s⁻¹ and slow maximum wind speed <10 m s⁻¹, resulting in the identification of robust cyclonic systems for a given time step (Fig. 5c). Similar to the AR example, all the described steps can be parallelized along the time dimension. Once potential cyclones are identified, they are tracked using a methodology similar to the AR tracking algorithm. However, in this case, the radius for search is limited to 1000 km, since we are using 6 h data, and the translation speeds of cyclones are notably slower than 150 km h⁻¹. A minimum duration of 48 h and a minimum total displacement of 500 km is applied to distinguish moving cyclonic circulations from stationary ones. An example of a tracked cyclone, commonly known as “Dorian” (Avila et al., 2020) is compared with the observed track from IBTrACS (Knapp et al., 2010, 2018) data set (Fig. 5d). In comparison to the observed track, SCAFET's track is much longer, due to the more relaxed conditions applied to cyclonic vorticity and wind speed. Additionally, SCAFET does not differentiate between tropical and extratropical cyclones, which can result in tracking the object throughout its transition from a tropical cyclone to a midlatitude storm. Despite this difference, the long-term averages for cyclone frequency and its seasonal variability calculated using SCAFET are comparable with other studies (e.g., Ullrich and Zarzycki, 2017). What sets SCAFET apart from other conventional cyclone detection algorithms is its approach to identifying cyclones not as point objects but as encompassing surfaces around the point of maximum $\left|\mathit{\zeta }\right|$. This enables a more comprehensive analysis of cyclone properties, including maximum and minimum values of wind speed and precipitation within the entire cyclone structure.

Figure 5 Major steps in the detection and tracking of cyclones. (a) Smoothed primary field of cyclonic relative vorticity ($\left|\mathit{\zeta }\right|$). The smoothing removes variability smaller than 750 km from $\left|\mathit{\zeta }\right|$. (b) Magnitude of the SI for the primary field. (c) Filtered cyclonic objects, with the background color representing unsmoothed values of ζ. (d) Track obtained for cyclone “Dorian” from SCAFET compared with the track from the IBTrACS data set. (e) Annual mean frequency of cyclone occurrence for the period 2000 to 2020. (f, g) Anomalous cyclone frequencies relative to the annual mean for (f) JJA (June–August) and (g) DJF (December–February).

3.3 Sea surface temperature fronts

SST fronts are regions where different water masses come together. They are typically characterized by strong horizontal gradients in temperature, salinity, density, and other properties (Bowman, 1978; Legeckis, 1978; Fedorov, 1986; Yoder et al., 1994). Unlike the larger meso- to synoptic-scale features discussed in this study, frontal structures are often observed in much smaller spatiotemporal scales. Accurate identification of SSTFs is essential because these features are frequently associated with strong upwelling and high levels of biogeochemical productivity (Clayton et al., 2014, 2021; Nagai and Clayton, 2017). Additionally, the detection of SSTFs serves as an example of how SCAFET can be applied to identify features in curvilinear grids.

Many prior SSTF detection algorithms rely on edge detection techniques and the gradient of sea surface temperature and/or height to identify these structures (Canny, 1986; Castelao et al., 2006). In our approach, we utilize the magnitude of the daily mean SST horizontal gradient as the primary field for detecting SST fronts. The SST data are obtained from a fully coupled, ultrahigh-resolution (≊10 km) Community Earth System Model (CESM) v1.2.2 simulation of present-day mean climate (Small et al., 2014; Chu et al., 2020; Nellikkattil et al., 2023). The data are processed by SCAFET in the tripolar (POP) grid. To illustrate the detection process, the analysis focuses on the Kuroshio frontal and extension region for the last 10 years of the simulation.

The extraction of frontal structures using the selected shapes of ridges, caps, and domes is similar to the method for the detection of ARs. Prior to extraction, a spatial smoothing of approximately 30 km is applied. From the extracted SSTF candidates, objects with a mean SST gradient lower than 10⁻⁴ K m⁻¹ are removed. Circular (eccentricity <0.5) and small (area <1000 km²) objects are also filtered out. It is worth noting that, in contrast to AR detection, frontal structures are not tracked. The detected frontal frequency exhibits general patterns and seasonality consistent with findings in previous studies (Xi et al., 2022).

Table 3 Same as in Table 1 but for parameters and values relevant to detecting sea surface temperature fronts (SSTFs).

Figure 6 Major steps in the detection of sea surface temperature fronts (SSTFs). (a) Magnitude of the SI as calculated from the smoothed primary field of the horizontal gradient of sea surface temperature (∇SST). Smoothing removes variability smaller than 15 km from ∇SST. (b) Filtered SSTF objects, in units of Kelvins per kilometer (K km⁻¹), where the background color represents unsmoothed values of ∇SST. (c) Annual mean frequency of SSTF occurrence across a 10-year period in the present climate simulation. (d, e) Anomalous frontal frequencies relative to the annual mean for (d) JJA and (e) DJF.

4 Application to 3D features

This section introduces the extension of SCAFET to detect features within three-dimensional (3D) primary fields. The process of scale–space selection involves applying Gaussian smoothing independently along each of the three dimensions. Notably, a 3D field yields three eigenvalues ($k_{\mathrm{1}}\ge k_{\mathrm{2}}\ge k_{\mathrm{3}}$) instead of the usual two. In this context, the SI can be calculated in three different ways by combining these eigenvalues.

For the extraction of jet streams, the SI calculated using k₁ and k₂ (the two largest eigenvalues) is used, as it provides a more conservative estimate for the jet-like structure (see Appendix A3 and Fig. S7). The decision to exclude the smallest eigenvalue, denoted as k₃, is based on empirical observations. Empirical evidence suggests that when dealing with regions exhibiting positive maxima (convex curvature), both SI(k₁,k₂) and SI(k₁,k₃) effectively capture the shape. Meanwhile, SI(k₂,k₃) has a trivial application (refer to Fig. A4). Conversely, for concave shapes, both SI(k₁,k₃) and SI(k₂,k₃) represent the shape, while the conditions for SI(k₁,k₂) become redundant, given that they are satisfied by SI(k₁,k₃) and SI(k₂,k₃).

4.1 Jet streams

Jet streams, regardless of the underlying dynamics, are narrow, high-wind-speed regions in the upper atmosphere, with faster wind speeds compared to their surroundings (Koch et al., 2006). These jet streams have a significant impact on aviation and strongly influence surface weather conditions. For example, a persistent jet stream in boreal summer can result in extreme heat and flooding events, while a meandering jet stream in winter leads to severe cold spells in the midlatitudes (Petoukhov et al., 2013; Coumou et al., 2014; Kretschmer et al., 2016). Additionally, the northward movement of jet streams due to greenhouse warming contributes to the poleward propagation of tropical cyclones (Studholme et al., 2021). Thus, accurately detecting and characterizing jet streams is crucial for predicting and projecting both climatology and extreme weather systems.

Table 4 Same as in Table 1 but for parameters and values relevant for detecting jet streams.

Much like the detection of other weather phenomena discussed in this study, previous research typically employs a physical threshold to identify jet streams. Furthermore, with the exceptions of Limbach et al. (2012) and Kern et al. (2018), most studies identify jet streams as either one- or two-dimensional features. However, it is important to emphasize that this section's focus is primarily on illustrating the method for detecting jet streams rather than validation of any analysis with published work. There is currently limited analysis available for comparing with a 3D perspective of jet streams, highlighting the need for such an approach. As a result, we present examples of jet stream detection in three selected time steps. A more comprehensive analysis and discussion regarding of the long-term characteristics of jet streams will be a topic for future research. For those interested, a video showcasing the results over an extended period can be found in the Supplement.

The primary field used in the extraction of jet streams is the 6 h, three-dimensional wind speeds obtained from ERA5 reanalysis data set, with a spatial resolution of 1^∘ with 37 vertical levels (Hersbach et al., 2020). The magnitude of wind speed is calculated as

$$\begin{array}{}\text{(8)}& W=\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}},\end{array}$$

where u and v are the zonal and meridional wind velocities.

Figure 7 The 3D jet streams extracted using SCAFET. Magnitude of 3D wind speed for (a) 28 August 2022 at 12:00 UTC, (c) 28 August 2022 at 18:00, and (e) 29 August 2022 at 00:00 UTC. Extracted 3D jet streams for corresponding periods are shown in panels (b), (d), and (f), respectively. The reader is encouraged to view the full video of these snapshots in the Supplement.

The detection process for jet streams begins similarly to the detection of 2D features. Gaussian smoothing is used to remove variability shorter than 3000 km in the horizontal dimensions. No smoothing is applied along the vertical dimension. Next, the SI is calculated using the two largest eigenvalues, k₁ and k₂. The vertical dimension for the three-dimensional wind speed is given in pressure coordinates. To calculate the gradient as change in wind speeds per kilometer, a rudimentary conversion from pressure to height coordinates is used (refer to Wallace and Hobbs, 1977, p. 60–61, and https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.pressure_to_height_std.html, last access: 20 December 2023, for further details).

Similar to the detection of ARs, regions characterized by the selected shapes of ridges, caps, and domes (SI >0.375) are isolated for filtering. Filtering is then applied to remove objects with a volume less than 10⁶ km³, a horizontal length less than 5000 km, and a maximum wind speed within each object less than 50 m s⁻¹. In the current version of SCAFET, the tracking algorithm is not applied to jet detection (see Fig. 7). The detailed list of parameters used in the detection of jet streams is given in Table 4.

5 Conclusions

In this study, we introduced a novel computational mathematical framework and an open-source Python package for extracting and tracking features from large climate data sets, called Scalable Feature Extraction and Tracking (SCAFET). The purpose of SCAFET is to tackle the challenges posed by the increasing volume and diversity of climate data by providing an alternative to traditional physical threshold-based feature detection methods. It enables the comparison of features between observational and model data with different mean states by attempting to remove the need for posterior data-specific assumptions. Furthermore, SCAFET introduces a novel shape-based approach to feature extraction, which helps uncover discrepancies in climate projections due to differences in detection methods and aims to help the community in building scientific consensus. To demonstrate SCAFET's capabilities and its potential in advancing these goals, we showcased its ability to detect various features, including two-dimensional features such as atmospheric rivers (ARs), tropical and extratropical cyclones, sea surface temperature fronts, and the detection of three-dimensional jet streams. Each application serves as an illustrative example from which users can customize SCAFET for their specific research needs.

SCAFET offers several significant advantages, including a more comprehensive framework and parallel computing implementation for efficiency. However, its most noteworthy contribution lies in offering a novel perspective on how we can relatively define various features within climate data sets that span extensive periods marked by significant changes in mean climate. Rather than relying on empirically derived, data-specific physical thresholds for feature extraction, SCAFET identifies features using shape-based absolute thresholds and the locally estimated shape within the field. This methodology offers a unique viewpoint, enabling us to observe the continuous changes in feature properties while accounting for shifts in the mean climate state. This approach is particularly valuable, as meso–synoptic-scale studies are highly sensitive to thresholds in a dynamically changing mean climate state. Consequently, the conclusions drawn from such studies can vary significantly, as demonstrated in research examining the response of ARs to greenhouse warming (Zhao, 2020; O'Brien et al., 2022; Nellikkattil et al., 2023). Thus, algorithms like SCAFET which are not influenced by data-specific conditions of various climate models play a crucial role in advancing scientific understanding and facilitating climate model development.

In conclusion, delving deeper into the principles of differential geometry to elucidate the physical interpretation of the relationship between the SI and local geometric shape has the potential to revolutionize our approach to feature extraction from large data sets. This avenue of research has the promise of significantly enhancing the algorithm's robustness and reliability. It is worth noting that, at present, SCAFET may not surpass the computational efficiency of other well-established feature extraction methods discussed above (see Sect. S2.2). However, ongoing efforts to optimize and streamline the algorithm for improved computational efficiency continue. One notable strength of SCAFET's design is its independence from data-set-specific posterior information when identifying features. Moreover, the shape-based thresholds used for detecting specific features remain consistent across various grids, data sets, and climatologies. Between these strengths and the full parallelization of the feature detection method, there are exciting possibilities for further development. This may eventually enable the algorithm to be used in operational feature identification and early-warning systems for extreme weather events.

Appendix A: Shape-based feature extraction on simple data sets

This section demonstrates how shape-based feature extraction can be performed on scalar fields represented by simple, idealized mathematical functions. It is intended to provide readers with more insights into the basic principles behind shape-based feature extraction and how it differs from other conventional methods. We have also tried to showcase some properties of shape-based feature extraction methods like its insensitivity to linear mean state trends.

A1 Application to 1D data sets

In this section, we draw an analogy between the use of SCAFET on a two-dimensional (2D) data set and shape-based feature extraction from a one-dimensional (1D) data set. Our intention is not to promote the use of shape-based extraction of features from 1D data sets but rather to provide readers with a fundamental understanding of this approach, along with its strengths and limitations.

For any differentiable curve C, the curvature is measured as the instantaneous rate of change in direction along the curve. Simply put, the curvature is measured as the rate of change in the unit tangent to the curve at any given point. An osculating circle can be used to intuitively represent the curvature of a surface or a curve (see Fig. A1). At any point P, the curvature, k is the reciprocal of the radius (R) of the circle. The sign of k determines if the curve has a concave or a convex curvature. More information and mathematical proof for these concepts can be found in any standard differential geometry textbook.

Following the derivation of the shape index (SI) for 2D data sets, we calculate the local shape of a function f using the shape parameter, defined as

$$\begin{array}{}\text{(A1)}& K={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}{\tan }^{-\mathrm{1}}\left(f^{\prime \prime }\right).\end{array}$$

Values of K closer to 1 are identified as regions of local minima, while K values closer to −1 are regions of local maxima (black curve in Fig. A2). Depending on the magnitude of the function, one could adjust the value of K to obtain regions of local maxima (red caps in Fig. A2) and local minima (green caps in Fig. A2). The curvature of the function is insensitive to linear trends and mean state changes. This is evident, as the application of identical shape thresholds identifies the same regions of the curves as local maxima and minima, whether on the base curve (blue curve in Fig. A2) or on the same curve with an added linear trend (orange curve in Fig. A2). The values of K for both curves are represented by the black line in Fig. A2. Thus, the shape parameter can be used to identify the local minima and maxima from a 1D data set, despite background state changes.

Figure A1 Schematic representation of curvature measurement of a curve C at point P. At P, the curvature is the reciprocal of the radius R of the osculating circle. In differential geometry, an osculating circle is defined as the circle passing through the point P and a pair of additional points infinitesimally close to P.

Figure A2 Comparison of shape extraction between a simple one-dimensional curve, given by $f=\sin \mathrm{2}x+\mathrm{3}\cos \mathrm{5}x$ (blue curve) and $f+\mathrm{0.5}x=\sin \mathrm{2}x+\mathrm{3}\cos \mathrm{5}x+\mathrm{0.5}x$ (orange; the blue curve with a linear trend). The left y axis shows the magnitude of both functions, and the right y axis indicates the values of the shape parameter (K). Note that value of K is the same for both functions. The green and red highlighting on the curves shows regions where K>0.99 and $K<-\mathrm{0.99}$, corresponding to regions of local maxima and minima, respectively.

A2 Application of SCAFET to simple geostrophic motion

In this section, we apply SCAFET to a basic geostrophic rotational motion. The goal of this discussion is to illustrate how the shape-based extraction of 2D features differs from conventional methods. The calculation of the SI involves the computation of the two eigenvalues, k₁ and k₂, of the Hessian matrix for any gridded data set. As discussed in the previous section, the curvature measurement provided by k₁ and k₂ can be visualized as the reciprocal of the radius of two osculating circles that intersect orthogonally at a point on the surface. Large negative eigenvalues signify surfaces with strong convex curvature, while positive values correspond to troughs or cups.

To demonstrate the characteristics and advantages of feature detection based on the SI, let us consider a simple rotational wind field (see Fig. A3a vectors) given by

$$\begin{array}{}\text{(A2)}& {\displaystyle }{u}_{\mathrm{g}}=-\mathrm{\Omega }y\text{(A3)}& {\displaystyle }{v}_{\mathrm{g}}=\mathrm{\Omega }x,\end{array}$$

where Ω is a constant ($\mathrm{\Omega }={\mathrm{10}}^{\mathrm{5}}{\mathrm{rad}\mathrm{s}}^{-\mathrm{1}}$), and x, y represents the grid. The geopotential height (h) of the field (see Fig. A3a shading) is used as our primary field in calculating the SI, which is computed as

$$\begin{array}{}\text{(A4)}& h={\displaystyle \frac{\mathrm{\Omega }f}{\mathrm{2}g}}(x^{\mathrm{2}}+y^{\mathrm{2}}),\end{array}$$

where f and g are the Coriolis parameter and the acceleration due to gravity, respectively. The SI is calculated from the eigenvalues of the Hessian of h, using the formula

$$\begin{array}{}\text{(A5)}& \text{SI}(k_{\mathrm{1}},k_{\mathrm{2}})={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}{\tan }^{-\mathrm{1}}\left[{\displaystyle \frac{k_{\mathrm{2}}+k_{\mathrm{1}}}{k_{\mathrm{2}}-k_{\mathrm{1}}}}\right],\end{array}$$

where the eigenvalues k₁ and k₂ are given by

$$\begin{array}{}\text{(A6)}& k_{\mathrm{12}}={\displaystyle \frac{f{\mathit{\zeta }}_{\mathrm{g}}}{\mathrm{2}g}}\pm \sqrt{{\left({\displaystyle \frac{f}{\mathrm{2}g}}\right)}^{\mathrm{2}}-{\left({\displaystyle \frac{f}{g}}\right)}^{\mathrm{2}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial y}}+{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial y}}},\end{array}$$

where ζ_g is the geostrophic vorticity. This then gives SI as

$$\begin{array}{}\text{(A7)}& \mathrm{\text{SI}}(k_{\mathrm{1}},k_{\mathrm{2}})={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}{\tan }^{-\mathrm{1}}\left[{\displaystyle \frac{{\mathit{\zeta }}_{\mathrm{g}}}{-\mathrm{2}\sqrt{{\left(\frac{{\mathit{\zeta }}_{\mathrm{g}}}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\partial v_{\mathrm{g}}}{\partial x}\frac{\partial u_{\mathrm{g}}}{\partial y}+\frac{\partial u_{\mathrm{g}}}{\partial x}\frac{\partial v_{\mathrm{g}}}{\partial y}}}}\right].\end{array}$$

A detailed derivation of the above equation can be found in Sect. B. By plugging in the values for the rotational motion, we get

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\mathit{\zeta }}_{\mathrm{g}}={\mathrm{\nabla }}^{\mathrm{2}}h=\mathrm{\Omega }f/g\text{(A9)}& {\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial x}}={\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial y}}=\mathrm{0}\text{(A10)}& {\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial y}}={\mathrm{\Omega }}^{\mathrm{2}}.\end{array}$$

Therefore,

$$\begin{array}{}\text{(A11)}& \text{SI}={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}{\tan }^{-\mathrm{1}}\left[{\displaystyle \frac{{\mathrm{\Omega }}^{\mathrm{2}}}{-\sqrt{{\mathrm{\Omega }}^{\mathrm{2}}-{\mathrm{\Omega }}^{\mathrm{2}}}}}\right]=-\mathrm{1}.\end{array}$$

Thus, SCAFET classifies the whole domain with counterclockwise rotational motion as a trough with the SI $\approxeq -\mathrm{1}$, regardless of the absolute value of the field or Ω. In contrast, traditional methods that rely on thresholding the geopotential height would identify regions based on the chosen threshold of h, which would need to be adjusted depending on the mean (time) and background (space) state. Another common approach is to establish a threshold on the smallest eigenvalue, thus aiming to identify extreme features based on the curvature strength rather than the field's actual value. TempestExtremes (Ullrich and Zarzycki, 2017), a feature extraction framework discussed in the main text, follows this method to detect atmospheric rivers from gridded data sets. In the current example, this approach would correspond to setting a threshold on $f\mathrm{\Omega }/g$. In other words, TempestExtremes would only identify a trough if the value of Ω exceeds the pre-determined threshold. SCAFET, on the other hand, identifies the trough region as a trough regardless of the specific value of the field or Ω. This illustrates how feature extraction, using SI and other published methods, can yield different results, depending on the input data, as they focus on distinct properties of the field.

Figure A3 Comparison between two feature extraction techniques on an idealized example of rotational wind field. (a) The geopotential height (h) (shading) of the rotational wind field (arrows). h is defined as $\mathrm{\Omega }f(x^{\mathrm{2}}+y^{\mathrm{2}})/\mathrm{2}g$, where $f={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, g=9.805 m s⁻¹, and Ω=10⁵ rad s⁻¹. (b) Magnitude of the smallest eigenvalue, derived from the equation as $f\mathrm{\Omega }/g=\mathrm{1.0199}$, thus illustrating a uniform field as expected. (c) Value of the SI where $\text{SI}=-\mathrm{1}$ throughout the domain, as expected.

Figure A4 Various approaches for extraction of a 3D cylinder from a scalar field. (a) A simple scalar field represented by $\sin \left(\mathrm{3}X\right)+\cos \left(\mathrm{4}Y\right)*\cos \left(Z\right)$ is shown. (b–d) The extracted cylinders by applying the conditions in panel (b) $\text{SI}(k_{\mathrm{1}},k_{\mathrm{2}})>\mathrm{0.375}$, while in panel (c) $\text{SI}(k_{\mathrm{1}},k_{\mathrm{3}})>\mathrm{0.375}$ and panel (d) $\text{SI}(k_{\mathrm{2}},k_{\mathrm{3}})>\mathrm{0.375}$ are shown. The values enclosed in parentheses within the figure titles indicate the percentage of data that satisfies the respective conditions applied in each case.

A3 Application of SCAFET to 3D fields

This section aims to demonstrate the detection of a cylindrical volume within a three-dimensional scalar field. To illustrate the effectiveness of the SI in identifying 3D structures embedded within scalar fields, we offer a straightforward example of how the SI can be used to isolate a cylinder embedded in a scalar field defined by $f=\sin \left(\mathrm{3}X\right)+\cos \left(\mathrm{4}Y\right)\cos \left(Z\right)$. It is worth mentioning that this specific problem bears significant similarities to the task of identifying 3D jet cores.

As explained in Sect. 4, a three-dimensional field provides us with three eigenvalues, thus satisfying the condition $k_{\mathrm{1}}\ge k_{\mathrm{2}}\ge k_{\mathrm{3}}$. The SI can be computed using SI(k₁,k₂), SI(k₁,k₃), or SI(k₂,k₃). Setting a threshold of the SI >0.375 effectively isolates the cylinder when using either SI(k₁,k₂) or SI(k₁,k₃) (see Fig. A4b–d). Between these two options, SI(k₁,k₂), which utilizes the two largest eigenvalues, imposes a more conservative criterion for identifying the embedded cylinder. The percentage of data identified as the cylinder is provided in the title of each plot in Fig. A4. Notably, employing SI(k₂,k₃) is not suitable, as it fails to isolate the desired cylinder shape effectively. The choice of using SI(k₁,k₂) is specifically tailored for extracting convex shapes or local maxima. Interestingly, to identify concave shapes or local minima, one should utilize the SI derived from the two smallest eigenvalues, namely SI(k₂,k₃).

While the simple example presented here may not provide a comprehensive illustration of 3D feature detection, we hope that it encourages further fundamental research into 3D feature extraction to expand the capabilities of analysis and increase precision.

Appendix B: Derivation of shape index for geostrophic motion

The complete derivation of the SI for geostrophic wind fields is shown in this section. The result from the derivation is used in Sect. A.

Let h be the geopotential height at a certain level. The Hessian of h is given by

$$\begin{array}{}\text{(B1)}& \mathcal{H}\left(h\right)=\left(\begin{array}{cc}\frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}& \frac{{\partial }^{\mathrm{2}}h}{\partial x\partial y}\\ \frac{{\partial }^{\mathrm{2}}h}{\partial y\partial x}& \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}\end{array}\right).\end{array}$$

The eigenvalues of the symmetric matrix ℋ are calculated by solving the quadratic equation.

$$\begin{array}{}\text{(B2)}& \left({\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}}-\mathit{\lambda }\right)\left({\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}}-\mathit{\lambda }\right)-{\left({\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x\partial y}}\right)}^{\mathrm{2}}=\mathrm{0},\end{array}$$

which can be expanded as

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\mathit{\lambda }}^{\mathrm{2}}-\mathit{\lambda }\left({\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}}\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}}.{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}}-{\left({\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x\partial y}}\right)}^{\mathrm{2}}& {\displaystyle }=\mathrm{0}\text{(B4)}& {\displaystyle }{\mathit{\lambda }}^{\mathrm{2}}-\mathit{\lambda }{\mathrm{\nabla }}^{\mathrm{2}}h+{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}}.{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}}-{\left({\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x\partial y}}\right)}^{\mathrm{2}}& {\displaystyle }=\mathrm{0}.\end{array}$$

Note that the geostrophic vorticity (ζ_g) is defined as

$$\begin{array}{}\text{(B5)}& {\mathit{\zeta }}_{\mathrm{g}}={\displaystyle \frac{g}{f}}{\mathrm{\nabla }}^{\mathrm{2}}h.\end{array}$$

The geostrophic velocities are defined as

$$\begin{array}{}\text{(B6)}& {\displaystyle }{u}_{\mathrm{g}}=\mathrm{-}{\displaystyle \frac{g}{f}}{\displaystyle \frac{\partial h}{\partial y}}=-{\displaystyle \frac{\partial \mathit{\psi }}{\partial y}}\text{(B7)}& {\displaystyle }{v}_{\mathrm{g}}={\displaystyle \frac{g}{f}}{\displaystyle \frac{\partial h}{\partial x}}={\displaystyle \frac{\partial \mathit{\psi }}{\partial x}},\end{array}$$

where ψ is the geostrophic stream function. This implies

$$\begin{array}{}\text{(B8)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial x^{\mathrm{2}}}}={\displaystyle \frac{f}{g}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial x^{\mathrm{2}}}}={\displaystyle \frac{f}{g}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial x}}\text{(B9)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial y^{\mathrm{2}}}}={\displaystyle \frac{f}{g}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial y^{\mathrm{2}}}}=-{\displaystyle \frac{f}{g}}{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial y}}.\end{array}$$

Adding the abovementioned relationships to Eq. (3) leads to the following:

$$\begin{array}{}\text{(B10)}& {\mathit{\lambda }}^{\mathrm{2}}-{\displaystyle \frac{\mathit{\lambda }f}{g}}{\mathit{\zeta }}_{\mathrm{g}}-{\displaystyle \frac{f^{\mathrm{2}}}{g^{\mathrm{2}}}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{g}}}{\partial y}}+{\displaystyle \frac{f^{\mathrm{2}}}{g^{\mathrm{2}}}}{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{g}}}{\partial y}}.\end{array}$$

By solving for λ, we get

$$\begin{array}{}\text{(B11)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{12}}={\displaystyle \frac{f{\mathit{\zeta }}_{\mathrm{g}}}{\mathrm{2}g}}\pm \sqrt{{\left({\displaystyle \frac{f}{\mathrm{2}g}}\right)}^{\mathrm{2}}-{\left({\displaystyle \frac{f}{g}}\right)}^{\mathrm{2}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial y}}+{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial y}}}\text{(B12)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{12}}={\displaystyle \frac{f}{g}}\left[{\displaystyle \frac{{\mathit{\zeta }}_{\mathrm{g}}}{\mathrm{2}}}\pm \sqrt{{\left({\displaystyle \frac{{\mathit{\zeta }}_{\mathrm{g}}}{\mathrm{2}}}\right)}^{\mathrm{2}}-{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial y}}+{\displaystyle \frac{\partial u_{\mathrm{g}}}{\partial x}}{\displaystyle \frac{\partial v_{\mathrm{g}}}{\partial y}}}\right].\end{array}$$

Thus, the shape index for h is

$$\begin{array}{}\text{(B13)}& \text{SI}={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\tan \left[{\displaystyle \frac{{\mathit{\zeta }}_{\mathrm{g}}}{-\mathrm{2}\sqrt{{\left(\frac{{\mathit{\zeta }}_{\mathrm{g}}}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\partial v_{\mathrm{g}}}{\partial x}\frac{\partial u_{\mathrm{g}}}{\partial y}+\frac{\partial u_{\mathrm{g}}}{\partial x}\frac{\partial v_{\mathrm{g}}}{\partial y}}}}\right].\end{array}$$