Algorithm for detecting cyclone and anticyclone centres from mean sea level pressure layer

ABSTRACT

Automatic methods for identifying and tracking pressure systems have been traditionally focused on cyclones (and particularly on tropical cyclones), but the question of anticyclone centre detection remained unsolved since they are usually not a source of turbulent weather, precipitation, etc. An algorithm for automatic detection of both, cyclones and anticyclones based on the mean sea-level pressure field, is presented. The main advantages of our solution are easy implementation based on two-dimensional raster data, sufficient performance of the algorithm, and especially the possibility of high-pressure systems detection. Moreover, the presented solution does not need direct terrain filtering.

KEYWORDS:

• Cyclone detection
• anticyclone detection
• autodetection
• image processing

1. Introduction

Pressure systems are the main features of atmospheric circulation (Ahrens and Henson 2018) . Cyclones and anticyclones (or lows and highs) affect daily weather over the whole Earth, and they are the main drivers of surface wind. Cyclones (anticyclones) are areas of relatively low (high) pressure compared to its surrounding, usually having enclosed isobars (on the surface pressure field) or isohypse (at the higher pressure level), they are, based on convention, divisible by $5$ $hPa$, or $4$ $gpdam$ (geopotential decametres). Cyclones are essentially important for transporting heat and energy between tropical and polar regions. Therefore, it is useful to have a profound knowledge of cyclones and anticyclones: their behaviour, occurrence and movement. Especially in the era of deepening climate change, it is important to study changes in their frequency in different parts of the world.

For this reason, we need to have a robust algorithm for detecting and tracking pressure systems. While we can find many papers on automatic tropical cyclone position detection, papers on the detection of extra-tropical cyclones are rare, especially those elaborating on methods used in operational praxis. The reason is caused by an urgent need to find the tropical cyclone centre connected with very pronounced weather and its changes. Last but not least, tropical cyclones can be detected more easily due to their regular structure with a quasi-circular closed centre compared to extra-tropical cyclones (e.g. Wang et al. 2020), but as our paper is focused on extra-tropical systems, tropical cyclones will not be discussed in detail.

Some decades ago, the detection of pressure systems had traditionally been based on a manual analysis of synoptic weather charts. Although this method is accurate, it is connected with a huge amount of time-consuming manual work. It has to be pointed out that in data-scarce areas (e.g. over oceans or desert regions), extrapolation does not necessarily lead to accurate detection of pressure system centres. This time-consuming method could also have another source of problems, e.g. inconsistencies that may have arisen in the identification of cyclone and anticyclone centres by various analysts from day to day. Later, with the onset of enhanced numerical models and better observational techniques, this lack-of-data disadvantage has been minimised, but the need for automated detection of low and high centres has appeared. Automated schemes started to be used from the 1990s and these schemes have overcome the above described drawbacks of manual analyses.

One of the first automatic schemes was introduced by Murray and Simmonds (1991) with a focus on mid-latitude cyclones in the Southern Hemisphere. Other automatic methods for identifying and tracking cyclones were constructed during the end of the 20th and the beginning of the 21^st century (e.g. Hodges 1994, Blender and Schubert 2000, Rudeva and Gulev 2007, Raible et al. 2008, Hanley and Caballero 2012, Lu 2017). These methods could benefit from new observation technologies as well as from the development of numerical weather models.

To detect cyclones, various meteorological parameters and criteria were used in the past, typically, pressure or vorticity parameters. Searching minima in the sea-level pressure and/or in the geopotential height at pressure levels can be found in more studies (e.g. Geng and Sugi 2001, Rudeva and Gulev 2007, Hanley and Caballero 2012). In addition, some papers focus on finding minima of a geopotential height of the $850$ $hPa$ level (Kew et al. 2010), using relative vorticity of $850$ $hPa$ at pressure levels (Flaounas et al. 2014), or a combination of these parameters (e.g. mean sea-level pressure minima and low-level vorticity maxima (Simmonds et al. 2008)). A very comprehensive overview comparing extra-tropical cyclone detection and tracking algorithms before 2013 can be found in the work of Neu et al. (2013). This overview poses many interesting questions and highlights problems that may have arisen when constructing the detection algorithms, e.g. detecting cyclones over elevated topography that has not been standardised. Various methods deal with this point differently. Another approach for detecting cyclone centres is represented by a deep learning approach, first presented by Hinton et al. (2006). This approach can be based either on classical convolutional neural networks alone or in combination with a region proposal network.

