CycloneDetector (v1.0) - Algorithm for detecting cyclone and anticyclone centers from mean sea level pressure layer

. Automatic methods for identifying and tracking cyclones were firstly constructed in 1990’s and since then there was a big increase in a precision and probability of detection. These methods have been traditionally focused on cyclones (and particularly on tropical cyclones), but the question of anticyclone centers detection remained unsolved since they are usually not a source of turbulent weather, precipitation etc. However, this issue can be important in the era of the climate change. In this paper, an algorithm for an automatic detection of both, cyclones and anticyclones based on mean sea level 5 pressure field, is presented. The algorithm uses two-dimensional raster data as an input and returns a list of detected pressure systems. The main advantages of our solution are easy implementation since it is based on the standard image processing algorithm, sufficient performance of the algorithm, and especially the possibility of high-pressure systems detection. Moreover, the presented solution does not need a direct terrain filtering needed for some algorithms to be done. To validate the quality of detection algorithm results, a comparison against manually prepared data by Met Office was used. It follows from the 10 comparison that the presented algorithm produces results similar to those by Met Office . The most significant differences can be found in the detection of cyclones at the beginning or the end of the lifespan stage. Met Office detects more cyclones in these stages than the presented solution.

In this paper, we focus on the surface pressure level field. As will be explained and discussed in details further, the main features of our proposed solution are: automatic detection of low-and high-pressure systems low number and predictable behavior of input parameters fast and easy to implement solution In Chapter 2 we provide a detailed explanation of our solution. Chapter 3 presents the algorithm results. Our conclusion and final discussion can be found in Chapter 4. The proposed algorithm detects low-and high-pressure systems directly from the mean sea level pressure layer (MSLP). This quantity is widely available as an output from numerical weather prediction models (NWP), such as ECMWF, GFS, ICON, and HRRR. The algorithm uses two-dimensional raster data as an input. It returns a list of detected pressure systems. Generally speaking, there is a need for 2D Earth projection to be used for the raster data. The projection provided directly by the generated 105 NWP is used in our solution because the algorithm does not depend on any particular projection directly. Typically, inputs in the equirectangular projection are used. Nevertheless, the Lambert-Conic projection for data from regional models such as HRRR can also be applied.
The proposed solution consists of two steps: the detection of candidates (see Subsection 2.1) and the pressure system center detection within these candidates (see Subsection 2.2).

