Annual Quantiles of Tornado Energy Dissipation

Get packages.

suppressMessages(library("tidyverse"))
suppressMessages(library("lubridate"))
suppressMessages(library("rgeos"))
suppressMessages(library("sf"))

Get tornado data. Set the domain and grid resolution. Data from the Storm Prediction Center (SPC) http://www.spc.noaa.gov/gis/svrgis/.

if(!file.exists("torn")){ 
    download.file(url = "http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
              destfile = "tornado.zip")
    unzip("tornado.zip")
  }

Load and prepare the tornado data.

Tor.sfdf <- sf::st_read(dsn = "torn", 
                      layer = "torn", 
                      stringsAsFactors = FALSE)

## Reading layer `torn' from data source `/Users/jelsner/Dropbox/Tornadoes/torn' using driver `ESRI Shapefile'
## Simple feature collection with 61092 features and 22 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -6279240 ymin: -1835231 xmax: 3405399 ymax: 3888106
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

Tor.sfdf <- Tor.sfdf %>%
  filter(yr >= 1970, mag != -9) %>%
  mutate(DateTime = as.POSIXct(paste(yr, mo, dy, time), format = "%Y%m%d%H:%M:%S"),
         Year = year(DateTime),
         Length = len * 1609.34,
         Width = wid * .9144,
         Width = ifelse(Width == 0, min(Width[Width > 0]), Width), #takes care of zero width
         Width = ifelse(Year >= 1995, Width * pi/4, Width), #takes care of change: avg to max
         cas = inj + fat,
         AreaPath = Length * Width,
         Ma = factor(month.abb[mo], levels = month.abb[1:12])) %>%
  sf::st_sf()
#TornL <- as(Tor.sfdf, "Spatial")

Compute energy dissipation.

perc <- c(1, 0, 0, 0, 0, 0, 
         .772, .228, 0, 0, 0, 0,
         .616, .268, .115, 0, 0, 0,
         .529, .271, .133, .067, 0, 0,
         .543, .238, .131, .056, .032, 0,
         .538, .223, .119, .07, .033, .017)
percM <- matrix(perc, ncol = 6, byrow = TRUE)
threshW <- c(29.06, 38.45, 49.62, 60.8, 74.21, 89.41)
midptW <- c(diff(threshW)/2 + threshW[-length(threshW)], threshW[length(threshW)] + 7.5)
ef <- Tor.sfdf$mag + 1
EW3 <- numeric()
for(i in 1:length(ef)) EW3[i] <- midptW^3 %*% percM[ef[i], ]
Tor.sfdf$e3 <- EW3
Tor.sfdf$ED <- Tor.sfdf$e3 * Tor.sfdf$AreaPath
mean(Tor.sfdf$ED)/10^9  # ED in units of gigawatts

## [1] 89.00028

Energy dissipation by EF rating.

as.data.frame(Tor.sfdf) %>%
  group_by(mag) %>%
  dplyr::summarize(mED = mean(ED)/10^9,
                   tED = sum(ED)/10^9)

## # A tibble: 6 x 3
##     mag     mED      tED
##   <int>   <dbl>    <dbl>
## 1     0    3.81   97527.
## 2     1   33.0   532872.
## 3     2  183.   1048498.
## 4     3  885.   1422205.
## 5     4 3051.   1095357.
## 6     5 5703.    205317.

Annual quantiles

df <- as.data.frame(Tor.sfdf) %>%
  filter(yr >= 1970) %>%
  group_by(yr) %>%
  dplyr::summarize(q25ED = quantile(ED, prob = .25),
                   q75ED = quantile(ED, prob = .75),
                   q90ED = quantile(ED, prob = .9),
                   q50ED = quantile(ED, prob = .5))

ggplot(df, aes(x = yr, y = q50ED/10^9)) +
  geom_errorbar(aes(ymin = q25ED/10^9, ymax = q75ED/10^9), color = "grey", width = .2) +
  geom_point() +
  geom_point(aes(x = yr, y = q90ED/10^9), data = df, color = "red") +
  scale_y_log10() +
  xlab("Year") + ylab("Energy Dissipation [GW]") +
  theme_minimal() +
  ggtitle("Annual Quantiles of Tornado Energy Dissipation",
          subtitle = "Median, Interquartile Range, 90th Percentile")