Solution based on convolutional neural network was presented by Matsuoka et al. (2018). They trained the network to track tropical cyclones and detect the centres.

Recently, Lu et al. (2020) presented a mask region-based convolutional neural network for detecting cyclones showing higher efficiency, especially for recognising shallow or moderately deep cyclones of subsynoptic scale.

Generally, the scheme of cyclone centre detection comprises several steps. Firstly, the centre of the anti/cyclone has to be detected. This can be done by several methods, such as finding the local pressure maximum or minimum or searching the Laplacian of the pressure or relative vorticity. The detection of the anti/cyclonic centre can be done by various methods. For example, Lu (2017) searched for a point lower than the surrounding grid points (in all directions) on a latitude-longitude mesh in the grid of a mean sea-level pressure or a geopotential height of the $850$ $hPa$ level.

Ullrich and Zarzycki (2017) proposed a solution based on multiple inputs (pressure, wind, and topographic height layers). They use contours to detect systems based on user-specified thresholds. However, their system is more suitable for tropical cyclones.

Jiang et al. (2020) introduced an eight-section slope detecting method, based on mean sea-level pressure or geopotential height. Usually, several iterations are required over an enclosed contour search stage to conduct the procedure for every cyclone-centre candidate. After the centre is detected, its pressure value and other parameters like the cyclone’s depth, size, or radius can be obtained.

Finally, it is important to know what is the cyclone centre detection accuracy. Generally speaking, no identification algorithm can guarantee $100\mathrm{\%}$ accuracy (Jiang et al. 2020). A large fraction of errors in detecting cyclone centres are caused by mistaking a trough for a closed cyclone. Other discrepancies between manually and automatically detected cyclones can be connected with the value taken as the interval for detecting enclosed contours. It can be expected that the number of identified lows increases as the interval decreases (Wernli and Schwierz 2006). Small and/or shallow cyclones are mostly sensitive to the choice of the interval. The denser the counters become, the more relatively small-scale weather systems will appear, while the original systems remain the same, despite the potential extensions of their original boundaries (Lu 2017). It is not easy to find the accuracy of detection methods. While some authors (e.g. Jiang et al. 2020) provide approximately $85\mathrm{\%}$ accuracy compared to manual analysis, these numbers are usually not confirmed by other studies. Qualitative conclusions are given instead. As Neu et al. (2013) states, the largest differences when comparing various methods and schemes can be found in the frequency distributions for short-lived, shallow, and slowly moving cyclones, especially in the Northern Hemisphere than in the Southern Hemisphere. This method’s sensitivity should be borne in mind when using a single method. This is especially true for total cyclone counts and number and role of weak cyclones in the statistics.

When dealing with mid- and high-latitude pressure systems, sometimes the very miscellaneous structure of extratropical cyclones and especially anticyclones need to be processed by algorithms for their automatic detection. Relatively regularly shaped pressure patterns are typical for very deep lows or intense highs, characteristically occurring in the initial development stages. However, complicated structures with more centres can be detected during the mature and dissipation stages. Algorithms for automatic detection have to deal with various sizes (hundreds to thousands of kilometres), depths (from $920$ to $1080$ $hPa$), or translational velocities ($0$ to $100$ $km/h$) of the pressure systems. The general problem to be solved by automatic algorithms is the accurate detection of an outer anti/cyclonic boundary because it is not easy to find a set of closed contours. The detection of short-lived systems (such as thermal lows) can also be quite challenging for some widely used algorithms. Some schemes use a minimum life cycle time in order to omit such systems (see e.g. Hanley and Caballero 2012). Nonetheless, such an approach can cause the omission of important weather features (e.g. heavy precipitation). For instance, this can be brought by short-lived cut-off cyclones. Another problem can arise with open systems, i.e. local pressure maxima/minima without any closed isobars/isohypse. These include lee cyclones or luv anticyclones. It can also happen that algorithms were designed to search for cyclones of specific types (e.g. polar lows, extremely deep cyclones). However, these are not able to detect all other cyclones virtually in various latitudes.

