Message Boards Message Boards

Hurricane Matthew and historical tropical storms

Background

Over the last couple of days hurricane Matthew has battered the East Coast of the United states. It has brought devastation and death to many communities. At the time of writing this post, the storm is still ongoing, but the Wolfram Language has already important data about the storm in its repository.

Grid[{CommonName[#], Entity["TropicalStorm", "Matthew2016"][#]} & /@ Evaluate[Entity["TropicalStorm", "Matthew2016"]["Properties"]], 
Frame -> All]

In this post I will present some very simple analyses and graphics to visualise some aspects of tropical storms. We will study historical paths of storms

enter image description here

and study some more of their characteristics.

Wolfram Language Data

The Wolfram Language holds an incredible amount of information about tropical storms. As of today he function

tropicalstorms = TropicalStormData["Entities"];

provides data of 12741 storms. The last 100 of which are

tropicalstorms[[-100 ;;]]

enter image description here

As you can see there are many storms in the list after Matthew. We will come back to some of that data later.

External data

For this post I will also need data on the position of the storm at different times - data which is not contained in the Wolfram Language. For Matthew this data can be found on this website. We can use

dataMatthew = 
Flatten[Import["https://www.wunderground.com/hurricane/atlantic/2016/Hurricane-Matthew?map=radar", "Data"][[3, 3, 1, 2 ;;]], 1];

to import and clean the data. It then looks a bit like this:

dataMatthew[[1 ;; 5]] // TableForm

enter image description here

Matthew's path

So now we have positions (3rd and 4th columns) and also windspeed (5th column in miles per hour). The max speed

Max[dataMatthew[[All, 5]]]

is 160 miles per hour so we can generate a coloured trail by normalising the windspeed.

GeoGraphics[{Yellow, Thickness[0.01], GeoPath[GeoPosition /@ dataMatthew[[All, {3, 4}]], VertexColors -> Evaluate[Hue[#/160.] & /@ dataMatthew[[All, 5]]]]},GeoBackground -> "Satellite"]

enter image description here

That image is still a bit static. I would like to display a photo of a (any) storm so that we can get a feeling for the movement. On google I found an image (which was marked as freely usable and distributable) - see attached. When I place that image on the map it has a squared shape. I did not like that - so I used the masking tool to create a mask (also attached.)

Importing the two

hurricatetemplate=Import["~/Desktop/HurricaneTemplate.jpeg"];
mask=Import["~/Desktop/HurricateTemplate.jpeg"];

we can create

hurrtemplate = RemoveBackground[ImageMultiply[hurrtemplate, mask]]

enter image description here

Now we can make a little video:

frames = 
Table[GeoGraphics[{Yellow, Thickness[0.01], GeoPath[GeoPosition /@ dataMatthew[[All, {3, 4}]], 
 VertexColors -> Evaluate[Hue[#/160.] & /@ dataMatthew[[All, 5]]]], Opacity[0.4], GeoMarker[GeoPosition[dataMatthew[[k, {3, 4}]]], hurrtemplate, 
 Scale -> 0.1]}, GeoBackground -> "Satellite", Epilog -> Inset[Text[Style[dataMatthew[[k, 1]] <> " " <> dataMatthew[[k, 2]], White,20]], {0.1, 0.2}]], {k, 1, Length[dataMatthew]}];

enter image description here

Ok. Now we can follow the path of Matthew. What about all the other storms?

More external data

I would need data for their paths, too. Luckily there is this website. If we download the data, we get a folder which contains several files. We will need to import two of them:

stormdata = Import["/Users/thiel/Desktop/Allstorms/Allstorms.ibtracs_all_lines.v03r09.shp"];
data = Import["/Users/thiel/Desktop/Allstorms/Allstorms.ibtracs_all_lines.v03r09.dbf"];

The shape file contains easily decipherable data about the paths. If we look at

(stormdata[[1, 2, 1, 2]][[1 ;; 32]])[[All, 1]]

We see a couple of geopositions for a couple of entries.

enter image description here

Preprocessing and cleaning the data

Unfortunately, they are not separated by storm. We can use the data file to separate the paths. Here is the first entry:

data[[1]]

enter image description here

The first entry is a unique identifier of the storm. It contains the year and day of the year encoded. The last five numbers give the year and time in a more standard form. It turns out that the 9th entry is the windspeed in knots. the 10th entry is the atmospheric pressure. Using the identifier for the name we can now separate the entries of each storm:

timeslength = Tally[data[[-150000 ;;, 1]]][[2 ;;]];

I only use the last 150000 entries. time length looks like this

timeslength[[1 ;; 3]]
{{"1967241N15170", 84}, {"1967241N16324", 28}, {"1967242N18253", 16}}

So it contains the identifier and the number of entries. We can now separate the storms:

separatedstorms = 
-Reverse /@ Join[{{1, (Reverse@timeslength[[All, 2]])[[1]]}}, {1, 0} + # & /@ Partition[Accumulate[Reverse@timeslength[[All, 2]]], 2, 1]];

Plotting larger numbers of storms

Plotting the paths of the last 1000 storms is not easy:

GeoGraphics[Join[{Green}, (GeoPath[GeoPosition[Reverse[#]] & /@ Flatten[(stormdata[[1, 2, 1, 2]][[#[[1]] ;; #[[2]]]])[[All, 1]], 1]]) & /@ separatedstorms[[1 ;; 1000]]], GeoRange -> "World", GeoBackground -> "Satellite", ImageSize -> Full]

enter image description here

The thing is that I would also like to represent the windspeed with colours along the paths. One needs to fiddle a bit with the speeds and normalisations but it is quite straight forward to come up with this:

nStorms = 300; 
GeoGraphics[(GeoPath[DeleteDuplicates[GeoPosition[Reverse[#]] & /@ 
Flatten[(stormdata[[1, 2, 1, 2]][[#[[1]] ;; #[[2]]]])[[All, 1]], 1]], 
VertexColors -> Evaluate[Hue[0.5 (1 - Sqrt[#])] & /@ ((#/Max[data[[separatedstorms[[nStorms, 1]] ;;, 9]] + 1]) & /@ ((data[[#[[1]] ;; #[[2]], 9]] /. {-999. -> -1}) + 1))]]) & /@ separatedstorms[[1 ;; nStorms]], GeoRange -> "World", GeoBackground -> "Satellite", ImageSize -> Full]

enter image description here

The nStorms variable allows us to choose the number of storms we want to visualise. nSTorms=1000 looks like this:

enter image description here

Geographical distribution of storms

We can clearly see different densities of storms in different regions and an equatorial gap. We now also have everything in place to visualise the density of storms in different regions. The new function GeoHistogram is quite useful for that:

nStorms3 = 3000; 
GeoHistogram[Flatten[DeleteDuplicates[GeoPosition[Reverse[{(Mod[180 + #[[1]], 360] - 180) /. -180. -> -179.9, #[[2]]}]] & /@ 
Flatten[(stormdata[[1, 2, 1, 2]][[#[[1]] ;; #[[2]]]])[[All, 1]],1]] & /@ separatedstorms[[1 ;; nStorms3]]], 
 Quantity[500, "Kilometers"], "PDF", PlotStyle -> Opacity[0.6], 
 ColorFunction -> "TemperatureMap", PerformanceGoal -> "Quality", 
 GeoRange -> "World", ImageSize -> Full, GeoBackground -> "Satellite"]

enter image description here

This gives some kind of birds eye view of different regions and their storm risk. I actually wanted to use different weights according to the wind speeds and GeoHistrogram allows us in principle to do that. I thought it should work like this:

nStorms2 = 2000; 
GeoHistogram[<|Rule @@@ Transpose@{Take[Flatten[DeleteDuplicates[GeoPosition[Reverse[#]] & /@ 
Flatten[(stormdata[[1, 2, 1, 2]][[#[[1]] ;; #[[2]]]])[[All,1]], 1]] & /@ separatedstorms[[1 ;; nStorms2]]], 
UpTo[Length[Flatten[(data[[#[[1]] ;; #[[2]], 9]]) & /@ 
separatedstorms[[1 ;; nStorms2]]]]]], Flatten[(#/Max[data[[separatedstorms[[nStorms2, 1]] ;;, 9]] + 1]) & /@ ((data[[#[[1]] ;; #[[2]], 9]] /. {-999. -> -1}) + 1) & /@ 
separatedstorms[[1 ;; nStorms2]]]}|>, 
Quantity[500, "Kilometers"], "Intensity", 
GeoBackground -> "Satellite", PlotStyle -> Opacity[0.3], 
ColorFunction -> "TemperatureMap", ImageSize -> Full]

But this obviously goes wrong:

enter image description here

I made a couple of observations of why this goes wrong and what causes it, but this is not the place to discuss this at length.

Northern vs southern hemisphere

Another interesting thing to plot is the paths of the storms for those on the northern hemisphere and those on the southern hemisphere. I will normalise all storms so that they all start at the origin. Here is the plot for north:

Monitor[ListLinePlot[
  Table[Reverse /@ (({#[[1]], Mod[#[[2]], 360]} - {Evaluate[Select[paths2, #[[1, 1]] > 0 &]][[k, 1]][[1]], 
  Mod[Evaluate[Select[paths2, #[[1, 1]] > 0 &]][[k, 1]][[2]], 360]}) & /@ Evaluate[Select[paths2, #[[1, 1]] > 0 &][[k]]]), {k, 1, 100}], 
  PlotRange -> All, Axes -> False, Background -> Black], k]

enter image description here

and there the same for south:

Monitor[ListLinePlot[
 Table[Reverse /@ (({#[[1]], Mod[#[[2]], 360]} - {Evaluate[Select[paths2, #[[1, 1]] < 0 &]][[k, 1]][[1]], 
 Mod[Evaluate[Select[paths2, #[[1, 1]] < 0 &]][[k, 1]][[2]], 360]}) & /@ 
 Evaluate[Select[paths2, #[[1, 1]] < 0 &][[k]]]), {k, 1, 100}], 
 PlotRange -> All, Axes -> False, Background -> Black], k]

enter image description here

The colours are quite meaningless, but we can draw everything together and obtain:

Monitor[southernstorms = 
  ListLinePlot[
   Table[Reverse /@ (({#[[1]], Mod[#[[2]], 360]} - {Evaluate[Select[paths2, #[[1, 1]] < 0 &]][[k, 1]][[1]], 
   Mod[Evaluate[Select[paths2, #[[1, 1]] < 0 &]][[k, 1]][[2]],360]}) & /@ 
   Evaluate[Select[paths2, #[[1, 1]] < 0 &][[k]]]), {k, 1, 100}], 
   PlotRange -> All, Axes -> False, Background -> Black, 
   PlotStyle -> Directive[Yellow, Thickness[0.001]]], k]

and

Monitor[northernStorms = 
  ListLinePlot[Table[Reverse /@ (({#[[1]], Mod[#[[2]], 360]} - {Evaluate[Select[paths2, #[[1, 1]] > 0 &]][[k, 1]][[1]], 
  Mod[Evaluate[Select[paths2, #[[1, 1]] > 0 &]][[k, 1]][[2]],360]}) & /@ 
   Evaluate[Select[paths2, #[[1, 1]] > 0 &][[k]]]), {k, 1, 100}], 
   PlotRange -> All, Axes -> False, Background -> Black, 
   PlotStyle -> Directive[Red, Thickness[0.001]]], k]

to give:

Show[southernstorms, northernStorms, PlotRange -> All, ImageSize -> Full]

enter image description here

One can clearly see the different directions of movement in the north and the south.

Maximal windspeeds of storms

Another thing we can try is look for the maximal windspeeds. First we calculate the max windspeed for each storm:

maxStormspeed = 
 DeleteCases[Max /@ (((data[[#[[1]] ;; #[[2]], 9]] /. {-999. -> -1}) + 1) & /@ separatedstorms[[1 ;; 2000]]), 0] ;

and then plot the corresponding histogram:

Plot[PDF[SmoothKernelDistribution[maxStormspeed], x], {x, 0, 320}, Filling -> Bottom, PlotTheme -> "Marketing", 
FrameLabel -> {"Windspeed", "Probability"}, LabelStyle -> Directive[Bold, Medium]]

enter image description here

So it is clear that every strong storms are very rare - as expected. We can plot that logarithmically to emphasise this result:

Plot[PDF[SmoothKernelDistribution[maxStormspeed], x], {x, 0, 320}, 
 Filling -> Bottom, PlotTheme -> "Marketing", 
 ScalingFunctions -> {"Linear", "Log"}, 
 FrameLabel -> {"Windspeed", "Log Probability"}, 
 LabelStyle -> Directive[Bold, Medium]]

enter image description here

Note that the windspeed is in knots. One of the strongest storms ever recorded appears to be Camille of 1969 with windspeeds up to 190 miles per hour. The corresponds to:

N[Quantity[1911096/11575, "Knots"]]

165 knots, which is more or less where this curve ebbs off. Using data from the Wolfram Language directly, we obtain a very similar picture:

mmawindspeeds = #[EntityProperty["TropicalStorm", "WindSpeed"]] & /@ tropicalstorms[[-500 ;;]];
Plot[PDF[SmoothKernelDistribution[
QuantityMagnitude /@ DeleteMissing[mmawindspeeds]], x], {x, 0, 320}, 
 Filling -> Bottom, PlotTheme -> "Marketing", 
 FrameLabel -> {"Windspeed", "Probability"}, 
 LabelStyle -> Directive[Bold, Medium]]

enter image description here

If we wish to do so, we can also ask Mathematica to "guess" the mathematical form of the distribution.The four most likely candidates are:

TableForm[FindDistribution[QuantityMagnitude /@ DeleteMissing[mmawindspeeds], 4]]

enter image description here

We can look at the quality of the different fits:

GraphicsGrid[
 Partition[
  Table[Plot[{PDF[
      SmoothKernelDistribution[
       QuantityMagnitude /@ DeleteMissing[mmawindspeeds]], x], 
     PDF[distribs[[k]], x]}, {x, 0, 320}, Filling -> Bottom, 
    PlotTheme -> "Marketing", 
    FrameLabel -> {"Windspeed", "Probability"}, 
    LabelStyle -> Directive[Bold, Medium], PlotRange -> All], {k, 1, 
    4}], 2]]

enter image description here

The discussion about windspeed distribution has come up before in this community.

How are storms called?

Using data built into the Wolfram Language we can also look for common names that are given to these storms:

WordCloud[Select[DeleteCases[CommonName[tropicalstorms], "Unnamed"], StringFreeQ[#, CharacterRange["0", "9"]] &]]

enter image description here

Back to Matthew

Regarding Matthew, I tried to look at the wind flow today. Mathematica has the data:

Monitor[windvectordata = 
Table[{GeoPosition[{i, j}], WindVectorData[GeoPosition[{i, j}], 
DateObject[{2016, 10, 9}, TimeObject[{6, 0}]]]}, {i, 12, 50, 2}, {j, -90, -55.8, 2}];, {i, j}]

which we can then plot:

ls = ListStreamPlot[
Cases[{Flatten[Reverse[#[[1, 1]]]], QuantityMagnitude[#[[2]]]} & /@
Flatten[windvectordata, 1], {_, {_?NumberQ, _?NumberQ}}]];
arrows = Cases[ls, Arrow[_], Infinity];
GeoGraphics[{Yellow, Arrowheads[Small], arrows}, 
 GeoBackground -> "Satellite"]

enter image description here

The eye of the storm becomes quite obvious but even in a ListVectorDensityPlot

ls2 = ListVectorDensityPlot[
  Cases[{Flatten[Reverse[#[[1, 1]]]], QuantityMagnitude[#[[2]]]} & /@ 
    Flatten[windvectordata, 1], {_, {_?NumberQ, _?NumberQ}}], 
  ColorFunction -> "TemperatureMap", StreamColorFunction -> Hue, 
  StreamPoints -> Fine, VectorPoints -> 0, Frame -> None, 
  ImagePadding -> False]

The current situation does not seem to be quite well reflected. Here is a screenshot from another program I have:

enter image description here

These images do not represent the exact same time, but there should be another storm out at the see. I can only suppose that this is because WindVectorData uses interpolation between weather stations, but I am not sure about this.

I will probably adds some more bits to this post over the coming days, so stay tuned.

Cheers,

Marco

Attachments:
POSTED BY: Marco Thiel
3 Replies

Marvelous exploration of tropical storms, @Marco, thank you for sharing! Normalized storm trajectory bundles for different hemispheres is a neat idea.

POSTED BY: Vitaliy Kaurov

Very nice analysis! Thanks for sharing! I had no idea this data was available :)

POSTED BY: Sander Huisman

enter image description here - another post of yours has been selected for the Staff Picks group, congratulations! We are happy to see you at the top of the "Featured Contributor" board. Thank you for your wonderful contributions, and please keep them coming!

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract