TempestExtremes : A Framework for Scale-Insensitive Pointwise Feature Tracking on Unstructured Grids

This paper describes a new open-source software framework for automated pointwise feature tracking that is applicable to a wide array of climate datasets using either structured or unstructured grids. Common climatological pointwise features include tropical cyclones, extratropical cyclones and tropical easterly waves. To enable support for a wide array of detection 5 schemes, a suite of algorithmic kernels have been developed that capture the core functionality of algorithmic tracking routines from throughout literature. A review of efforts related to pointwise feature tracking from the past three decades is included. Selected results using both reanalysis datasets and unstructured grid simulations are provided.


Introduction
Automated pointwise feature tracking is an algorithmic technique for objective identification and tracking of meteorological features, such as extratropical cyclones, tropical cyclones and tropical easterly waves, and has emerged as an important and desirable data processing capability in climate science.Software tools for feature tracking -typically referred to as "trackers" -have been employed to evaluate model performance and answer pressing scientific questions regarding anticipated changes in atmospheric features under climate change.Exploration of tracker literature has exposed a breadth of potential techniques that have been applied to climate datasets with varied spatial resolution and temporal frequency (a comprehensive review of the tracking literature can be found in Appendices A, B and C).Nonetheless, the definition of an optimal objective criteria for key atmospheric features has eluded development, and ambiguity in the formal definition of these features suggests that there may be no singular criteria capable of both perfect detection and zero false alarm rate.Further, as observed by Walsh et al. (2007) and Horn et al. (2014) for tropical cyclones and Neu et al. (2013) for extratropical cyclones, feature tracking schemes can produce wildly varying results depending on the specific choice of threshold variables and values.Therefore, we argue that uncertainties associated with objective tracking criteria should be addressed with an ensemble of detection thresholds and variables, whereas blind application of singular tracking formulations should be avoided.To this end, it is the goal of this paper to review the vast literature on trackers and use this information to inform the development of a unified framework (TempestExtremes) that enables a variety of tracking procedures to be quickly and easily applied across arbitrary spatial resolution and temporal frequency.This paper focuses largely on the technical aspects of pointwise feature tracking, but sets the stage for future studies on parameter sensitivity and optimization.
Most algorithmic Lagrangian trackers of pointwise features (such as cyclones and eddies) share a common procedure: 1.They must identify an initial set of candidate points by searching for local extrema.Local extrema can be further specified, for instance, by requiring that they be sufficiently anomalous when contrasted with their neighbors.For most cyclonic structures, either minima in the sea level pressure field or maxima in the absolute value of the relative vorticity are used.
P. A. Ullrich and C. M. Zarzycki: A framework for scale-insensitive pointwise feature tracking 2. They must eliminate candidate points that do not satisfy a prescribed set of thresholds.For instance, tropical cyclones typically require the presence of an upper-level warm core that is sufficiently near the sea level pressure minima that defines the storm center.Additional criteria, such as a minimum threshold on relative vorticity, can be used to eliminate spurious detections.
3. They must connect candidate points together in time (referred to as "stitching") to generate feature paths, eliminating paths that are of insufficient length or do not meet additional criteria.
Although the actual implementations of this procedure does vary throughout the literature, a review of this material reveals several core algorithms (kernels) that are common across implementations.Based on our analysis, the five most commonly employed kernels are computation of anomalies in a data field from a spatially averaged mean; identification of local extrema in a given 2-D data field (for instance, sea level pressure minima); determination of whether a closed contour exists in a data field around a particular point; determination of whether, in the neighborhood of a particular point, a data field satisfies a given threshold; and stitching of candidates from sequential time slices to build feature tracks.
The development of a robust implementation of these five kernels will be the focus of the remainder of this paper.Feature tracking that is robust across essentially arbitrary datasets requires some additional considerations.Detection criteria and thresholds for tracking are often tuned based on the characteristics of a particular dataset, such as temporal resolution, spatial resolution and regional coverage.Unfortunately, this has led to an abundance of schemes that often cannot be directly compared, or applied in a more general context.To this end, we focus on kernels that are insensitive to the characteristics of the input data.For instance, averaging or searching over a discrete number of grid points around each candidate (a common approach) is incompatible with scale insensitivity since the physical search radius would be dependent on the spatial resolution of the data.Identification of local extrema is also a resolution-sensitive procedure, since the number of extrema will often scale with the number of spatial data points; however, a closed contour criteria based on a physical distance is largely resolution insensitive.To achieve robust applicability, a general framework should use great-circle arcs for all distance calculations (this avoids issues associated with latitude-longitude distance that emerges near the poles); support structured and unstructured grids (this eliminates the need for post-processing of large native grid output files and enables detection and characterization simultaneous with the model execution); not contain hard-coded variable names, so as to ensure robust applicability across reanalysis datasets and applicability to a variety of problems; and allow for easy intercomparison of detection schemes by enabling detection criteria and thresholds that are compactly specified on the command line.
Well-known automated software trackers include TRACK (Hodges, 1994(Hodges, , 1995(Hodges, , 2015) ) and the Geophysical Fluid Dynamics Laboratory (GFDL) TSTORMS package (Vitart et al., 1997;Zhao et al., 2009).Both of these software packages have been used extensively to examine pointwise features in the atmosphere, but do not completely satisfy the four requirements above.The remainder of the paper is organized as follows: Sect. 2 describes the algorithms and kernels that have been implemented in the TempestExtremes software framework.Selected examples of tropical cyclone, extratropical cyclone and tropical easterly wave detection are then provided in Sect.3, followed by conclusions in Sect. 4. The appendices provide a review of relevant literature on pointwise feature trackers of extratropical cyclones (Appendix A), tropical cyclones (Appendix B) and tropical easterly waves (Appendix C).A technical guide to the use of the Tempes-tExtremes tools DetectCyclonesUnstructured and StitchNodes is provided in Appendices D and E. Additional examples and updates are available as part of the software package.

TempestExtremes algorithms and kernels
This section describes the key building blocks that have been developed in constructing our detection and characterization framework.Pseudocode is utilized throughout to describe the structure of each algorithm.

Unstructured grid specification
For purposes of determining connectivity of the unstructured grid, we require the specification of a node graph that enumerates the neighbors of each node (one such node graph is depicted in Fig. 1).The connectivity information is stored textually as an adjacency list via a variable-length commaseparated variable file.The total number of nodes (N ) is specified at the top of the file, followed by N lines containing the longitude (lon), latitude (lat), associated area, number of adjacent nodes and finally a 1-indexed list of all adjacent nodes, such as depicted below: <total number of nodes>, <lon>,<lat>, <area>, <# adj.nodes>, <first adj.node>,..,<last adj.node> ...

Great-circle distance
As mentioned earlier, in order to avoid sensitivity of the detection scheme to grid resolution, great-circle distance has been employed throughout.In terms of regular latitudelongitude coordinates, the great-circle distance (r), for a sphere of radius a, between points (λ 1 , ϕ 1 ) and (λ 2 , ϕ 2 ), is defined via the symmetric operation: (1) Algorithmically, this calculation is implemented as gcdist(p,q) for given graph nodes p and q.The difference between great-circle distance and latitudelongitude distance is striking at high latitudes, as depicted in Fig. 2. Whereas latitude-longitude distance is generally sufficient for tropical cyclone detection, tracking of high-latitude phenomena such as extratropical cyclones is expected to be more consistent when great-circle distance is employed.

Efficient neighbor search using k-d trees
Three-dimensional (k = 3) k-d trees (Bentley, 1975) are used throughout our detection code using the implementation of Tsiombikas (2015).instead of great-circle distance, we utilize the observation that straight-line and great-circle distance maintain the same ordering for points confined to the surface of the sphere.
Three key functions made available by the k-d tree structure are used: K = build_kd_tree(P) constructs a k-d tree K from a point set P. q = kd_tree_nearest_neighbor(K, p) locates the nearest neighbor q to point p using the k-d tree K. S = kd_tree_all_neighbors(K, p, dist) locates all points that are within a distance dist of a point p within the k-d tree K.
An example k-d tree in two dimensions is depicted in Fig. 3, along with a brief description of how the nearest-neighbor algorithm is performed.

Computing a spatially averaged mean
Many existing tracking algorithms use either a spatially averaged mean field (over a given distance) or an anomaly field computed against the spatially averaged mean (Haarsma et al., 1993;Bengtsson et al., 1995).Algorithm 1 Compute the spatial mean value of a field G over a region of radius dist using graph search on an unstructured grid.
field F = mean (field G, dist) for each node p total_area = 0 while tovisit is not empty q = remove node from tovisit add q to visited F for each neighbor s of q if (gcd(p,s) < dist) and (s is not in visited) then add s to tovisit Algorithm 2 Locate the set of all nodes P that are local minima for a field G (for instance, SLP) defined on an unstructured grid.The procedure for locating maxima is analogous.
Algorithm 3 Given a field G defined on an unstructured grid and a set of candidate points P, remove candidate minima that are within a distance dist of a more extreme minimum, and return the new candidate set Q.

Extrema detection
For purposes of computational efficiency, candidate points are initially located by identifying local extrema in a given field (for instance, sea level pressure (SLP)) via find_all_minima (Algorithm 2).This algorithm compares all nodes against their associated neighbors and only tags points that are less/greater than those neighbors.Candidates are then eliminated if they are "too close" to stronger extrema (Algorithm 3) (e.g., Pinto et al., 2005).This algorithm is performed by eliminating candidate nodes that are within a given distance of a stronger extrema.The initial search field is specified to TempestExtremes either via the --searchbymin or --searchbymax command line argument.The merge distance used in merge_candidates_minima is specified via the --mergedist command line argument.for each node p total_area = 0 while tovisit is not empty q = remove node from tovisit add q to visited F[p] = F[p] + G[q] * area[q] total_area = total_area + area[q] for each neighbor s of q if (gcd(p,s) < dist) and (s is not in visited) then add s to tovisit F[p] = F[p] / total_area Algorithm 2 Locate the set of all nodes P that are local minima for a field G (for instance, SLP) defined on an unstructured grid.The procedure for locating maxima is analogous.
Algorithm 3 Given a field G defined on an unstructured grid and a set of candidate points P, remove candidate minima that are within a distance dist of a more extreme minimum, and return the new candidate set Q.
set Q = merge_candidates_minima(field G, set P, dist) Taylor, M. A., and Ullrich, P. A.: Aquaplanet experiments using CAM's variable-resolution dynamical core, J. Climate, 27, 5481-5503, doi:10.1175/JCLI-D-14-00004.1, 2014b. Zhao, M., Held, I. M., Lin, S.-J., andVecchi, G. A.: Simulations 40 of global hurricane climatology, interannual variability, and response to global warming using a 50-km resolution GCM, J. Climate, 22, 6653-6678, 2009. Zolina, O. and  That is, the distance between a node and its neighbors decreases proportionally to the local grid spacing, and so does not define a "physical" criterion.Consequently, we instead advocate for a "closed contour criteria" to define candidate nodes.Closed contours were first employed by Bell and Bosart (1989), who used a 30 m 500 hPa geopotential height contour to identify closed circulation centers.Their approach used radial arms generated at 15 • intervals over a great-circle distance of 2 • and required that geopotential heights rise by at least 30 m along each arm.Unfortunately, the use of radial arms to define the closed contour is again sensitive to model resolution, since it has the potential to only sample as many neighbors as radial arms employed.
Here, we propose an alternative closed contour criteria that is largely insensitive to model resolution, using graph search to ensure that all paths along the unstructured grid from an initial location p0 lead to a sufficiently large decrease (or increase) in a given field G before reaching a specified radius.This criteria is illustrated in Fig. 4, and is implemented in Algorithms 4 and 5 (for closed contours around local maxima).The closed contour criteria is implemented in TempestExtremes via the command line argument --closedcontourcmd.An analogous command line argument, --noclosedcontourcmd, is also provided, which has similar functionality but discards candidates that satisfy the closed contour criteria (this may be desirable, for instance, to identify cyclonic structures that do not have a warm core).
Geosci node pmax = find_max_near(node p, field G, maxdist) set visited = {} set tovisit = {p} pmax = p while tovisit is not empty q = remove node from tovisit if (q in visited) then continue add q to visited if (gcdist(p,q) > maxdist) then continue if (G[q] > G[pmax]) then pmax = q Algorithm 5 Determine if there is a closed contour in field G of magnitude thresh around the point p0, defined by p0 = find_max_near(p, G, maxdist), within distance dist.That is, along all paths away from p0, the field G must drop by at least thresh within distance dist.
The closed contour criteria is depicted in Fig. 54.An analogous procedure is defined for closed contours around minima.
closed_contour_max(point p, field G, dist, maxdist, thresh) p0 = find_max_near(p, G, maxdist) set visited = {} set tovisit = {p0} while tovisit is not empty q = remove point from tovisit if (q in visited) then continue add q to visited if (gcdist(p0,q) > dist) then return false if (G[p0] -G[q] < thresh) then add all neighbors of q to tovisit return true Algorithm 6 Determine if a candidate node p satisfies the requirement that there exists another node p0 within distance dist of p with G[p] > thresh.
) and (gcdist(p,q) < dist) then add q into s remove q from P[u] p = q else if (gap < maxgap) then gap = gap + 1 else break add s into S node pmax = find_max_near(node p, field G, maxdist) set visited = {} set tovisit = {p} pmax = p while tovisit is not empty q = remove node from tovisit if (q in visited) then continue add q to visited if (gcdist(p,q) > maxdist Algorithm 5 Determine if there is a closed contour in field G of magnitude thresh around the point p0, defined by p0 = find_max_near(p, G, maxdist), within distance dist.That is, along all paths away from p0, the field G must drop by at least thresh within distance dist.
The closed contour criteria is depicted in Fig. 4.An analogous procedure is defined for closed contours around minima.
closed_contour_max(point p, field G, dist, maxdist, thresh) p0 = find_max_near(p, G, maxdist) set visited = {} set tovisit = {p0} while tovisit is not empty q = remove point from tovisit if (q in visited) then continue add q to visited if (gcdist(p0,q) > dist) then return false if (G[p0] -G[q] < thresh) then add all neighbors of q to tovisit return true Algorithm 6 Determine if a candidate node p satisfies the requirement that there exists another node p0 within distance dist of p with G[p] > thresh.
) and (gcdist(p,q) < dist) then add q into s remove q from P[u] p = q else if (gap < maxgap) then gap = gap + 1 else break add s into S