The main reason for automatic detection of pressure patterns is to have an instrument for a great number of (re)analyses: detecting lows and highs manually would be too time-consuming. However, to have this automatic detection available, climatology of various anti/cyclone characteristics can be prepared. This makes it possible to analyse parameters such as spatial distribution, evolution characteristics or frequency of occurrences in various regions (e.g. Dacre and Gray 2009, Lu 2017). This is especially valuable when studying long-term data with respect to the climate change.

In this paper, we focus on the surface pressure-level field.

Currently, neural networks gain a strong interest in various research topics. They can be adopted to detect cyclone and anticyclone centres as well as their tracks (e.g. Matsuoka et al. 2018, Lu et al. 2020). However, the main problem is the dependency on the input data. If the input changes its parameters (projection, resolution, etc.), the entire network must be trained again.

We propose a traditional algorithm based on simple input parameters that does not require any training data or training phase. Input is only a single mean sea-level pressure layer image. The proposed solution does not rely on terrain heights. The problem of this approach would be that we would have to generate a special height map aligned with the data on a pixel basis. If the resolution is low, the height map would accumulate many different heights within a single pixel.

As will be explained and discussed in further detail, the main features of our proposed solution are:

• automatic detection of low- and high-pressure systems from a single image

• low number and predictable behaviour of input parameters

• detect systems only from the mean sea-level pressure layer (MSLP)

• fast and easy to implement solution

In Section 2 we provide a detailed explanation of our solution. Section 3 presents the algorithm results. Our conclusion and final discussion can be found in Section 4.

2. Proposed algorithm

The proposed algorithm detects low- and high-pressure systems directly from the mean sea-level pressure layer (MSLP). This layer 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 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 workflow of the main steps is visualised in Figure A1. In the visualisations, we use a simple greyscale input image for simplicity. However, the algorithm can use any type of an input raster data. The proposed solution consists of two main parts: the detection of candidates (see Section 2.1, first line of the workflow diagram) and the pressure system centre detection within these candidates (see Section 2.2, second line of the workflow diagram).

In the following sections, we describe each part in more detail.

2.1. Candidate 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$ $hPa$ (Rudeva 2008). Therefore, the isobars are usually drawn with steps of $1$, $2$, $4$, $5$, or $10$ $hPa$. We use the same steps for low- and high-pressure systems. With lower values, we can detect more systems. Based on our experiments we recommend using $2$ or $4$ $hPa$ for global models. For local models, capturing a smaller area, we recommend using $1$ $hPa$. From the field of isobars, we select a list of candidates, i.e. the potential centres of pressure systems.

Detecting isobars from raster data of size $W\times 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 polygons and have a fine, pixel-based, resolution. For further computations, these polygons can be simplified (e.g. using polyline simplification from Douglas and Peucker 1973) to improve the performance of the following steps.

The pressure system candidates are selected from the generated, closed isobars. Each candidate must not have any other closed isobars inside.

Secondly, we use filtering based on the area size to remove incorrect candidates. This step also helps to remove unwanted small areas that can originate from terrain influence on the data (mostly in mountain regions). The area can be directly calculated from the polygon. Instead of the precise area calculation, we can use a rough approximation based on the axis-aligned bounding box (AABB).