Candidates detection
In the first step of the algorithm, the isobars with a given step, C, between pressures are found. The size of the step C can be set by the user. Most of the low-pressure systems are in an intensity range 1 − 15 hP a (see Rudeva, 2008). Therefore, the isobars are usually drawn with steps of 1, 2, 4, 5, or 10 hP a. We use the same steps for low-and high-pressure systems. With lower values, we can detect more systems. Our recommendation for global models is to use 2 or 4 hP a. For local models, capturing 115 a smaller area, we recommend to use 1 hP a. From the field of isobars, we select a list of candidates, i.e., the potential centers of pressure systems.
Detecting isobars from raster data of size W ×H is a straightforward process. We use the simple Marching Squares (Lorensen and Cline, 1987) algorithm. It is a 2D version of the well-known Marching Cubes algorithm for iso-surface extraction. The detected isobars form enclosed polylines and have a fine, pixel-based, resolution. For further computations, these polylines can 120 be simplified (e.g., using Douglas and Peucker, 1973) to improve the performance of the following steps.
From the generated isobars, we select the pressure system candidates. They do not have any other closed isobars. An example can be seen in Fig. 1. There are two main candidates (#1 and #2) and several small areas that are clearly not supposed to be indicated as a pressure system centers. Therefore, these small areas need to be removed. To remove unwanted areas, we use filtering based on the area of the contour area. For each candidate, we calculate its axis-125 aligned bounding box (AABB) that is used as a rough approximation of the contour area. However, the area size cannot be easily represented in image pixels because of the projection distortions that can be visualized with Tissot's indicatrix (Goldberg and Gott, 2007). Based on this observation, we use the area size in km 2 . For this, the projection of the input image must be known and the inverse projection formulas are used to convert pixels positions [x, y] to GPS coordinates [lon, lat]. From the GPS coordinates of AABB corners, the area is approximated with equations based on Chamberlain and Duquette (2007).

130
For filtering, the area size threshold, T , is used. The threshold value can be adjusted by the user. The value depends on the average area that the pressure system occupies. We recommend to use the values between 10, 000km 2 and 20, 000km 2 .
If the area size is smaller than the selected threshold, we have two possibilities: remove the area.
change the candidate to the nearest larger detected contour (called the "parent").

135
The second option is done if the "parent" area is smaller than the M multiple of the threshold T (based on our experiments, we selected M = 100).
The process is repeated until there is a change in the candidates. The result of this filtration step, applied to Fig. 1, can be seen in Fig. 2. The small areas are entirely removed. Candidate #2 has changed from the inner contour to its "parent". After this change, small areas with significant pressure changes are preserved that would otherwise be removed.  Candidates for isobars (red thick lines) after filtration. Small areas were either removed, or the candidate was moved to the "parent" contour (compare candidate #2 in Fig. 1 and in Fig. 2). Data are from a segment over Europe taken from ICON-EU model, 2021-10-24 12:00 UTC.

Center detection
The detection of pressure system centers is based on finding local extrema in input raster data. Only discrete data are available for our purposes. A possible solution would be to create an interpolant and search for extrema directly from the interpolated function analytically. However, this is a slow and complex task. Instead, we apply convolution-based operations directly on discrete data.

145
Using the Sobel operator (see Kanopoulos et al., 1988), we calculate the first derivative dx and dy of the input raster data.
To identify the extrema, we only need to find the signs of the derivatives (signDx calculated from dx and signDy calculated form dy).

Low-pressure system
To determine if a pixel represents a low-pressure system, we use its neighborhood and compare every pixel against the mask 150 with the radius r (therefore, the width and height are 2r+1). For simplicity, we use square masks with patterns shown in Table 1.
The size of the radius r is an input parameter that can be adjusted by the user. However, based on our experiments, the default size of the radius, r, can be auto-calculated. The size is the average square root of the pixel size from the pressure system candidates AABB's (see Subsection 2.1). For high-resolution inputs, this can lead to large radius values and consequently to longer computations. If this is an issue, the radius can be set manually. Usually, based on our experiments, the size of 155 radius selected from range (5, 50) pixels is sufficient based on the input resolution. For a higher resolution, a larger radius is recommended.
For every pixel [x, y] that is inside the candidate obtained in Subsection 2.1, we apply a mask to signDx and signDy. If the signum value corresponds to the mask value, the pixel is marked as correct (the counter okX and/or okY is increased). Finally, the ratio of correctly marked pixels okX and okY against the area size of the candidate is calculated. If the ratio of correct 160 pixels is above 60%, the pixel [x, y] is marked as extreme. The process can be seen in Algorithm 1.
Algorithm 1 Pseudo-code for testing whether the pixel (x,y) is an extrema candidate. Variable r is the mask radius. The code does not handle image borders 1: maskArea ← (2 * r + 1) * (2 * r + 1) Running the algorithm 1, extrema areas within candidate contours are obtained as can be seen in Fig. 3. Sometimes, there can be multiple isolated areas detected inside a single candidate. This is caused by numerical problems.
To overcome this, we select the largest area from all isolated areas inside a single candidate and find its centroid. This centroid is used as the location of the pressure system. 165

High-pressure system
High-pressure systems are detected with the same algorithm as low-pressure systems in 2.2.1. The only change is in the masks.
We need to find local maxima and. Therefore, we use a mask with swapped signs (cf., Table 2). Table 2. Example masks for local maximum detection. The left mask for dx, the right mask for dy. The mask radius r is set to be 2.

Filtering results
In some cases, it may happen that low-and high-pressure systems are detected together inside a single candidate. This is often caused by multiple shallow systems inside a single candidate. In this case the algorithm is set to select a larger system. We can also use weighting based on the distance of the centers of the systems from the center position of the candidate contour.
Centers that are too close to each other might cause another problem. If this is a problem, we may apply the optional, userdefined threshold distance D in km. In this case, if two systems are closer than the threshold D, the system with a smaller area 175 is removed.

Performance
The performance of the algorithm depends on the resolution of input raster data and on the number of detected pressure systems. The more systems have been detected, the slower the algorithm is.
Algorithm performance can be influenced by detection if a pixel is inside a candidate (see Subsection 2.2.1). In the naive 180 solution, we have to test if a pixel is inside the polygon for every candidate. This is rather slow. However, if a bounding volume hierarchy created from candidates of AABBs is used (e.g., Pharr et al., 2016), we can improve the speed of the algorithm considerably. In this case, the point in AABB is tested first. Therefore, a computationally expensive test to determine whether a pixel is inside the polygon is performed only for several contours.

185
The proposed solution has its limitations in the input projection. It follows from the projection of the input raster data, that not all the systems may be detected. For example, in the case of equirectangular projection, systems near the poles may be detected incorrectly. Hence, instead of one center, several or no centers may be obtained. This problem is caused by a distortion near the poles. Data reprojection might lead to a solution if one projection for low-distorted latitudes is used, which may be from the interval (−85 • , 85 • ). A different projection is only applied for polar regions.

190
There might be another minor limitation in the selection of the input parameters, mainly the area size threshold T . However, once the parameters have been set, they can be reused for the same NWP model in most cases.

Experiments and results
In this section, a comparison of the proposed solution with a manually created list of pressure systems is provided.

Algorithm setup 195
To evaluate the proposed solution, we have used historical data obtained from the reanalysis of the ECMWF model (dataset ERA5, Hersbach et al., 2021). We have run the proposed algorithm with the following settings as shown in Table 3.  The relation between input parameters, area threshold T and isobars step C, is depicted in Fig. 4  We recommend to use auto-calculated values of mask radius r from Table 3. For comparison, we present some autocalculated values in Table 4. The values are obtained from multiple runs over different times. Results are averaged and rounded to whole numbers to represent pixels.