Thresholding
Additional threshold criteria may be applied at the detection stage in order to further eliminate undesirable candidates.For example, a common threshold criteria requires that a field G satisfy some minimum value within a distance dist of the candidate, as implemented in Algorithm 6. TempestExtremes implements thresholding via the command line argument --thresholdcmd and includes thresholds for a particular field at candidate nodes to be greater than, less than, equal to or not equal to a specified value.
if (gcdist(p0,q) > dist) then return false if (G[p0] -G[q] < thresh) then add all neighbors of q to tovisit return true Algorithm 6 Determine if a candidate node p satisfies the requirement that there exists another node p0 within distance dist of p with G[p] > thresh.The basic track stitching procedure (which represents the Reduce() stage in MapReduce) is implemented in Algorithm 7 using the output from the detection procedure at each time level (stored in set array P[1..T]).It requires additional parameters to specify a maximum great-circle distance between in-sequence nodes (dist), and a maximum gap size (maxgap).Here, gap size refers to the maximum number of sequential non-detections that can occur before a path is considered terminated.This argument is useful, for instance, for tracking tropical storms that temporarily weaken below acceptable criteria before restrengthening.node pmax = find_max_near(node p, field G, maxdist) set visited = {} set tovisit = {p} pmax = p while tovisit is not empty q = remove node from tovisit if (q in visited) then continue add q to visited if (gcdist(p,q) > maxdist) then continue if Algorithm 5 Determine if there is a closed contour in field G of magnitude thresh around the point p0, defined by p0 = find_max_near(p, G, maxdist), within distance dist.That is, along all paths away from p0, the field G must drop by at least thresh within distance dist.
The closed contour criteria is depicted in Fig. 54.An analogous procedure is defined for closed contours around minima.
.T] and maximum great-circle distance between nodes at subsequent time levels dist.
path set S = stitch_nodes(set array P[1..T], dist, maxgap) for each time level t = 1..T K[t] = build_kd_tree(P[t]) for each time level t = 1..T while P[t] is not empty initialize empty path s p = remove next candidate from P[t] add p into s gap = 0 for time level u = t+1..T q = kd_tree_nearest_neighbor(K[u], p) if (q in P[u]) and (gcdist(p,q) < dist) then add q into s remove q from P[u] p = q else if (gap < maxgap .This gives rise to a tree structure (right) that, in conjunction with an input node, can then be searched recursively for a corresponding rectangular domain in physical space.The last leaf node is labeled as the best candidate for nearest neighbor and the tree is "unwound" to test other potential candidates.The number of nodes that need to be examined is limited to domains that overlap a hypersphere with origin at the input node and with distance to the current candidate.candidates.Once paths have been constructed, additional criteria can be applied -for instance, minimum path length or additional criteria based on minimum path length or minimum distance between the start and endpoints of the path (see Appendix E).Thresholds based on field values may also be applied; e.g., wind speed must be greater than a particular value for at least eight time steps of each track.

Parallelization considerations
Feature tracking fits well into a general framework known as MapReduce (Dean and Ghemawat, 2008), which is a combination of a Map(), an embarrassingly parallel candidate identification procedure applied to individual time slices, and a Reduce(), which stitches candidates across time to build feature tracks.A key advantage of employing this framework is that substantial work has been undertaken to understand optimal strategies for parallelization of MapReduce-type algorithms (e.g., Prabhat et al., 2012)  ./DetectCyclonesUnstructured--in_data "$uvfile;$tpfile;$hfile" --out $outf --searchbymin PRMSL_L101 --mergedist 2.0 --closedcontourcmd "PRMSL_L101,200.,4,0;TMP_L100,-0.4,8.0,1.1"--outputcmd "PRMSL_L101,max,0; _VECMAG(U_GRD_L100,V_GRD_L100),max,4; HGT_L1,max,0" All outputs from DetectCyclonesUnstructured are then concatenated into a single file containing candidates at all times (pgbhnl.dcu_tc_all.dat).Candidates are then stitched in time to form paths, with a maximum distance between candidates of 8.0 • (great-circle distance), consisting of at least eight candidates per path, and with a maximum gap size of two (most consecutive time steps with no associated candidate).Because localized shallow low-pressure regions that are unrelated to tropical cyclones can form as a consequence of topographic forcing, we also require that for at least eight time steps the underlying topographic height (zs) be at most 100 m.The associated command line for StitchNodes is ./StitchNodes--in pgbhnl.dcu_tc_all.dat--out pgbhnl.dcu_tc_stitch.dat--format "i,j,lon,lat,psl,maxu,zs" --range 8.0 --minlength 8 --maxgap 2 --threshold "zs,<=,100.0,8" Once the complete set of tropical cyclone paths has been computed, total tropical cyclone counts over each 2 • grid cell are plotted in Fig. 5.The results show very good agreement with reference fields (Gray, 1968;Knapp et al., 2010).

Extratropical cyclones in CFSR
For our second example, we are interested in tracking extratropical cyclone features (defined by a cyclonic structure with no distinct warm core).The command line we have used to detect cyclonic features without the characteristic warm core of tropical cyclones (here referred to as extratropical cyclones) is given below.The command is identical to the tropical cyclone (TC) detection configuration specified in Sect.3.1, except it requires that the feature does not possess a closed contour in the 400 hPa temperature field (no warm core).

Tropical easterly waves in CFSR
Tropical easterly waves are our third example of a pointwise feature that has been assessed in the tracking literature.In this example, Northern Hemisphere easterly waves (associated with positive relative vorticity) are tracked separately from Southern Hemisphere easterly waves (associated with negative relative vorticity).DetectCyclonesUnstructured and StitchNodes are executed separately in both hemispheres and the resultant track files concatenated.All tracking is performed on the 600 hPa relative vorticity field, using relative vorticity maxima for Northern Hemisphere waves and relative vorticity minima for Southern Hemisphere waves.Since CFSR only provides absolute vorticity, relative vorticity must first be extracted by taking the difference between absolute vorticity and the planetary vorticity (the Coriolis parameter).This is done on the command line via _DIFF(ABS_V_L100,_F()), where ABS_V_L100 is the CFSR absolute vorticity variable and _F() is a builtin function for computing the Coriolis parameter (defined by f = 2 sin φ).In the Northern Hemisphere, we follow Thorncroft and Hodges (2001) and isolate tropical easterly wave features by requiring a drop of relative vorticity equal to 5 × 10 −5 s −1 .The command line used is as follows: ./DetectCyclonesUnstructured --in_data "$uvfile;$hfile" --out $outf --searchbymax "_DIFF(ABS_V_L100(0),_F())" --mergedist 2.0 --closedcontourcmd "_DIFF(ABS_V_L100(0), _F()),-5.e-5,4,0"--outputcmd "ABS_V_L100(0),max,0 ;_DIFF(ABS_V_L100(0),_F()),max,0; HGT_L1,max,0" --minlat -35.0 --maxlat 35.0 Tropical easterly wave paths are constructed using a maximum distance of 3 • great-circle distance between subsequent detections, a minimum path length equal to eight sequential detections, no allowed gaps, and a distance of at least 10 • between track start and endpoint.Northern (Southern) Hemisphere waves must also be present in the Northern (Southern) Hemisphere for at least eight time steps (2 days).The command line for Northern Hemisphere waves is as follows: ./StitchNodes --in pgbhnl.dcu_aew_nh_all.dat--out pgbhnl.dcu_aew_nh_stitch.dat--format "i,j,lon,lat,relv,zs" --range 3.0 --minlength 8 --maxgap 0 --min_endpoint_dist 10.0 --threshold "lat,<=,25.0,8;lat,>=,0.0,8" An analogous procedure is applied in the Southern Hemisphere, except it searches for minima in the relative vorticity field and limits the latitudinal range in StitchNodes to [25S,0] for at least eight time steps.Counts of total (Northern Hemisphere plus Southern Hemisphere) tropical easterly waves within each 2 • grid volume are given in Fig. 7, showing heavy wave activity throughout the Atlantic and Pacific basins.These results are very similar to other reported easterly wave densities, such as in Belanger et al. (2014) and Thorncroft and Hodges (2001), except for (a) the substantially enhanced tropical easterly wave count reported over eastern Africa (which could be eliminated by filtering over topography) and (b) essentially no observed wave activity off of the western coast of South America.Nonetheless, it is well known that easterly wave climatology varies strongly across reanalysis datasets and exhibits sensitivity to the choice of tracking scheme (Hodges et al., 2003).

Tropical cyclone forecast trajectories
For our third example, we have used TempestExtremes to track forecasted tropical cyclones in numerical weather prediction simulations.Here, we show two deterministic forecasts (initialized at 00Z on 21 and 22 October 2012) for Hurricane Sandy using the Community Atmosphere Model (CAM) with a 0.125 • (14 km) grid spacing over the North Atlantic basin.Details about the model setup, forecast skill of CAM, and a case study of Hurricane Sandy results can be found in Zarzycki and Jablonowski (2015).Both forecasts were initialized prior to the National Hurricane Center declaring Sandy as a tropical depression, which occurred at 12Z on 22 October (Blake et al., 2013).In Fig. 8, black dots indicate Sandy's forecast trajectory when applying the operational tracker used by the National Centers for Environmental Geosci.Model Dev., 10, 1069-1090, 2017 www.geosci-model-dev.net/10/1069/2017/Prediction (NCEP) (Marchok, 2002) while red dots show the same using a sample configuration of TempestExtremes.This configuration finds local minima in the 6-hourly SLP (slp) fields, which cannot lie within a great-circle distance of 10.0 • of another.An increase in SLP of at least 0.5 hPa within 5 • of the candidate node is required (closed contour) as is a decrease in 300 hPa air temperature (tm) of 0.1 K within 5 • of the node, with a 1.0 • offset permitted between the upper-level warm core maximum and sea level pressure minimum (relevant for sheared TCs where the vortex may be tilted).Note that no "best guess initial location" of the cyclone is defined, as is the case with many operational tracking systems.The tracker command line is as follows: ./DetectCyclonesUnstructured --in_data $ffile --out $outf --mergedist 10.0 --closedcontourcmd "slp,0.5,5.0,0;tm,-0.1,5.0,1.0"--outputcmd slp,min,0;_VECMAG(u850,v850), max,2;_VECMAG(u_ref,v_ref),max,2" ./StitchNodes --in cand.cyc--out forecast.traj--format "i,j,lon,lat,slp,wind850,windbot" --range 6.0 --minlength 8 --maxgap 2 The results here demonstrate good agreement with the NCEP vortex tracker, highlighting the capability of this framework to track even pre-genesis storm features, although the sensitivity (and associated potential noise) required to find weak, shallow or sheared storms depends on the thresholds defined in DetectCyclonesUnstructured.Some differences between tracked storm centers are noted, particularly at the beginning of the forecasts, where the storm's SLP is greater (weaker) than 1005 hPa.This is due to the fact that the pre-genesis vortex is naturally somewhat disorganized, and the NCEP tracker uses an average of multiple primary fixes (e.g., 700 and 850 hPa relative vorticity, sea level pressure, 700 and 850 hPa geopotential heights) to define the cyclone center, whereas this configuration of TempestExtremes defines storm location based on sea level pressure minimum only.

Tropical cyclones in a simulation with VR-CAM
For our final example, we assess the differences in tropical cyclone character obtained from native and regridded datasets.Using the variable-resolution spectral element option in CAM (VR-CAM-SE; Neale et al., 2012;Zarzycki et al., 2014b) to refine the Northern Hemisphere to 0.25 • resolution, a simulation of a hurricane season (June-November) has been performed.With the high-order spectral element dynamical core used to solve the fluid equations in the atmosphere, VR-CAM-SE has been demonstrated to be effective in simulating tropical cyclone-like features (Zarzycki and Jablonowski, 2014; www.geosci-model-dev.net/10/1069/2017/Geosci.Model Dev., 10, 1069-1090, 2017 Figure 10.Tropical cyclone trajectories and associated intensities as obtained from the simulation of a single hurricane season in CAM 3.5 using (top) native spectral element grid data and (bottom) data regridded to a regular latitude-longitude grid with 0.25 • grid spacing.Zarzycki et al., 2014a;Zarzycki and Jablonowski, 2015).
Since VR-CAM-SE uses an unstructured mesh with degrees of freedom stored at spectral element Gauss-Lobatto (GL) nodes, data are typically analyzed only after being regridded to a regular latitude-longitude mesh of approximately equal resolution.The regridding procedure has the potential to clip local extrema and smear out grid-scale features.
For this example, we use the high-order regridding package TempestRemap (Ullrich and Taylor, 2015;Ullrich et al., 2016) for remapping the native spectral element output to a regular latitude-longitude grid with 0.25 • grid spacing.For purposes of determining connectivity on the variableresolution spectral element mesh, we connect GL nodes along the coordinate axis of each quadrilateral element (see Fig. 9).DetectCyclonesUnstructured is then applied to both the native grid data and the regridded data on the regular latitude-longitude mesh (using the configuration specified in Sect.3.1) and tropical cyclones are categorized (color-coded) by maximum surface wind speed as defined by the Saffir-Simpson scale (Simpson, 1974), such that orange and red trajectories represent the strongest classifications of storms.The results of this analysis are depicted in Fig. 10.As expected, the native grid output produces essentially identical tracks, but an increase in tropical cyclone intensity in some cases (with some tropical cyclones dropping down by a full category as a consequence of the remapping procedure and discrete nature of binning storm strength).

Conclusions
Automated pointwise feature trackers have been frequently and successfully employed over the past several decades to extract useful information from large climate datasets.With spatial and temporal resolution increasing rapidly in response to enhanced computational resources, climate datasets have grown increasingly unwieldy and so there has been a growing need for such large dataset processing tools.This paper has outlined a framework for pointwise feature tracking (TempestExtremes) that exposes a suite of generalized kernels drawn from the literature on trackers of the past several decades.This framework is sufficiently robust to be applicable to many climate and weather datasets, including data on unstructured grids.We expect such a framework would be useful for isolating uncertainties that emerge from particular parameter choices in tracking schemes, or to compute optimal threshold values for detecting pointwise features in, e.g., reanalysis data.Future development plans in Tempes-tExtremes include the construction of analogous kernels for tracking areal features (blobs), such as clouds or atmospheric rivers.

Code availability
The open-source software described in this paper has been released as part of the TempestExtremes software package, and is available for use under the Lesser GNU Public License (LGPL).Software and examples can be obtained from GitHub at https://github.com/ClimateGlobalChange/tempestextremes.
Geosci.Model Dev., 10, 1069-1090, 2017 www.geosci-model-dev.net/10/1069/2017/This appendix reviews the existing literature on extratropical cyclone tracking, one of the earliest and most common instances of both manual and automated feature tracking.Manual counts of cyclones were performed by Petterssen (1956) in the Northern Hemisphere from 1899 to 1939, and latter binned by Klein (1957) to determine the spatial distribution of such storms.These techniques were later refined by Whittaker and Horn (1982) by accounting for cyclone trajectories.A similar survey in the Southern Hemisphere was performed by Taljaard (1967) for July 1957-December 1958.
Manual tracking and characterization of cyclones was also performed by Akyildiz (1985) using ECMWF forecast data for the 1981/1982 winter.One of the first automated detection and tracking for extratropical cyclones was developed by Williamson (1981) using nonlinear optimization to fit cyclonic profiles to anomalies in the 500 mb geopotential height field.Storms were then tracked over a short forecast period using the best fit to the cyclone's center point.Counts of cyclones neglecting the cyclone trajectory were automatically generated from climate model output for both hemispheres by Lambert (1988) using local minima in 1000 hPa geopotential height.This method had some shortcomings, including mischaracterization of local lows due to Gibbs' ringing and topographically driven lows.To overcome these problems, Alpert et al. (1990) proposed an additional minimum threshold on the local pressure gradient.Similarly, Le Treut and Kalnay (1990) detected cyclones in ECMWF pressure data using a local minima in the sea level pressure that must also be 4 hPa below the average sea level pressure of neighboring grid points, and must persist for three successive 6 or 12 h intervals.Murray and Simmonds (1991) extracted low pressure centers from interpolated general circulation model (GCM) data using local optimization, based on earlier work in Rice (1982).These original papers primarily sought minima in the SLP field or looked for maxima in the Laplacian of the SLP field.
Several modern extratropical cyclone detection algorithms remain in use, having built on this earlier work.Short descriptions of many of these schemes are given here.Some of these algorithms use the notion of a local neighborhood or periphery, as defined in Fig. A1.
- Serreze et al. (1993); Serreze (1995) tracked cyclones in a ∼ 381-400 km Arctic dataset.Candidates were identified using a local minimum SLP at least 2 hPa higher than immediate neighbors.Tracking was performed with a maximum search distance of 1400 km per 12 h period.
- Sinclair (1994Sinclair ( , 1997) ) tracked cyclones in a 2.5 • ECMWF dataset over the Southern Hemisphere.Candidates were identified using a local minimum in the 1000 hPa geostrophic vorticity field ζ g (computed from the Laplacian of the 1000 hPa geopotential), adjusted for topography and the presence of heat lows (see paper for details), which further satisfied ζ g < −2 × 10 −5 s −1 .
- Blender et al. (1997) tracked cyclones in T106 (∼ 125 km) ECMWF analyses.Candidates were identified using a local minimum in the 1000 hPa geopotential height field, and required to have a positive mean gradient in the 1000 hPa geopotential height field in a 1000 × 1000 km 2 area around each candidate.Tracking was performed using a nearest-neighbor search with a maximum displacement velocity of 80 km h −1 , eliminating cyclones with tracks shorter than 3 days.
- Lionello et al. (2002) tracked cyclones in a T106 (∼ 125 km) ECHAM-4 dataset.Candidates were identified using a local minimum in the SLP field.Tracking was performed using previous cyclone velocity to demarcate a prediction region, and candidates were discarded if they do not continue into the prediction region.
- Zolina and Gulev (2002) tracked cyclones in a T106 (∼ 125 km) and a T42 (∼ 300 km) dataset.Candidates were identified using a local minimum in the SLP field.
- - Pinto et al. (2005) tracked cyclones in T42 (∼ 300 km) NCEP reanalysis, regridded onto a 0.75 • grid by cubic spline interpolation.Candidates were identified using local minima in the pressure field that were within 1200 km of a maximum in the quasi-geostrophic relative vorticity, which was computed from the Laplacian of pressure.Candidates that were over a topography of above 1500 m were removed.They further required that the quasi-geostrophic relative vorticity was greater than 0.1 hPa / ( • lat) and only the strongest candidates within 3 • were retained.Cyclone tracking required a prediction velocity and search following Murray and Simmonds (1991).
- Benestad and Chen (2006) tracked cyclones in 2.5 • ERA40 data.Candidates were identified using multiple least-squares regression to estimate the values of the coefficients of a Fourier approximation to the SLP field (in effect a smoothing operator), followed by a 1-D search in the north-south and east-west directions.
- Simmonds et al. (2008) tracked cyclones in several 2.5 • datasets over the Arctic.Candidates were identified using local minima in the Laplacian of pressure, rejecting cyclones over topography above 1000 m and requiring the presence of a nearby pressure minimum.
Identified lows must satisfy a Laplacian with value > 0.2 hPa / ( • lat) 2 over a radius of 2 • .Tracking used a probability estimate using a predicted position.

Appendix B: A review of tropical cyclone tracking algorithms
More recently, and as higher-resolution climate data have become available, extratropical cyclone tracking techniques have been modified in order to support tropical cyclone tracking.To eliminate "false positives" associated with extratropical cyclones and weak cyclonic depressions, many schemes require that the candidate be associated with a nearby warm core and be associated with a minimum threshold on surface winds for at least 1-3 days.The definition of a "warm core" varies between modeling centers, including such options as air temperature anomaly on pressure surfaces (Vitart et al., 1997;Zhao et al., 2009;Murakami et al., 2012), geopotential thickness (Tsutsui and Kasahara, 1996) and decay of vorticity with height (Bengtsson et al., 2007a;Strachan et al., 2013).Additional filtering of candidate storms over topography or within a specified latitudinal range may be required.To better match observations, additional geographical-, model-or feature-dependent criteria may be applied (Camargo and Zebiak, 2002;Walsh et al., 2007;Murakami and Sugi, 2010a;Murakami et al., 2012).
It is widely acknowledged that weaker tropical storms are difficult to track, and the observational record of these lessintense, short-lived storms is questionable (Landsea et al., 2010).
A tabulated overview of the thresholds utilized by many of these schemes can be found in Walsh et al. (2007), along with several proposed guidelines on detection schemes.We extend this tabulation with the following short descriptions of many published schemes.
-Broccoli and Manabe (1990) tracked tropical cyclones in a R15 (∼ 600 km) dataset and a R30 (∼ 300 km) dataset.Candidates were identified from PSL that had a 1.5 hPa local min (R15) or 0.75 hPa local min (R30), with local surface wind velocity > 17 m s −1 , and latitude < 30 • .Tracking was performed using nearest-neighbor search with a maximum velocity of 1200 km day −1 .
- Wu and Lau (1992) tracked tropical cyclones in a 7.5 • longitude × 4.5 • latitude dataset.Candidates were identified by local minima in 1000 hPa geopotential height with a positive 950 hPa relative vorticity, negative 950 hPa divergence, positive 500 hPa vertical velocity, latitude < 40.5 • , locally maximal 200 h minus 1000 hPa layer thickness that exceeded the average layer thickness within 1500 km west to east by 60 m, and 950 hPa wind greater than 17.2 m s −1 locally.Tracking imposed a maximum tropical cyclone velocity of 7.5 • longitude or 9 • latitude per day.
- Haarsma et al. (1993) tracked tropical cyclones in a ∼ 300 km dataset.Candidates were identified from local minimum PSL, with 850 hPa relative vorticity > 3.5 × 10 −5 s −1 .Temperature anomaly T was required to satisfy T 250 > 0.5 K at 250 hPa, T 500 > −0.5 K at 500 hPa, and T 250 − T 850 > −1.0 K, where the anomaly is computed against a 15 • × 15 • spatial mean around the center of the storm.Tracking required tropical cyclones to be persistent for a minimum of 3 days.
- Tsutsui and Kasahara (1996) tracked tropical cyclones in a T42 (2.8 • , ∼ 300 km) dataset.Candidates were identified as minima in the 1000 hPa geopotential height field, with at least an average drop of 20 m among neighboring points, and a further 20 m drop of average among neighboring points from periphery.Candidates were further required to satisfy that the average local 900 hPa vorticity was cyclonic, average local 900 hPa divergence was negative, average local 500 hPa vertical velocity was upward, 200 hPa minus 1000 hPa layer thickness maximum among neighbors was greater than any value in periphery, and average local 200 hPa zonal wind velocity was less than 10 m s −1 or local points contained at least one point with easterly velocity.The latitude of the candidate was require to be less than 40 • , graphic height underlying candidates was less than 400 m, one local point had a 900 hPa wind speed of at least 17.2 m s −1 , and one local point exceeded 100 mm d −1 over at least 1 day.Tracking required tropical cyclones to be persistent for a minimum of 2 days.
- Vitart et al. (1997Vitart et al. ( , 1999Vitart et al. ( , 2001Vitart et al. ( , 2003) ) tracked tropical cyclones in a T42 (2.8 • , ∼ 300 km) dataset.Candidates were identified as 850 hPa relative vorticity maxima > 3.5 × 10 −5 s −1 with a nearby PSL minimum.They were required to possess a warm core within 2 • latitude, defined as a local average 500 hPa to 200 hPa temperature maximum with a decrease of 0.5 K in all directions within 8 • , and a local maximum in the 200-1000 hPa layer thickness with a decrease of 50 m in all directions within 8 • .Tracking required the minimum distance between storms to be 800 km day −1 , that tropical cyclones lasted at least 2 days and that the maximum wind velocity within 8 • of the storm center must be 17 m s −1 for at least 2 (not necessarily consecutive) days.
-Walsh (1997); Walsh and Watterson (1997); Walsh and Katzfey (2000) tracked tropical cyclones in a 125 km regional climate dataset over Australia.Candidates were identified as points with 850 hPa relative vorticity > 2.0 × 10 −5 s −1 , temperature anomaly sum T 700 + T 500 + T 300 > 0 K and T 300 > T 850, with anomaly computed against the mean over a region 2 grid points north-south and 13 grid points east-west.Candidates were also required to have 10 m surface wind > 10 m s −1 and 850 hPa tangential wind speed > 300 hPa tangential wind speed.Tracking required tropical cyclones to be persistent for a minimum of 2 days.
- Krishnamurti et al. (1998) tracked tropical cyclones in a T42 (∼ 300 km) climate dataset.Their approach was similar to Bengtsson et al. (1995Bengtsson et al. ( , 1996)), except using a 4 × 4 grid-point region for the 850 hPa wind maximum, the SLP minimum and the temperature mean.Tracking required tropical cyclones to be persistent for a minimum of 1 day.
-Nguyen and Walsh ( 2001) assessed a 125 km regional dataset over Australia using a similar approach to Walsh and Watterson (1997).The vorticity requirement was changed to 850 hPa relative vorticity > 1.0 × 10 −5 s −1 with a PSL minimum within 250 km.Candidates also must possess a mean wind speed in a 500 km × 500 km region at 850 hPa that was larger than mean wind speed at 300 hPa, and a mean tangental wind speed within a radius of 1 • and 2.5 • greater than 5 m s −1 .Tracking required tropical cyclones to be persistent for a minimum of 1 day, with relaxed criteria after this time (see paper for further information).
- Sugi et al. (2002) tracked tropical cyclones in a T106 (∼ 125 km) climate dataset, using tracking criteria similar to Bengtsson et al. (1995).Candidates were identified by local PSL minima that was at least < 1020 hPa.
Tracking required tropical cyclones to be persistent for a minimum of 2 days.
- Camargo and Zebiak (2002) tracked tropical cyclones in a T42 (∼ 300 km) climate dataset.Their approach was similar to Bengtsson et al. (1995Bengtsson et al. ( , 1996)), except with basin-specific thresholds are applied for 850 hPa relative vorticity, 850 hPa wind speed, and temperature anomaly sum T 700 + T 500 + T 300.Thresholds were determined by sampling the tails of probability density functions for relevant variables in each ocean basin.Following candidate identification, a relaxed 850 hPa relative vorticity threshold (> 1.5 × 10 −5 s −1 ) in an area of 3 × 3 grid points around prior detections was applied to construct trajectories.Tracking required tropical cyclones to be persistent for a minimum of 2 (1.5) days in daily (6-hourly) output.
- Tsutsui (2002) tracked tropical cyclones in a T42 (2.8 • , ∼ 300 km) dataset.Their approach was similar to Tsutsui and Kasahara (1996), but with simplified criteria.Candidate PSL was required to be less than the local average minus 2 hPa, and local average PSL must be less than the periphery average minus 2 hPa.Layer thickness between 200 and 700 hPa, denoted by Z, was required to satisfy Z 0 + max(Z ±1 ) > 2max(Z ±2 ), where Z ±1 denotes immediate neighbors and Z ±2 denotes the periphery.
-Cheung and Elsberry (2002); Halperin et al. (2013) tracked tropical cyclones in weather forecast models using an approach similar to Walsh (1997).Candidates required a grid-point maximum in 850 hPa relative vorticity larger than all surrounding grid points within 4  Halperin et al. (2013), tropical cyclones were required to persist for at least 24 h and have a 925 hPa horizontal wind speed greater than a model-specific threshold within 5 • of the PSL center.
- Walsh et al. (2004) tracked tropical cyclones in a 30 km dataset using a similar tracking strategy to Nguyen and Walsh (2001).Candidates' temperature anomaly was computed against a 1200 km longitude × 400 km latitude region, and the mean wind speed was computed over a 800 km × 800 km region around the storm.Candidates were further required to have local 10 m meridional velocity ≥ 17 m s −1 .
- McDonald et al. (2005) tracked tropical cyclones in a 2.5 • latitude by 3.75 • longitude dataset.Candidates were identified as local maxima in the 850 hPa relative vorticity field with magnitude greater than 5×10 −5 s −1 , with temperature anomaly T 300 > 0 along the track, T 300 > 0.5 K for any two points along the track and T 300 > T 850 for any two points along the track, where the anomaly was computed against a 15 • × 15 • mean.Tracking required that the storm's initial latitude was < 30 • and the tropical cyclones persisted for a minimum of 2 days.
- Chauvin et al. (2006) tracked tropical cyclones in a T319 (∼ 50 km) climate time-slice simulation.Candidates were identified by local minimum PSL with 850 hPa relative vorticity > 1.4 × 10 −4 s −1 , 850 hPa wind > 15 m s −1 , mean 700-300 hPa temperature anomaly T 700 − T 300 > 3K, T 300 > T 850, and 850 hPa wind > 300 hPa wind.Anomalies were computed against environmental values 500 km from the cyclone center.A similar relaxation technique to Camargo and Zebiak (2002) was applied to eliminate split trajectories.Tracking required tropical cyclones to be persistent for a minimum of 1.5 days.
- Oouchi et al. (2006) tracked tropical cyclones in a 20 km dataset using a similar technique to Sugi et al. (2002).
Candidate PSL was required to be 2 hPa lower than mean over 7 × 7 grid box and require a storm center of latitude < 45 • with an initial position of < 30 • .Near the candidate it was further required that the relative vorticity at 850 hPa was > 3.5 × 10 −5 s −1 , the maximum wind speed at 850 hPa was > 15 m s −1 , the maximum wind speed at 850 hPa was larger than at 300 hPa, and the temperature anomaly sum satisfied T 700 + T 500 + T 300 > 2 K near the candidate.Tracking required tropical cyclones to be persistent for a minimum of 1.5 days.
Tracking required tropical cyclones to be persistent for a minimum of 1 day.
- Knutson et al. (2007); Zhao et al. (2009) tracked tropical cyclones in a ∼ 50 km dataset using a technique similar to Vitart et al. (1997Vitart et al. ( , 2003)).Candidates were required to have an absolute 850 hPa relative vorticity maxima > 1.6×10 −4 s −1 within a 6 • ×6 • area, a local minimum in SLP within 2 • of the detection, and a maximum in average 300 hPa and 500 hPa layer temperature within 2 • that was 1 K warmer than the local mean.Tracking required storms to be persistent for a minimum of 3 days, with a maximum search radius of 400 km per 6 h, and at least 3 days with a maximum surface wind speed greater than 17 m s −1 .
-Murakami and Sugi (2010b) tracked tropical cyclones in four datasets with resolutions from TL95 (∼ 180 km) to TL959 (∼ 20 km) using a procedure similar to Oouchi et al. (2006) with a resolution-dependent relative vorticity criteria.
- Caron et al. (2011Caron et al. ( , 2013) ) tracked tropical cyclones in a 0.3 • (∼ 35 km) climate model.Candidates were required to have a local SLP minimum (with candidates merged within 2 • ), 850 hPa relative vorticity > 4.0 × 10 −5 s −1 , temperature anomalies T 500 > 1 K and T 250 > 0 K (calculated relative to 5 • radial mean around the TC center), 850 hPa relative vorticity > 250 hPa relative vorticity, and a resolution-specific surface wind threshold as in Walsh et al. (2007).Tracking required storms to be persistent for a minimum of 1 day.Tracks with a relaxed set of thresholds were calculated in parallel and applied to the main tracking to minimize broken trajectories, similar to Camargo and Zebiak (2002).
- Murakami et al. (2012) tracked tropical cyclones in four datasets with resolutions from 20 to 60 km using a procedure similar to Oouchi et al. (2006) with a resolutiondependent relative vorticity and temperature anomaly criteria.Temperature anomalies were computed against a 10  (Walsh et al., 2015) (∼ 60 km to ∼ 110 km).Tracking was performed similar to Walsh et al. (2004), except a resolution-dependent value for surface winds was applied based on Walsh et al. (2007).Tracking required storms to originate equatorward of extratropical ridges.
-Zarzycki and Jablonowski (2014) tracked tropical cyclones in a ∼ 28 km dataset.Candidates were selected using absolute 850 hPa relative vorticity maxima > 1.0 × 10 −4 s −1 with latitude < 45 • and SLP minimum within 4 • .Candidates were further required to have a local maximum 500-200 hPa average temperature within 2 • of the storm center which decreases by at least 0.8 K at a radius of 5 • in all directions.Tracking required storms to be persistent for a minimum of 2 days, with a maximum velocity of 400 km over 6 h, and required that the maximum surface wind speed within 4 • of the candidate was greater than 17 m s −1 for at least 2 days.
Tracking also allowed a maximum gap size of 1 (a single time-step failure).
- Reed and Chavas (2015) tracked tropical cyclones on a planet in radiative-convective equilibrium.Candidates were identified using SLP minima followed by a closed contour criteria that requires a pressure increase of at least 4 hPa in all directions within 5 • (great-circle distance), and using an early release of TempestExtremes.

Appendix D: Software documentation: DetectCyclonesUnstructured
This section contains the software documentation for the executable DetectCyclonesUnstructured from the TempestExtremes package.The software is provided for use within a command-line shell, such as bash.--in_data <string> is a list of input data files in NetCDF format, separated by semicolons.
--in_data_list <string> is a file containing the --in_data argument for a sequence of processing operations (one per line).
--in_connect <string> is a connectivity file, which uses a vertex list to describe the graph structure of the input grid.This parameter is not required if the data are on a latitude-longitude grid.
--out <string> is the output file containing the filtered list of candidates in plain text format.
--out_file_list <string> is a file containing the --out argument for a sequence of processing operations (one per line).
--searchbymin <string> is the input variable to use for initially selecting candidate points (defined as local minima).By default, this is "PSL", representing detection of surface pressure minima.Only one of searchbymin and searchbymax may be set.
--searchbymax <string> is the input variable to use for initially selecting candidate points (defined as local maxima).Only one of searchbymin and searchbymax may be set.
--minlon <double> is the minimum longitude for candidate points.
--maxlon <double> is the maximum longitude for candidate points.
--minlat <double> is the minimum latitude for candidate points.
--maxlat <double> is the maximum latitude for candidate points.
--minabslat <double> is the minimum absolute latitude for candidate points.
--mergedist <double> merges candidate points with distance (in degrees) shorter than the specified value.Among two candidates within the merge distance, only the candidate with lowest searchbymin or highest searchbymax value will be retained.
--closedcontourcmd <cmd1>;<cmd2>;... Eliminate candidates if they do not have a closed contour.Closed contour commands are separated by a semicolon.Each closed contour command takes the form var,delta,dist,minmaxdist. These arguments are as follows.
var <variable> is the variable used for the contour search.
dist <double> is the great-circle distance (in degrees) from the pivot within which the closed contour criteria must be satisfied.
delta <double> is the amount by which the field must change from the pivot value.If positive (negative) the field must increase (decrease) by this value along the contour.minmaxdist <double> is the distance away from the candidate to search for the minima/maxima.If delta is positive (negative), the pivot is a local minimum (maximum).
--thresholdcmd <cmd1>;<cmd2>;... eliminates candidates that do not satisfy a threshold criteria (there must exist a point within a given distance of the candidate that satisfies a given equality or inequality).Threshold commands are separated by a semicolon.Each threshold command takes the form var,op,value,dist.These arguments are as follows.
var <variable> is the variable used for the contour search.
value <double> is the value on the right-hand side (RHS) of the comparison.
dist <double> is the great-circle distance away from the candidate to search for a point that satisfies the threshold (in degrees).
op <string> is an operator that is applied over all points within the specified distance of the candidate (options include max, min, avg, maxdist, mindist).
dist <double> is the great-circle distance away from the candidate wherein the operator is applied (in degrees).
--regional is used when a latitude-longitude grid is employed, to not assume longitudinal boundaries to be periodic.
--out_header outputs a header describing the columns of the data file.

D1 Variable specification
Quantities of type <variable> include both NetCDF variables in the input file (for example, "Z850") and simple operations performed on those variables.By default, it is assumed that NetCDF variables are specified in the .ncfile as float Z850(time, lat, lon) or float Z850(time, ncol) for structured latitude-longitude grids and unstructured grids, respectively.If variables have no time variable, they have the related specification float Z850(lat, lon) or float Z850(ncol).If variables include an additional dimension, for instance, float Z(time, lev, lat, lon) or float Z(time, lev, ncol) they may be specified on the command line as Z(<lev>), where the integer index <lev> corresponds to the first dimension (or the dimension after time, if present).

D2 MPI support
The DetectCyclonesUnstructured executable supports parallelization via MPI when the --in_data_list argument is specified.When enabled, the parallelization procedure simply distributes the processing operations evenly among available MPI threads.
--min_path_dist <double> is the minimum path length, defined as the sum of all great-circle distances between candidate nodes (in degrees).
--maxgap <integer> is the largest gap (missing candidate nodes) along the path (in discrete time points).
--threshold <cmd1>;<cmd2>;... eliminates paths that do not satisfy a threshold criteria (a specified number of candidates along path must satisfy an equality or inequality).Threshold commands are separated by a semicolon.Each threshold command takes the form col,op,value,count.These arguments are as follows.
col <integer> is the column in the input file to use in the threshold criteria.
value <double> is the value on the right-hand side of the operator.count <integer> is the minimum number of candidates along the path that must satisfy this criteria.

Figure 1 .
Figure 1.An example node graph describing an unstructured grid (blue lines), where nodes are co-located with volume center-point locations (solid circles) and edges connect adjacent volumes.
k-d trees are a data structure that enable O(log N ) average time for nearest neighbor search, where N is the total number of nodes, while also requiring only O(N log N ) construction time and a O(N ) data storage requirement.Although k-d trees use 3-D straight-line distance

Figure 2 .
Figure 2. Grid cells on a latitude-longitude mesh whose centroids are within 30 • of a cell at 68 • N latitude using latitude-longitude distance (left) and great-circle distance (right).
Determine all feature paths S, given array of candidate nodes P[1..T] and maximum great-circle distance between nodes at subsequent time levels dist.path set S = stitch_nodes(set array P[1..T], dist, maxgap) for each time level t = 1..T K[t] = build_kd_tree(P[t]) for each time level t = 1..T while P[t] is not empty initialize empty path s p = remove next candidate from P[t] add p into s gap = 0 for time level u = t+1..T

Figure 51 .
Figure 51.An example node graph describing an unstructured grid (blue lines), where nodes are co-located with volume centerpoint locations (solid circles) and edges connect adjacent volumes.
Determine all feature paths S, given array of candidate nodes P[1..T] and maximum great-circle distance between nodes at subsequent time levels dist.path set S = stitch_nodes(set array P[1..T], dist, maxgap) for each time level t = 1..T K[t] = build_kd_tree(P[t]) for each time level t = 1..T while P[t] is not empty initialize empty path s p = remove next candidate from P[t] add p into s gap = 0 for time level u = t+1..T

Figure 51 .
Figure 51.An example node graph describing an unstructured grid (blue lines), where nodes are co-located with volume centerpoint locations (solid circles) and edges connect adjacent volumes.
Figure (blue l locatio 18 P. A. Ullrich and C. M. Zarzycki: A Framework for Scale-Insensitive Pointwise Feature Tracking Algorithm 4 Find the node pmax containing the maximal value of the field G within a distance maxdist of the node p.An analogous procedure find_min_near is provided for locating nodes containing minimal values of the field.

Figure 3 .
Figure3.An example two-dimensional k-d tree (k = 2) built from nodes a through h.Dividing planes are constructed by cycling through each coordinate and determining the median node (left).This gives rise to a tree structure (right) that, in conjunction with an input node, can then be searched recursively for a corresponding rectangular domain in physical space.The last leaf node is labeled as the best candidate for nearest neighbor and the tree is "unwound" to test other potential candidates.The number of nodes that need to be examined is limited to domains that overlap a hypersphere with origin at the input node and with distance to the current candidate.

Figure 4 .
Figure 4.An illustration of the closed contour criteria.Nodes shaded in white (gray) satisfy (do not satisfy) the threshold of the field value at p0.Since only edge neighbors are included, B constitutes a boundary to the interior of the closed contour.Because A lays outside the solid circle, the contour with distance d 0 is not a closed contour, whereas the dashed contour with distance d 1 does satisfy the closed contour criteria.

Figure 8 .
Figure 8. Forecast CAM trajectories for Hurricane Sandy initialized at 00Z on (a) 21 and (b) 22 October 2012.Black dots indicate trajectories defined using the NCEP operational vortex tracker with red dots denoting trajectories defined using a sample configuration of TempestExtremes.

Figure 9 .
Figure 9.An illustration of how connectivity is defined in this work for nodes on a spectral element mesh.Arrows indicate connectivity for nodes A and B.

Figure A1 .
Figure A1.The local neighborhood of a central node (shaded) typically refers to the surrounding eight nodes (diagonal hatching).The periphery, used byTsutsui and Kasahara (1996), refers to the set of nodes that surround the local neighborhood (unshaded nodes).
Gulev, S. K.: Improving the accuracy of mapping cyclone numbers and frequencies,Mon.Weather Rev., 130, . Model Dev., 10, 1069-1090, 2017 www.geosci-model-dev.net/10/1069/2017/P. A. Ullrich and C. M. Zarzycki: A Framework for Scale-Insensitive Pointwise Feature Tracking Algorithm 4 Find the node pmax containing the maximal value of the field G within a distance maxdist of the node p.An analogous procedure find_min_near is provided for locating nodes containing minimal values of the field.
Our first example employs TempestExtremes for tropical cyclones (defined here as a cyclonic structure with a distinct warm core).The command line we use to detect tropical cyclone-like features in CFSR is provided below.Climate data are drawn from three files denoted $uvfile (contain- (Saha et al., 2010)e bottlenecks associated with I/O and load balancing.TempestExtremes currently implements a simple parallelization strategy via MPI, although future work on this issue is forthcoming.As a timing example, TempestExtremes with MPI (16 tasks) finds and tracks tropical cyclones in 10 years of 6-hourly climate data on a 0.5 • latitude-longitude grid in an average of 3.8 min on the National Center for Atmospheric Research's (NCAR's) Yellowstone supercomputer.3SelectedexamplesSeveralselectedexamples are now provided.The first three examples use data from the NCEP Climate Forecast System Reanalysis (CFSR), available at 0.5 • global resolution with 6-hourly output from 1979 to the present(Saha et al., 2010).The remaining example uses a custom variable-resolution simulation (Zarzycki and Jablonowski, 2014) (6-hourly output on a 110 km base domain that is refined to 28 km in the northern Atlantic and Pacific ocean basins) on both the native grid data and the regridded latitude-longitude grid data.Geosci.Model Dev., 10, 1069-1090, 2017www.geosci-model-dev.net/10/1069/2017/3.1 Tropical cyclones in CFSR • of the candidate with maximum air temperature.Since CFSR is on a structured latitude-longitude grid, the output format is i,j,lon,lat,psl,maxu,zs, where i,j are the longitude and latitude coordinates within the dataset; lon,lat are the actual longitude and latitude of the candidate; psl is the SLP at the candidate point (equal to the maximum SLP within 0 • of the candidate); maxu is the vector magnitude of the maximum 850 hPa wind within 4 • of the candidate; and zs is the topographic height at the candidate point.