The area size cannot be easily represented in image pixels because of the projection distortions that can be visualised with Tissot’s indicatrix (Goldberg and Gott 2007). Based on this observation, we use the area size in $k{m}^2$ which makes the algorithm independent of the used projection. However, for this approach, the projection of the input image must be known since the inverse projection formulas are used to convert pixel 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).

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. The smallest tropical cyclones have sizes of the effective radius from $300$ to $400km$ in winter and $200$ to $300km$ in summer (e.g. Rudeva and Gulev 2007, Frame et al. 2017). However, according to Jelenak et al. (2012), extratropical cyclones can be even smaller with a radius less than $50km$. Based on these findings, our recommended value for $T$ is between $10,000k{m}^2$ and $20,000k{m}^2$ to detect tropical and extratropical cyclones. However, the value can be increased to detect only larger and significant systems.

If the area size is smaller than the selected threshold, we check the enclosing isobar (called the ‘parent’). If there is no other isobar inside the parent, the parent becomes a new pressure system candidate. Otherwise, the candidate is rejected.

The filtering process is repeated until there is a change in the candidates.

2.2. Pressure systems detection

The detection of pressure system centres 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.

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 from $dy$).

2.2.1. Low-pressure system

To determine if a pixel represents a low-pressure system, we use its neighbourhood and compare every pixel against the mask with radius $r$ (therefore, the width and height are $2r+1$). A pattern is based on Prewitt operator (see Gonzalez and Woods 2008) with the size of mask itself being larger (full pattern is shown in Tables 1 and 2 a supplementary file). The size of the radius is an input parameter that can be adjusted by the user.

Table 1. Experiment settings. Mask radius $r$ is auto-calculated from area threshold $T$. Minimal centre distance is disabled.

Parameter	Value
Isobars step $C$	2 $hPa$
Area threshold $T$	10ʹ000 $k{m}^2$
Mask radius $r$	automatic
Minimal centre distance $D$	0 $km$

Table 2. Auto-calculated values of mask radius $r$ in pixels for some NWP models with a standard deviation $\sigma$ in [pixels].

Model	Coverage	Resolution	$r$ [pixels]	$\sigma$ [pixels]
ECMWF	Global	$1440\times 721$	9	1.4
ICON	Global	$2880\times 1441$	14	2.1
ICON	Europe	$1094\times 657$	28	4.3
HRRR	USA	$1799\times 1059$	55	7.7

However, based on our experiments, the default size of the radius 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, auto-calculation can lead to large radius and consequently to longer computations. If this is an issue, the radius can be set manually. Usually, based on our experiments, the size of radius selected from range $\left(5,50\right)$ 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 pixels is above $60\mathrm{\%}$, 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\leftarrow (2\ast r+1)\ast (2\ast r+1)$

2: $okX\leftarrow 0$

3: $okY\leftarrow 0$

4: for (yy = y – r; yy ≤ y + r; yy++) do

5: for (xx = x – r; xx ≤ x + r; xx++) do

6: $maskSign\leftarrow mask[yy-(y-r)][xx-(x-r)]$

7: if $maskSign\le signDx[yy][xx]$ then

8: $okX\leftarrow okX+1$

9: end if

10: if $maskSign\le signDy[yy][xx]$ then

11: $okY\leftarrow okY+1$

12: end if

13: end for

14: end for

15: if $okX/maskArea\ge 0.6\mathbf{a}\mathbf{n}\mathbf{d}okY/maskArea\ge 0.6$ then

16: pixel (x, y) belongs to extrema

17: end if

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.

2.2.2. 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. Therefore, we use a mask with swapped signs.

2.2.3. Centre detection

In the previous subsections, we have detected areas of extremas. However, we need to obtain the centre of each system. This process is straightforward. For a low-pressure system, we found the position of minimal value over the system area pixels. On the other hand, if the area corresponds to a high-pressure system, we are looking for a maximal value. The obtained position is used as the centre of the system.

