Using Persistent Homology to Quantify a Diurnal Cycle in Hurricane Felix
Sarah Tymochko, Elizabeth Munch, Jason Dunion, Kristen Corbosiero, Ryan Torn
UUsing persistent homology to quantify a diurnal cycle in hurricanes
Sarah Tymochko , Elizabeth Munch , Jason Dunion , Kristen Corbosiero , and Ryan Torn Michigan State University, Dept. of Computational Mathematics, Science and Engineering Michigan State University, Dept. of Mathematics Cooperative Institute for Marine and Atmospheric Studies, University of Miami Hurricane Research Division, NOAA/Atlantic Oceanographic and Meteorological Laboratory University at Albany - SUNY Albany, Dept. of Atmospheric and Environmental Sciences
Abstract
The diurnal cycle of tropical cyclones (TCs) is a daily cycle in clouds that appears in satellite images and mayhave implications for TC structure and intensity. The diurnal pattern can be seen in infrared (IR) satellite imageryas cyclical pulses in the cloud field that propagate radially outward from the center of nearly all Atlantic-basinTCs. These diurnal pulses, a distinguishing characteristic of this diurnal cycle, begin forming in the storm’s innercore near sunset each day, appearing as a region of cooling cloud-top temperatures. The area of cooling takeson a ring-like appearance as cloud-top warming occurs on its inside edge and the cooling moves away from thestorm overnight, reaching several hundred kilometers from the circulation center by the following afternoon. Thestate-of-the-art TC diurnal cycle measurement in IR satellite imagery has a limited ability to analyze the behaviorbeyond qualitative observations. We present a method for quantifying the TC diurnal cycle using one-dimensionalpersistent homology, a tool from Topological Data Analysis, by tracking maximum persistence and quantifying thecycle using the discrete Fourier transform. Using Geostationary Operational Environmental Satellite IR imageryfrom Hurricanes Felix and Ivan, our method is able to detect an approximate daily cycle.
Keywords:
Topological Data Analysis, Atmospheric Science, Diurnal Cycle, Image Processing
The field of atmospheric science has numerous observa-tion platforms that provide high space and time resolutiondata, but has yet to find methods which can quantify theintuitive patterns explicitly. Meanwhile, the young fieldof Topological Data Analysis (TDA) encompasses meth-ods for quantifying exactly these sorts of structural intu-itions seen by atmospheric scientists. This paper mergesthese two fields by using persistent homology, a now well-established tool in TDA, to quantify a diurnal cycle ob-served in a hurricane using Geostationary OperationalEnvironmental Satellite (GOES) infrared (IR) satellitedata.Persistent homology, and more generally TDA meth-ods, has found significant success in rather disparate ap-plications by finding structure in data and using this in-sight to answer questions from the domain of interest.For instance, Giusti et al. used the homology of randomsimplicial complexes to investigate the geometric organi-zation of neurons in rat brains ( [9]). Persistent homologyhas also been used to understand periodicity in time seriesarising from biological ( [3, 22]) and engineering applica-tions ( [12]), as well as for image analysis ( [2, 23, 24, 27]).This paper presents an application of the use of bothtime series and image analysis using TDA to study a daily cycle in hurricanes that is of interest in the fieldof atmospheric science. The diurnal cycle of tropicalcyclones (TCs) has been described in previous studies( [4–7, 14, 16, 20, 21, 29, 30]) that provide evidence of theregularity of this cycle as well as its potential impacts.This diurnal pattern can be seen in GOES IR imagery ascyclical pulses in the cloud field that propagate radiallyoutward from TCs at speeds of 5-10 m s − ( [4, 6, 7]).These diurnal pulses, a distinguishing characteristic ofthe TC diurnal cycle, begin forming in the TC’s corenear the time of sunset each day and appear as a regionof cooling cloud-top temperatures. The area of coolingthen takes on a ring-like appearance as marked cloud-topwarming occurs on its inside edge and it moves away fromthe storm overnight, reaching several hundred kilometersfrom the TC center by the following afternoon. Observa-tions and numerical model simulations indicate that TCdiurnal pulses propagate through a deep layer of the TCenvironment, suggesting that they may have implicationsfor TC structure and intensity ( [4–7, 16]).The current state of the art TC diurnal cycle mea-surement has a limited ability to analyze the behaviorbeyond qualitative observations. This paper presents amore advanced mathematical method for quantifying theTC diurnal cycle using tools from TDA, namely one-dimensional persistent homology to analyze the holes in1 a r X i v : . [ c s . C V ] J a n space. This research aims to detect the presence of thediurnal cycle in GOES IR satellite imagery and to trackthe changes through a time series.In this paper, we present a method of automaticallydetecting and quantifying periodic circular structure insatellite imagery as well as the results of our method ap-plied to GOES IR hurricane imagery for Hurricane Felixin 2007 and Hurricane Ivan in 2004. The naive combina-tion of persistent homology with the IR imagery did notshow a recurring pattern due to drastically variable valuesin IR brightness temperature data. Despite this, lookingat the data, there is a clear circular feature that is visiblein the imagery. Thus we develop a more sophisticatedmethod using tools including the distance transform andone-dimensional persistent homology to detect the TCdiurnal cycle quantitatively using maximum persistence.Using our method, we are able to detect this 24-hour cy-cle in both hurricanes automatically, improving upon theexisting qualitative methods. Previous research has documented a clear diurnal cy-cle of cloudiness and rainfall in TCs: enhanced convec-tion (i.e., thunderstorms) occurs overnight, precipitationpeaks near sunrise, and upper-level cloudiness (i.e., thecirrus canopy) expands radially outward throughout theday, reaching its maximum areal coverage in the earlyevening hours ( [4–7,14,16,20,21,29,30]). To quantify theexpansion and contraction of the cirrus canopy, Dunionet al. used GOES satellite IR imagery to examine thesix-hour cloud-top temperature differences of major hurri-canes in the Atlantic basin from 2001 to 2010 ( [7]). Theyfound that an area of colder cloud tops propagated out-ward around 5-10 m s − over the course of the day, withwarming temperatures on its inner edge. More recently,in ( [4]), Ditchek et al. expanded Dunion et al.’s workto include all tropical cyclones in the Atlantic basin from1982 to 2017 and found that the diurnal pulse is nearlyubiquitous, with 88% of TC days featuring an outwardlypropagating pulse.Despite the consistent signature and documentationof this diurnal cloud signature, open questions remain asto how the diurnal cycle is linked to inner-core convectiveprocesses and whether it is a column-deep phenomenonor mainly tied to upper-level TC cloud dynamics relatedto incoming solar radiation ( [4, 20, 21, 30]). Investigat-ing these questions is relevant to TC forecasting as thediurnal cycle of clouds and rainfall has implications forforecasting storm structure and intensity, as evidencedby the diurnal cycle in objective measures of TC intensityand the extent of the 50-kt wind radius documented byDunion et al. Additionally, and especially relevant to thecurrent work, most of the papers above have identified thepulse using subjective measures of cloud-top temperaturechange and timing ( [4, 7]). The current work seeks to Figure 1: An example matrix, M (top left) and corre-sponding persistence diagram (top right). Second row:The black portions are sublevel sets, M r , where r =0 . , . , and 0 .
7. The existence of a point far fromthe diagonal in the persistence diagram shows that thereis a prominent circular structure; while the other pointsare caused by the smaller circles.quantify the pulse to determine its true periodicity usingpersistent homology, a topological tool that is particu-larly effective at capturing the type of patterns visible inthe pulse.
Persistent homology is a tool from the field of TDA whichmeasures structure in data. This data can start in manyforms, including as point clouds or, as in the case of thiswork, as a function on a domain. In this section, we willbriefly review the necessary background to understandcubical complexes and persistent homology, and refer theinterested reader to [8, 10, 11, 18] for a more complete in-troduction.
In this section, we largely follow Chapter 2 of [11] withthe caveat for the informed reader that because we usehomology with Z coefficients, we can be lazy about ori-entations of cubes. In addition, our data consists of 2Dimages, so we need only define cubes up to dimension 2.An elementary interval is a closed interval I ⊂ R of theform [ (cid:96), (cid:96) + 1] or [ (cid:96) ] for (cid:96) ∈ Z , which are called nondegen-erate and degenerate respectively. An elementary cube Q ∈ R is a product of elementary intervals Q = I × I .The dimension of Q , dim( Q ), is the number of nondegen-erate components of Q . Note that 0-dimensional cubesare just vertices at the points on the lattice Z × Z in R , 1-dimensional cubes are edges connecting these ver-tices, and 2-dimensional cubes are squares. Let K de-note the set of all elementary cubes in R and K d ⊂ K the set of d -dimensional cubes. A set X ⊂ R is cubi-cal if it can be written as a finite union of elementary2ubes. Then we denote the associated cubical complex as K ( X ) = { Q ∈ K | Q ⊂ X } , with the d -dimensional subsetdenoted K d ( X ) = { Q ∈ K ( X ) | dim( Q ) = d } . If Q ⊆ P ,then we say Q is a face of P , denoted Q ≤ P . If Q (cid:40) P ,then Q is a proper face of P , denoted Q < P , and isadditionally a primary face of P if dim( Q ) = dim( P ) − m × n matrix,can be viewed as a function M : D → R where D = { ( i, j ) | ≤ i < m, ≤ j < n } . We will model thisas a function defined on a particularly simple cubical set K = K ([0 , m ] × [0 , n ]). For simplicity, we denote by s i,j the square [ i, i + 1] × [ j, j + 1]. So, given a matrix M , weequivalently think of this data as a function M : K → R where we set M ( s i,j ) equal to the matrix entry M i,j andset M ( P ) = min s i,j >P M ( s i,j ) for all lower dimensionalcubes P . Note that we will abuse notation and use M todenote both the original matrix and representation as afunction with domain K . Homology ( [10]) is a standard tool in algebraic topol-ogy which provides a vector space H k ( X ) for each di-mension k = 0 , , , . . . for a given topological space X .The different dimensions measure different properties ofthe space. In particular, for this work we are interestedin 1-dimensional homology; i.e. when k = 1. The 1-dimensional homology group measures the number of loopsin the space; equivalently, we can think of this as the num-ber of holes in the space. In particular, if we look at theblack region in each of the examples in Fig. 1, the rankof the first homology for each is (1,2,1).The exact definition of homology is as follows. For anycubical set L (which for the purposes of this discussionwill always be a subset of K ), we have sets giving thecubes of different dimensions: K i ( L ) for i = 0 , ,
2. An i -chain is a formal linear combination of i -simplices in L , c = (cid:80) Q j ∈K i ( L ) a j Q j , with coefficients a j ∈ Z . Wecan of course add these objects by setting ( (cid:80) a j Q j ) +( (cid:80) b j Q j ) = (cid:80) ( a j + b j ) Q j and multiply by a constant.Thus, the collection of all i -chains forms a vector space C i ( L ).We define a linear transformation δ i : C i ( L ) → C i − ( L )called the boundary map, by setting δ i ( Q ) = (cid:80) P wherethe sum is over the primary faces P < Q . The ker-nel of δ , Ker( δ i ), (that is, the set of elements of C ( L )which map to 0) is generated by closed loops in L . Theimage of δ , Im( δ ), is generated by boundaries of 2-cells. Then the 1-dimensional homology group is definedto be H ( K ) = Ker( δ ) / Im( δ ). An element of this group γ ∈ H ( L ), represents an equivalence class of loops whichcan differ by collections of 2-cells. Normally a group, however, we are working with field coeffi-cients. Again, notice that because we are working with Z coefficients,the book-keeping normally needed for orientation is unnecessary. For a static space L , H ( L ) measures information aboutthe number of loops. Persistent homology takes as inputa changing topological space, and summarizes the infor-mation about how the homology changes.Let an m × n R -valued matrix M be given. Fix a func-tion value r ∈ R and let M r = f − ( −∞ , r ]. That is, M r isthe subset of squares in K m × n which have value at most r in the matrix, along with all edges and vertices which arefaces of any included square. M r is often called a sublevelset of M . See Fig. 1 for M r regions corresponding to anexample image.These spaces have the property that M r ⊆ M s for r ≤ s , thus we can consider the sequence M r ⊆ M r ⊆· · · ⊆ M r k for any set of numbers r < r < · · · r k . Thissequence of spaces is called a filtration. For each of thespaces, we can compute the homology group H p ( M r i ).The inclusion maps give rise to linear maps H p ( M r ) → H p ( M r ) → · · · → H p ( M r k ) . It is these maps that westudy to understand how the space changes. In partic-ular, when we are focused on 1-dimensional homology( k = 1) as in this study, a loop is represented by anelement γ ∈ H ( M r i ). We say that this loop is bornat r i if it is not in the image from the previous space;that is, γ (cid:54)∈ Im( H ( M i − ) → H ( M i )). This same loopdies at r j if it merges with this image in M r j ; that is, γ ∈ Im( H ( M i − ) → H ( M j )) where we abuse notationby using γ to both refer to the class in H ( M r i ) andthe image of this class under the sequence of maps in H ( M r j ). We refer to r j − r i as the lifetime of the class.A persistence diagram, as seen in the right of Fig. 1,is a collection of points where for each class which is bornat r i and dies at r j is represented by a point at ( r i , r j ).The intuition is that a class which has a long lifetime isfar from the diagonal while a class with a short lifetimeis close. In many cases, a long lifetime loop implies thatthere is some sort of inherent topological feature beingfound, and thus that this point far from the diagonal isimportant, while short lifetime loops are likely caused bytopological noise due to sampling or other errors in thesystem. In the example of Fig. 1, there is a prominentoff-diagonal point which shows that the function definedby the matrix has a circular feature. Thus, a commonmeasure for looking at the persistence diagram when in-vestigating a single, circular structure is the maximumpersistence, defined asMaxPers( D ) = max ( r i ,r j ) ∈ D r j − r i (1)for a given persistence diagram, D . The data was given in the form of storm-centered GOESIR (10.7- µ m) satellite imagery. The 10.7- µ m long chan-3igure 2: Original satellite imagery from Felix GOES-12data set (left) and Felix GridSat-GOES data set (right)at approximately the same time.nel detects IR energy emitted from the Earth and is notstrongly affected by atmospheric water vapor. Thus, thisparticular channel is useful for detecting clouds at alltimes of the day and night and is ideal for tracing thediurnal evolution of the TC cloud fields. For our analy-sis, we used two different types of data sets. These twodata sets have the same native spatial resolution, but dif-fer in temporal resolution. The first type (hereafter theGOES-12 data sets), utilizes brightness temperatures de-rived directly from GOES-12 4-km IR satellite imageryand consists of data in hourly increments, with the excep-tion of 0415 and 0515 UTC each day in Hurricane Felix,and 0445 and 0545 UTC each day in Hurricane Ivan (dueto the GOES-12 satellite eclipse period). Imagery wasremapped such that each pixel has a spatial resolution of2 km and each image covers a total area of approximately1500 km × ×
752 matrix.This remapping was performed using the Man ComputerInteractive Data Access System (McIDAS; [15]) in orderto generate storm-centered satellite images that were fo-cused on the relevant TC environment. The McIDAS 4km to 2 km remapping procedure replicates the original8-bit grayscale values such that none of the original pixelinformation is lost. The second type is the GridSat-GOES( [13]) data set and consists of data in 3-hour incrementswith the exception of 0600 UTC each day. Each pixel hasa resolution of approximately 8 km and each image coversa total area of approximately 2400 km × ×
301 matrix. This data is cropped to a191 ×
191 matrix to approximately match the area cov-ered by the first set of data. The cropped version coversa total area of approximately 1530 km × · . . − . / . . Some images in the GridSat-GOES data set also contain missing values where thebrightness temperature for certain pixels was not recordedand is instead assigned a fill value. In order to preventthese values from impacting our results, we interpolatevalues for these pixels. For a given pixel with a missingvalue, we compute the average value of a 5 × M ( t ),from the Felix GOES-12 data set; Top right: thresholdedsubset, M ( t ) µ where µ = 80; Bottom left: distance trans-form function; Bottom right: corresponding persistencediagram.tered at the pixel, not including the pixels in this rangethat also have missing values.For Hurricane Felix, we studied both types of datasets to test the flexibility of the method across spatialand temporal resolution. The Felix GOES-12 data setspans 2 to 4 September 2007, while the Felix GridSat-GOES dataset spans spanning 31 August to 6 September2007. For Hurricane Ivan, we used the GOES-12 data set,which spans 30 August to 1 September 2004. We present a method of detecting and quantifying peri-odic circular structure representing the TC diurnal cyclein IR satellite imagery. Our method combines existingmethods from the fields of image processing, topologicaldata analysis, and signal processing.Initially, we have a time series of IR satellite images,represented as a matrix of pixel values S ( t ) for time t .The TC diurnal pulse is propagating outward throughthe day; thus, in order to see the movement and changesin the GOES satellite brightness temperature, we con-sider the difference in matrices six hours apart ( [4, 7]).For all times t , given the original brightness tempera-ture image S ( t ), we compute the six-hour differences, M ( t ) = S ( t + 6) − S ( t ). While circular features are vi-sually prominent in the data, simply using persistence onthe difference data did not show any relevant features.This discrepancy is due to the extreme differences in thefunction values between the circular sections which pre-4ents the sublevel sets from containing the full circularstructure until very late in the filtration. Thus, we areunable to detect the circular structure from M ( t ) and wedefine a new function on the difference matrix using thefollowing method. Fix a threshold µ and let M ( t ) µ bethe subset of M ( t ) which has function value less than µ .This method results in a binary matrix defined entry-wise,with M ( t ) µ [ i, j ] = 1 if M ( t )[ i, j ] < µ , and M ( t ) µ [ i, j ] = 0otherwise. We will address this choice of threshold inSec. 5.2; however, we will focus on the case µ = 80 de-grees for most of our analysis. Note that because M ( t ) isa difference of two images, the threshold is not isolatingall pixels above a certain temperature, but rather thosepixels that increase in value by at least 80 degrees overthe six hours.After thresholding, we now have binary images; how-ever, the persistent homology is uninteresting, as thereare only two possible sublevel sets. In order to createa greyscale image that maintains the visually apparenttopological structure of the image, we apply the distancetransform, a method from the field of image processing( [25, 26]), which gives a new matrix D ( t ). The dis-tance transform is calculated as follows: given any pixel s i,j in an image M ( t ) represented as a matrix of pix-els, D ( t ) i,j = min d ( s i,j , x ) where x is a 0-valued pixeland d is any distance metric. In this application, wespecifically use the L ∞ distance, also known as the chess-board distance; however, the distance transform is definedfor any distance metric, the Euclidean distance trans-form being the most commonly used. Given two pixels, s i ,j , s i ,j the L ∞ distance between them is calculatedas d ( s i ,j , s i ,j ) = max {| i − i | , | j − j |} . This definesa distance on the pixels, which are the 2-cells in the cubi-cal complex. The distance can be extended to the lowerdimensional cells in the same manor as described in Sec.3.1.To calculate the distance transform, we use the pythonsubmodule scipy.ndimage , specifically the function distance transform cdt with the chessboard metric.Therefore, each entry in D ( t ) corresponds to the minimaldistance to an entry where M ( t )[ i, j ] ≥ µ . The distancetransform D ( t ) is then scaled by the resolution for eachdata set in order to convert the distance units to kilome-ters instead of pixels. For the GridSat-GOES data set,we scale by a factor of 8 km/pixel while for the GOES-12data sets, we scale by a factor of 2 km/pixel. We thencompute sublevel set persistence on the function D ( t ) us-ing the cubtop method in Perseus ( [17, 19]), which cal-culates persistent homology for cubical complexes usingconcepts from discrete Morse theory. Note that Perseusrequires an integer value filtration function on the cubi-cal complex; thus, we chose to use the chessboard metricfor the distance transform. Figure 3 shows an example ofeach step described so far.For each six-hour difference in each data set, we ap-ply the steps described above, then calculate maximumpersistence as defined in Eqn. 1. By plotting maximum persistence over time, we can see how the most promi-nent circular feature changes through the progression ofthe day and life of the TC. This plot should show an oscil-latory pattern, detecting the change in the diurnal cyclethroughout the day. In order to quantify this oscillatorypattern, we use the Fourier transform. In general, theFourier transform is a commonly used method for inves-tigating periodicity of time series ( [1]) by decomposinga wave into a sum of sinusoids with different frequencies.Since we are working with discrete data, we will workwith the discrete Fourier transform (DFT). Let T be thetime between discrete samples, then let t k = kT where k = 1 , . . . , N −
1. Then, the discrete Fourier transformis F n = (cid:80) N − k =0 f ( t k ) e − πink/N . This converts a functionfrom the time domain to the frequency domain. Thepower spectrum of F n can be estimated by calculatingthe square of the absolute value of the discrete Fouriertransform, | F n | .Using the DFT, we calculate the most prominent fre-quency in the data in order to determine how often thecyclic behavior repeats. Note, to use the discrete Fouriertransform the time steps must be equal; however, becauseof the missing times in our data, this is not the case.Therefore, we approximate the maximum persistence atthese values by adding a point along the line between thetimes immediately before and after the missing time. Ad-ditionally, we must truncate the maximum persistence toonly include the days where we have the data for the en-tire day. This means truncating the Felix GridSat-GOESdata set to include only 1-4 September 2007, the FelixGOES-12 data set to include only 1-2 September 2007,and the Ivan GOES-12 data set to include only 30-31 Au-gust 2004. The discrete Fourier transform was calculatedusing the python submodule numpy.fft .We first calculate the Fourier transform using the func-tion fft , then calculate the frequency bins using fftfreq .Using this information, we plot the approximate powerspectrum for each data set. Note, if working with realdata as we are in this application, the power spectrumwill be symmetric for positive and negative frequencies;therefore, we only need to look at the positive frequen-cies. Picking the frequencies corresponding to the highestpeaks in each power spectrum gives us the frequency ofthe most prominent periodic signal in the data. Usingthis frequency, we can calculate the period of the oscil-latory pattern and quantify the signal we are detecting.To verify the detected signal matches the visual oscilla-tory pattern we see, we can reconstruct the sinusoid cor-responding to the most prominent periodic signal usingthe inverse discrete Fourier transform, using ifft . Thenplotting this reconstructed sinusoid over the original data,we can see how closely this signal matches the patternsin the data.5igure 4: Maximum persistence plotted over time for alldata sets using threshold µ = 80 in addition to the recon-structed versions, created using inverse Fourier transform.Gray vertical lines separate days according to UTC. After the data is prepared, we apply the steps describedin Sec. 4.2 to each data set. For the two Hurricane Fe-lix data sets, as they are from the same hurricane, wewould expect the results to be similar despite the tem-poral and spatial resolution differences. Plotting the cal-culated maximum persistence over time, we get the timeseries plotted as solid lines in Fig. 4. The plots show anoscillatory pattern for all three data sets which appearsto repeat approximately daily.To verify the periodicity of the oscillatory pattern, weapply the discrete Fourier transform and calculate thepower spectrum for each data set. Each power spectrumis shown in Fig. 5. Picking the frequencies correspond-ing to the highest peaks in each power spectrum givesus the frequency corresponding to the most prominentsinusoidal signal in the data. The maximum peaks inFig. 5 give a frequency of 0.976 cycles per day for theFelix GridSat-GOES data set, 0.979 cycles per day forthe Felix GOES-12 data set, and 1.0 cycles per day forthe Ivan GOES-12 data set. We use this frequency, f , tocalculate the periodicity of the cycle by calculating using1 /f , giving the period of the sinusoid in days per cycle,then multiply by 24 to rescale to hours per cycle. Doingso gives the result that the cycle is repeating every 24.6hours for the Felix GridSat-GOES data set, every 24.5hours for the Felix GOES-12 data set, and 24.0 hours forthe Ivan GOES-12 data set. Figure 5: Power spectrum for each data set. The high-est peak on the Felix GridSat-GOES power spectrum oc-curs at a frequency of approximately 0.976 cycles/day,the peak on the Felix GOES-12 power spectrum occursat a frequency of approximately 0.979 cycles/day, and thepeak on the the Ivan GOES-12 power spectrum occurs ata frequency of 1.0 cycles/day.Using the most prominent frequency for each data set,we calculate the inverse Fourier transform, and plot thesereconstructed sinusoids over the original data. These si-nusoids, plotted as the lighter dashed lines in Fig. 4,closely resemble the patterns exhibited by the originalmaximum persistence versus time plots; therefore, theseapproximately 24 hour patterns visible in the plots arealso detected mathematically, which verifies the claimthat our method is detecting a daily cycle in each dataset. Additionally, since for both Hurricane Felix data sets,the plots of maximum persistence against time seem tomatch and both have similar detected periodicity fromthe DFT, our method seems robust to the temporal andspatial resolution differences in these two data sets. The method described involves a choice of threshold, sowe used a variety of thresholds, µ ∈ { , , . . . , } ,to test the sensitivity of our method to the parameterchoice. For both data sets of Hurricane Felix, our methodis very robust to the choice of threshold. In Fig. 6, the toprow are plots that represent maximum persistence versustime for the Hurricane Felix data sets using a variety ofthresholds. There is a clear periodic pattern for bothdata sets across all the thresholds shown. In fact, forall thresholds tested µ ∈ { , . . . , } , the period isconsistent at 24.6 hours for the Felix GridSat-GOES dataset and 24.5 hours for the Felix GOES-12 data set. For µ <
35 and µ >
90 the Fourier transform is unable to pickup the daily pattern in the Felix GridSat-GOES data set.For Hurricane Ivan, the plot is shown on the bottomrow of Fig. 6 for thresholds µ ∈ { , , . . . , } . For6igure 6: Maximum persistence vs time plot for Hurri-cane Felix (top row) and Hurricane Ivan (bottom row).Hurricane Felix results are shown for all thresholds µ ∈{ , , . . . , } while Hurricane Ivan results are shownfor µ ∈ { , , . . . , } .Figure 7: A example comparison of distance transformimages with and without opening for the Felix GOES-12data set. At left, the noise in the center of the hurri-cane causes the distance transform to fill in. At right,performing opening gets rid of the small noisy point, andthe distance transform does not get filled in.all of the threshold value shown, our method consistentlydetects a 24.0 hour period. This is a smaller range ofthreshold values than those that detect a daily cycle inHurricane Felix, but for thresholds µ ∈ { , , } , ourmethod detects a daily cycle in all three data sets. Thus,the method may require some parameter tuning, but ouranalysis of these three data sets gives a range of values tostart with when testing new data sets. While the above method detects a daily cycle, there aresome instances where the six-hour differencing introducesnoise because of varying behavior in the center of the hur-ricane. The left image of in Fig. 7 shows an example ofhow this noise can appear in the distance transform forthe Felix GOES-12 data set. A small area of pixels abovethe threshold cause the distance transform to fill in thecenter of the circular region, thus potentially changing thevalue of maximum persistence. Therefore, before apply-ing the distance transform, we use a method from mathe- Figure 8: Maximum persistence plotted over time for alldata sets using threshold µ = 80 in addition to the ver-sions using opening to remove noise. Gray vertical linesseparate days according to UTC.matical morphology ( [28]) called opening to de-noise theimage and see how this impacts the detected periodicity.Opening is the combination of two tools from mathe-matical morphology, erosion and dilation. Both involve akernel moving through a binary image. In erosion, a pixelin the original image will remain a 1 only if all pixels un-der the kernel are 1’s, otherwise it becomes a 0. Dilationis the opposite of erosion, A kernel moves through thebinary image and a pixel is assigned a 1 if at least onepixel under the kernel is a 1, otherwise it is assigned a 0.Opening is erosion followed by dilation, which will removenoise and rebuild the area around the boundary.We apply opening to the binary thresholded image us-ing a 8 × × cv2 for these computations. Opening is specifi-cally implemented using the function cv2.morphologyEx using cv2.MORPH OPEN as the second input. The right im-age in Fig. 7 show the result when opening is used on thethresholded matrix and then the distance transform is ap-plied. Since the distance transform is no longer filled in,the opening process has removed the noisy pixels causingthe issue.Using this extra step in the method, we recalculatemaximum persistence for all times and compute the es-timated period of the new maximum persistence valuesusing Fourier transforms. Figure 8 shows maximum per-7istence plotted versus time for using our original method,and the method including the additional opening step.While the new maximum persistence values vary a littlefrom the originals, the general oscillatory behavior seemssimilar. For both the Felix and Ivan GOES-12 data sets,the Fourier transform still detects a 24.5 and 24.0 hourcycle respectively. Thus the presence of noise in thesedata sets is not impacting the results. However, for theFelix GridSat-GOES data set, the Fourier transform nowdetects a 15.375 hour cycle, likely due to the differencein spatial resolution. The GOES-12 data has higher spa-tial resolution, so applying opening to remove noise doesnot impact the global circular structure. The GridSat-GOES data has lower spatial resolution, and is thereforemore sensitive to noise in the image. Thus, our methodis more reliable when applied to higher spatial resolutiondata, and should be used with caution on lower qualitydata. This paper presents a novel method for detecting and ana-lyzing the diurnal cycle of tropical cyclones using methodsfrom TDA. Current state of the art TC diurnal cycle mea-surement in the satellite imagery is mostly qualitative;our method provides a mathematically advanced methodfor automatic detection and measurement. While ourmethod involves a choice of a parameter for the thresh-old, we present evidence that a range of threshold valuesyield the same results.Here, we show that using two sets of GOES satellitedata for Hurricane Felix and one set of data for Hurri-cane Ivan, our method is able to detect almost identicalpatterns across all three data sets. Our method performsmore consistently and robustly on higher spatial resolu-tion data sets, represented by the GOES-12 data for bothHurricanes Felix and Ivan. While the method does detecta daily cycle in the lower resolution Felix GridSat-GOESdata, when the images are blurred to remove noise, thecycle is no longer detected by the Fourier transform.This is a novel application of methods from TDA andimage processing to the TC diurnal cycle. We believe thismethod could be used to study additional atmosphericphenomena exhibiting circular structure. A future direc-tion of this project is to apply this analysis to more TCs,other satellite channels and other atmospheric data tofurther test our method.
Acknowledgments
This manuscript is under consideration at Pattern Recog-nition Letters. The authors would like to thank WilliamDong, who was instrumental in creating code for the firstiteration of this project. The authors also thank thetwo anonymous reviewers whose suggestions improved thequality of the paper. ST and EM were supported in part by NSF grants DMS-1800446 and CMMI-1800466; EMwas also supported in part by NSF CCF-1907591. JD wassupported in part by the Office of Naval Research Pro-gram Element (PE)0601153N Grant N000141410132 and by NOAA GrantNA14OAR4830172 from the Unmanned Aircraft Systems(UAS) Program. KC was supported by NSF Grant AGS-1636799.
References [1] E. O. Brigham.
The Fast Fourier Transform andIts Applications . Prentice-Hall, Inc., Upper SaddleRiver, NJ, USA, 1988.[2] G. Carlsson, T. Ishkhanov, V. de Silva, andA. Zomorodian. On the local behavior of spaces ofnatural images.
International Journal of ComputerVision , 76:1–12, 2008.[3] A. Deckard, R. C. Anafi, J. B. Hogenesch, S. B.Haase, and J. Harer. Design and analysis of large-scale biological rhythm studies: a comparison of al-gorithms for detecting periodic signals in biologicaldata.
Bioinformatics , 29(24):3174–3180, 2013.[4] S. D. Ditchek, K. L. Corbosiero, R. G. Fovell, andJ. Molinari. Electrically active tropical cyclone di-urnal pulses in the atlantic basin.
Monthly WeatherReview , 147(10):3595–3607, 2019.[5] S. D. Ditchek, J. Molinari, K. L. Corbosiero, andR. G. Fovell. An objective climatology of tropicalcyclone diurnal pulses in the atlantic basin.
MonthlyWeather Review , 147(2):591–605, 2019.[6] J. P. Dunion, C. D. Thorncroft, and D. S. Nolan.Tropical cyclone diurnal cycle signals in a hurricanenature run.
Monthly Weather Review , 147(1):363–388, 2019.[7] J. P. Dunion, C. D. Thorncroft, and C. S. Velden.The tropical cyclone diurnal cycle of mature hurri-canes.
Monthly Weather Review , 142(10):3900–3919,2014.[8] H. Edelsbrunner and J. Harer.
Computational Topol-ogy: An Introduction . American Mathematical Soci-ety, 2010.[9] C. Giusti, E. Pastalkova, C. Curto, and V. Itskov.Clique topology reveals intrinsic geometric structurein neural correlations.
Proceedings of the NationalAcademy of Sciences , 2015.[10] A. Hatcher.
Algebraic Topology . Cambridge Univer-sity Press, 2002.[11] T. Kaczynski, K. Mischaikow, and M. Mrozek.
Com-putational homology , volume 157. Springer Science& Business Media, 2006.812] F. A. Khasawneh and E. Munch. Chatter detectionin turning using persistent homology.
MechanicalSystems and Signal Processing , 70-71:527–541, 2016.[13] K. R. Knapp and S. L. Wilkins. Gridded satel-lite (GridSat) GOES and CONUS data.
Earth Sys-tem Science Data , 10(3):1417–1425, aug 2018.[14] J. P. Kossin. Daily hurricane variability inferredfrom goes infrared imagery.
Monthly Weather Re-view , 130(9):2260–2270, 2002.[15] M. A. Lazzara, J. M. Benson, R. J. Fox, D. J.Laitsch, J. P. Rueden, D. A. Santek, D. M. Wade,T. M. Whittaker, and J. Young. The man computerinteractive data access system: 25 years of interac-tive processing.
Bulletin of the American Meteoro-logical Society , 80(2):271–284, 1999.[16] K. D. Leppert and D. J. Cecil. Tropical cyclone di-urnal cycle as observed by trmm.
Monthly weatherreview , 144(8):2793–2808, 2016.[17] K. Mischaikow and V. Nanda. Morse theory forfiltrations and efficient computation of persistenthomology.
Discrete & Computational Geometry ,50(2):330–353, 2013.[18] E. Munch. A user’s guide to topological data analy-sis.
Journal of Learning Analytics , 4(2), 2017.[19] V. Nanda. Perseus. , July 2015.[20] E. L. Navarro and G. J. Hakim. Idealized numericalmodeling of the diurnal cycle of tropical cyclones.
Journal of the Atmospheric Sciences , 73(10):4189–4201, 2016.[21] M. E. O’Neill, D. Perez-Betancourt, and A. A. Wing.Accessible environments for diurnal-period waves insimulated tropical cyclones.
Journal of the Atmo-spheric Sciences , 74(8):2489–2502, 2017.[22] J. A. Perea and J. Harer. Sliding windows and per-sistence: An application of topological methods tosignal analysis.
Foundations of Computational Math-ematics , pages 1–40, 2015.[23] V. Robins, M. Saadatfar, O. Delgado-Friedrichs, andA. P. Sheppard. Percolating length scales fromtopological persistence analysis of micro-CT imagesof porous materials.
Water Resources Research ,52(1):315–329, jan 2016.[24] V. Robins, P. J. Wood, and A. P. Sheppard. Theoryand algorithms for constructing discrete morse com-plexes from grayscale digital images.
IEEE Transac-tions on Pattern Analysis and Machine Intelligence ,33(8):1646–1658, Aug. 2011. [25] A. Rosenfeld and J. Pfaltz. Distance functions ondigital pictures.
Pattern Recognition , 1(1):33 – 61,1968.[26] A. Rosenfeld and J. L. Pfaltz. Sequential operationsin digital picture processing.
J. ACM , 13(4):471–494,Oct. 1966.[27] M. Saadatfar, H. Takeuchi, V. Robins, N. Francois,and Y. Hiraoka. Pore configuration landscape ofgranular crystallization.
Nature Communications ,8(15082), 2017.[28] J. Serra.
Image Analysis and Mathematical Morphol-ogy . Number v. 1 in Image Analysis and Mathemat-ical Morphology. Academic Press, 1984.[29] J. Steranka, E. B. Rodgers, and R. C. Gentry. Thediurnal variation of atlantic ocean tropical cycloneaoud distributioninferred from geostationary satel-lite infrared measurements.
Monthly weather review ,112(11):2338–2344, 1984.[30] X. Tang and F. Zhang. Impacts of the diurnal radia-tion cycle on the formation, intensity, and structureof hurricane edouard (2014).