However, in some cases, usually if the input image resolution is small (step is larger than approx. $15km$), the area inside the candidate may not be detected. This usually happens if there is a small area candidate with a large pressure change. If this is the case, we use the centre of the candidate, search its neighbourhood and calculate maximal pressure change. If the change is larger than the isobars step, we mark the candidate centre as a real pressure system centre.

2.3. 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 centres of the systems from the centre position of the candidate contour.

Centres that are too close to each other might cause another problem. If this is a problem, we may apply the optional, user-defined threshold distance $D$ in $km$. In this case, if two systems are closer than the threshold $D$, the system with a smaller area is removed.

2.4. 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 that have been detected, the slower the algorithm.

Algorithm performance can be influenced by detection if a pixel is inside a candidate (see Subsection 2.2.1). In the naive 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.

2.5. Limitations

The proposed solution has its limitations in the input projection. For example, in the case of equirectangular projection, systems near the poles may be detected incorrectly. Hence, instead of one centre, several or no centres may be obtained. This problem is caused by distortion near the poles. Data reprojection might lead to a solution. For low-distorted latitudes (e.g. in interval $\left(-{85}^{\circ }{,85}^{\circ }\right)$), one projection system can be used and a different projection is only applied for polar regions.

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.

3. Experiments and results

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

3.1. Algorithm setup

To evaluate the proposed solution, we have used historical data obtained from the reanalysis of the ECMWF model (dataset ERA5, Hersbach et al. 2021). As an input, several parameters need to be set. Once the parameters are fine-tuned for a given NWP model, they can be used across multiple images. To keep time consistency, there is no need to change parameters for every time input. In our experiments, we have run the proposed algorithm with the settings presented in Table 1.

To show the effect of parameter changes on the overall number of detected systems, we have run experiments with multiple setups. The relation between input parameters, area threshold $T$ and isobars step $C$, is depicted in Figure A2. The results are for a global model (ICON) from a single time, 2021–03-18 12:00 UTC. It can be observed that for a majority of parameter combinations we obtain a similar count of pressure systems in the interval $(75,125)$. On the other hand, for the small values of area threshold and isobars step we obtain too many pressure systems with a lot of false detections.

We recommend to use auto-calculated values of mask radius $r$ from Table 1. For comparison, we present some auto-calculated values in Table 2. The values are obtained from multiple runs over different times. Results are averaged and rounded to whole numbers to represent pixels.

3.2. Evaluation

The evaluation of the proposed solution is conducted on selected datasets from 1998 to 2020. As an input to the algorithm, we have used data from ECMWF reanalysis.

We have used datasets from Europe and the North Atlantic Ocean region provided by Met Office and available at MetOffice (2021), Müller and Floors (2021). There are pressure centres (highs and lows) available as being created and checked by the National Meteorological Institute in the UK; therefore, we consider this data as manually checked. Unfortunately, the results are only available as images. Thus, we have to process the data manually and create a list of the centres of pressure systems. We have selected this dataset as a ground-truth, since it is manually curated.

Another dataset Serreze (2009), we have tested, covers the Northern Hemisphere. It contains only low-pressure systems. Data are obtained via reanalysis using the algorithm from Serreze et al. (1997). However, the spatial resolution is only $2.5deg$ (projected to grid with step of $250km$). Some smaller systems therefore may be lost.

Data covering the North American region can be obtained from Prediction (2022). It contains low- and high-pressure systems and several other parameters that we do not use in our evaluation. However, to our best knowledge, we were not able to find the methodology used to create data. We do not know if they are manually or automatically created or based on some combination of both processes.

To evaluate experiments, we use two metrics: precision and recall. To calculate them, we need to evaluate True Positives ($TP$), False Positives ($FP$) and False Negatives ($FN$).

• $TP$ – systems that are correctly found by the algorithm.

• $FP$ – systems found by the algorithm but not present in the original dataset.

• $FN$ – systems in the original dataset but missing in the results from the proposed solution.

Precision presents percentage of correct systems from obtained set of system. It is calculated as $TP/(TP+FP)$. Recall express percentage of how many systems are correctly predicted from all systems. It is calculated as $TP/(TP+FN)$. In our experiments, we use macro versions of the metrics. They are calculated independently for each time image and averaged at the end.

Since it is considerably time-consuming to process all images from MetOffice (2021) manually, we have created only a random subset. Its size was determined based on Cochran (1977). For a minimal $95\mathrm{\%}$ confidence level, we have to determine a sample of at least $385$ images. Based on this observation, to determine the precision and recall of the proposed solution we have randomly selected $385$ pressure lows and highs (with the same ratio) from 1998 to 2020. We have used the same number of images for other datasets as well.

Our detected pressure system is marked same ($TP$) as in the dataset if both systems are of the same type and their centres are within a given radius to each other (the limit is set higher because of the differences in the source and the origin of the datasets). If our solution detected a system that is not in the vicinity of a system from the dataset, we marked this system as $FP$. Otherwise, the missing systems are marked as $FN$. In the experiments, we have calculated the values for low- and high-pressure systems.

The results are shown in Table 3. For Serreze (2009) there are no results for a high-pressure system, because it was not presented in the input dataset. Average value is also omitted.

Table 3. Comparison of the proposed solution for different datasets: Met Office analysis (MO), NWS Weather Prediction (2022) (NWS) and Serreze (2009) reanalysis (S).

		Precision (%)			Recall (%)
Radius (km)	Type	MO	NWS	S	MO	NWS	S
500	low	84	52	56	75	44	72
500	high	88	63	–	78	52	-
500	average	86	58	–	77	48	-
1000	low	88	65	63	80	61	82
1000	high	91	77	–	84	74	-
1000	average	89	73	–	81	67	-

The detection of highs (anticyclone) is slightly more similar to the systems detected in datasets with which we compared the proposed solution. Pressure highs occupy larger areas and are generally easier to detect.

Several main reasons for differences between the systems detected from our solution and provided datasets can be recognised.

First, the proposed solution is universal. Hence, it can be used for the entire world, whereas Met Office has traditionally focused on Europe and North Atlantic. This, for example, enables identifying cyclones at a very early stage.

We have also compared the quality of our solution with two other datasets – broader coverage with Serreze (2009) and local NWS Weather Prediction (2022). However, these datasets are mainly obtained as outputs from algorithms. Apart from this, the source of Serreze (2009) has a lower resolution and therefore some smaller systems are not present at all. Our experiments shows that our proposed solution is mainly capable to detect the same important systems. However, the overall quality is lower than in the case of ground-truth data from Met Office.

The second reason is time continuity. The proposed solution aims at detecting low- and high-pressure systems from a single input. Met Office takes the previous and the following time steps into consideration (time consecution is conserved). For example, filling up low-pressure systems can be detected manually for a longer period of time by meteorologists from Met Office.

The third reason is the data itself. We use reanalysis from the ECMWF that can be slightly different from the provided datasets.

3.2.1. Comparison with other methods

Comparison between the proposed solution and Met Office shows similarity $92\mathrm{\%}$ as we stated above. However, it is also important to demonstrate how accurate are other automated algorithms.

One of the latest algorithms for the detection of cyclones was published by Jiang et al. (2020). They use the mean sea-level pressure (MSLP) or geopotential height to identify extra-tropical cyclones. The accuracy is measured on a sample of $500$ cyclones, which are determined by the algorithm and compared with manually determined data. The average match is around $85\mathrm{\%}$. Thus, the results are similar to our proposed algorithm.

Further, Jiang et al. (2020) mentions that their results outperform three other solutions called M01 (König et al. 1993, Geng and Sugi 2001, Rudeva and Gulev 2007), M02 (Murray and Simmonds 1991, Simmonds and Keay 2000, Lim and Simmonds 2007) and M04 (Flaounas et al. 2014) which were published earlier. Other published algorithms do not provide their accuracy against manually created data sets. They usually demonstrate the accuracy of individual cases or measure the overall coherence of selected data.

3.2.2. Visual comparison

A sample comparison with data from Met Office is in Figure A3. Triangles (red) indicate high-pressure systems, whereas dots (blue) are used for low-pressure systems. The left part shows the results of the proposed solution in the equirectangular projection (the projection was adopted directly from the source of the dataset). The right part contains the Met Office data. They can be considered as manually checked by a meteorologist. The projection for this image is different, but the borders of land and the range of latitude and longitude area are the same for both images.

In Figure A3, our algorithm detected fewer systems: the missing systems are observed in the Atlantic Ocean ($\mathrm{\#}1$), near Scotland ($\mathrm{\#}2$), and over Northern Africa ($\mathrm{\#}3$). Systems $\mathrm{\#}1$ and $\mathrm{\#}2$ are near the end of their lifetime. Thus, they are not recognised completely by the proposed algorithm. The difference in System $\mathrm{\#}3$ can be caused by the difference in the input dataset.

Some additional results may be found in a supplementary file.

We have tested our solution for hurricane tracking. We have run the algorithm for consecutive images obtained from the ICON model. Each image was processed independently with no knowledge about the other images. To track the hurricane, a user must select one detected point belonging to the hurricane (this point does not have to be the starting point). The next points are found in consecutive images by simply finding the nearest point based on the previous or next image in the sequence.

The result for two hurricanes near Japan in September 2021 can be seen in Figure A4. There is no post-processing applied to the detected tracks.

3.3. Performance

The performance of the algorithm is also important. If the algorithm runs fast, it can be used to reanalyse a large amount of data in a short time. In addition, it also makes it possible to modify the input parameters and test quickly various scenarios.

We have performed several tests with different resolutions and input coverage (global and local data). We have used an older CPU with 4.0 GHz (Core i7-4790 K) and 32GB RAM.

The averaged run-times with the standard deviation are available in Table 4. We can see that the performance of the algorithm is affected mainly by the number of detected pressure systems. The regional models contain fewer systems. In addition, their performance is better even though the resolution is similar to that of global models. The resolution also affects the performance. With the higher resolution, the run time is longer even though there is a similar number of detected systems.

Table 4. Measured average run-times in [ms] with a standard deviation $\sigma$ in [ms]. The column Systems denotes the average number of detected systems.

Resolution	Coverage	Systems	Time [ms]	$\sigma$ [ms]
$1440\times 721$	Global	84	1483	148
$2880\times 1441$	Global	95	1960	163
$1094\times 657$	Europe	13	249	9
$1799\times 1059$	USA	21	503	18

4. Conclusions

An algorithm for automatic detection of low- and high-pressure systems has been presented in the paper.

The proposed solution uses only the MSLP layer for detection. This approach is similar to many existing solutions (Neu et al. 2013). Sometimes vorticity is used (e.g. Simmonds et al. 2008), which could improve identifying cyclones at a very early stage of development. However, vorticity is not suitable for determining high-pressure system centres. Therefore, it cannot be deployed in our algorithm. Moreover, the presented solution does not need direct terrain filtering used by some other algorithms. Our algorithm automatically discards many small systems that are usually scattered over higher mountains (e.g. the Himalayas, the Rocky Mountains, and the Andes).

Our solution provides advantages in several areas. In general, it is easy to implement because it is based only on standard image processing algorithms. In addition, the performance of the proposed algorithm is sufficient for general use. Since we do not use neural networks, there is no training phase and no need for training data. Last but not least, not only can it detect low-pressure systems, as is usually standard by other algorithms, but also it can localise the centres of high-pressure systems.

However, it is challenging to validate the quality of results. The correct position of the pressure system cannot be determined exactly (see Neu et al. 2013). There is no universal definition that would precisely determine the centre of a low- or high-pressure system. If systems are created manually, each author may prefer a slightly different position and properties of the system. This results in smaller or larger differences, which may partly be influenced by different data interpolations of the analysed field. Automated solutions, based on algorithms, share similar problems. Each solution is adjusted to detect systems that the authors of the algorithm consider important.

Some studies (e.g. Neu et al. 2013) compare several algorithms. However, they conclude that there are no ‘correct’ and ‘incorrect’ results. All algorithms can be considered ‘correct’, depending on the requirements that the system must meet. Based on these findings, our solution was not directly compared with the existing algorithms. Instead, as a main comparison, the selection of $385$ manually created image-analyses by Met Office, which are widely accepted by the meteorologist, was compared. The comparison shows $83\mathrm{\%}$ through $92\mathrm{\%}$ agreement between the two approaches. This is quite similar success to other methods (e.g. Dacre and Gray 2009, Hanley and Caballero 2012, Neu et al. 2013, Jiang et al. 2020). The most significant differences can be found in the detection of cyclones at the beginning or end of the lifespan stage. Met Office meteorologists detect cyclones in these stages more often than does the presented solution.

Overall, the problem with experiments is the absence of large, manually created datasets that can be used for extensive evaluation. If we compare results with another algorithm, we are still missing the comparison with a real situation.

The proposed solution requires to set up two parameters, while the others can be auto-calculated. However, the need for input parameters can be found in other solutions as well (Raible et al. 2008). Based on our experiments, the algorithm is only a little sensitive to the values of input parameters. Most of the time, the number of detected systems is in a stable interval regardless of the values of the parameters. The problems with the detection of too many systems arise if the values are set too small. If this is an issue, the number of detected systems can be adjusted by the change of input parameter values. Due to the fast performance, the result of parameters change can be quickly evaluated.

We provide a simple demo application of the proposed solution. Its source code is available at https://github.com/MartinPerry/cyclone-detector.

Currently, the algorithm is designed to detect systems only from a single image without the knowledge of the previous and next time. This can sometimes lead to a missing system detection, particularly during the initial and decaying stages. In our future work, we would like to include time consistency that could further improve results. As presented in Figure A4, we have tested single image time consistency on hurricane tracking. The proposed algorithm shows consistent behaviour and the majority of centres is detected.

In the perspective of our further research, we will use the presented solution for a larger reanalysis of historical data in order to track changes of cyclones and anticyclones over the last decades. This will help us to analyse how climate change affects the trajectories and frequencies of pressure systems.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/14498596.2022.2134221

Additional information

Funding

This work was supported by the Department of Atmospheric Physics, Faculty of Mathematics and Physics, Charles University, and partially supported by the Czech Ministry of Education, Youth and Sports, project [PUNTIS (LO1506)].

Appendix A.

Code and data availability

The code of the algorithms in C++ can be downloaded from https://github.com/MartinPerry/cyclone-detector . The data that support the findings of this study are ERA5 hourly data at single levels from 1979 to the present. They are openly available from Copernicus at http://doi.org10.24381/cds.adbb2d47 (Hersbach et al. 2021). Data for comparisons are obtained from Met Office and archived by WetterZentrale. They are publicly available at https://www.wetterzentrale.de/de/reanalysis.php?map=1model = bra var = 45 (Müller and Floors 2021)

Figure A1. Workflow diagram of the algorithm’s most important steps.

Figure A2. Relation between input parameters, area threshold and isobars step, with respect to the number of detected pressure systems.

Figure A3. Comparison of detected pressure systems from the proposed solution (left) and Met Office (right). The triangles (red) indicate high-pressure systems, and the dots (blue) represent low-pressure systems. The numbers denote systems missing in our solution. Data from 2021–01-27 00:00 UTC.

Figure A4. Tracking of two hurricanes near Japan. Left one is Chanthu (2021–09-07 to 2021–09-18), right one is Mindulle (2021–09-23 to 2021–10-01). Ground truth data obtained from U.S. Naval Research Laboratory and plotted with (WikiProject Tropical cyclones 2021).