Code in support of our paper submitted to the Electronic Journal of Severe Storms Meteorology: http://www.ejssm.org/
Holly Widen, James B. Elsner, Cameron Amrine, Rizalino Cruz, Erik Fraza, Laura Michaels, Loury Migliorelli, Brendan Mulholland, Michael Patterson, Sarah Strazzo, Guang Xing
date()
## [1] "Sun Aug 04 15:38:48 2013"
Set working directory and load packages.
# setwd('~/Dropbox/Tornadoes/Holly')
library(rgdal)
library(plyr)
library(rgeos)
library(xtable)
library(spatstat)
library(maptools)
library(ggplot2)
Downloaded from SPC at http://www.spc.noaa.gov/gis/svrgis/
tmp = download.file("http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
"tornado.zip", mode = "wb")
unzip("tornado.zip")
Torn = readOGR(dsn = ".", layer = "tornado")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "tornado"
## with 56221 features and 28 fields
## Feature type: wkbPoint with 2 dimensions
Tornadoes through 2012. CSV file only. State segments. NOT RUN.
con = "http://www.spc.noaa.gov/wcm/data/1950-2012_torn.csv"
TornCSV = read.csv(con, header = FALSE)
cn = c("OM", "YEAR", "MONTH", "DAY", "DATE", "TIME", "TIMEZONE", "STATE", "FIPS",
"STATENUMBE", "FSCALE", "INJURIES", "FATALITIES", "LOSS", "CROPLOSS", "SLAT",
"SLON", "ELAT", "ELON", "LENGTH", "WIDTH", "NS", "SN", "SG", "F1", "F2",
"F3", "F4")
names(TornCSV) = cn
TornA = TornA[-anyDuplicated(TornA), ]
TornA$ELON = -abs(TornA$ELON)
Torn2 = subset(Torn)
Torn2$Category = factor(Torn2$FSCALE, ordered = TRUE)
Torn2$Length = Torn2$LENGTH * 1609.34
Torn2$Width = Torn2$WIDTH * 0.9144
Torn2$Area = Torn2$Width * Torn2$Length
Illustrate method using statewide tornado data of KS
tmp = download.file("http://www2.census.gov/geo/tiger/PREVGENZ/st/st00shp/st20_d00_shp.zip",
"st20_d00_shp.zip", mode = "wb")
unzip("st20_d00_shp.zip")
STATE = readOGR(dsn = ".", layer = "st20_d00")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st20_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
ll = "+proj=longlat +ellps=GRS80 +datum=NAD83"
proj4string(STATE) = CRS(ll)
STATE = spTransform(STATE, CRS(proj4string(Torn2)))
stArea = gArea(STATE)
st = over(Torn2, STATE)
st = !is.na(st$AREA)
TornA = Torn2[st, ]
TornS = subset(TornA, FSCALE >= 2)
TornV = subset(TornA, FSCALE >= 4)
nYears = length(1950:2011)
AnnualDensityA = dim(TornA@data)[1]/stArea * 1e+10/nYears
AnnualDensityS = dim(TornS@data)[1]/stArea * 1e+10/nYears
There are 3713 Kansas tornado reports over the period 1950–2011, inclusive. Table~\ref{tab:KSFscaleDist} shows the number of tornadoes by F-scale category. A total 266 (r I(table(TornA@data$FSCALE)[1]/dim(TornA@data)[1]*100) \%) tornado reports have a missing F-scale category.
Table 1:
table(TornA@data$FSCALE)
##
## -9 0 1 2 3 4 5
## 266 1969 862 419 156 35 6
table(TornA@data$FSCALE[TornA@data$YEAR >= 1991 & TornA@data$YEAR <= 2010])
##
## 0 1 2 3 4 5
## 1370 387 92 43 15 2
Quantile path lengths and widths by FSCALE:
qLW = ddply(TornA@data, .(FSCALE), summarize, qL = quantile(Length), qW = quantile(Width))
The median path length is 853 m and the median path width is 46 m. The median path area is 4.4 × 104 m\( ^2 \).
On average stronger tornadoes travel longer distances. In Kansas, the subset of significant tornadoes (F2 or higher), has a median path length is 9576 m and the median width is 91 m. The median path area is 6.5338 × 105 m\( ^2 \).
There are 41 violent tornadoes (F4–F5) (1.1042% of the total number of reports) in Kansas over the 62-year period. The median path length and width are 3.8624 × 104 and 457 m, respectively. The median path area is 1.8395 × 107 m\( ^2 \).
cor.test(TornA@data$Length, TornA@data$Width)
##
## Pearson's product-moment correlation
##
## data: TornA@data$Length and TornA@data$Width
## t = 22.23, df = 3711, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3141 0.3709
## sample estimates:
## cor
## 0.3428
cor.test(TornS@data$Length, TornS@data$Width)
##
## Pearson's product-moment correlation
##
## data: TornS@data$Length and TornS@data$Width
## t = 6.298, df = 614, p-value = 5.741e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1707 0.3191
## sample estimates:
## cor
## 0.2463
cor.test(TornV@data$Length, TornV@data$Width)
##
## Pearson's product-moment correlation
##
## data: TornV@data$Length and TornV@data$Width
## t = 1.336, df = 39, p-value = 0.1893
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1052 0.4856
## sample estimates:
## cor
## 0.2092
cor(TornA@data$Area, TornA@data$FSCALE, method = "s")
cor(TornS@data$Area, TornS@data$FSCALE, method = "s")
cor(TornV@data$Area, TornV@data$FSCALE, method = "s")
The annual statewide density of Kansas is 2.84 reports per 10\( ^4 \) square km. For the set of significant tornadoes the density is 0.47 reports per 10\( ^4 \) square km.
With the exception of the EF-4 to EF-5 reports there is no obvious spatial trend in the report densities.
par(mfrow = c(2, 2), mex = 0.9, mgp = c(2, 0.4, 0))
plot(STATE, col = "grey90")
plot(TornA, add = TRUE, pch = 16, cex = 0.5, col = "red")
mtext("a", side = 3, line = 1, adj = 0, cex = 1)
plot(STATE, col = "grey90")
plot(TornS, add = TRUE, pch = 16, cex = 0.5, col = "red")
mtext("b", side = 3, line = 1, adj = 0, cex = 1)
plot(STATE, col = "grey90")
plot(TornA[TornA$FSCALE >= 3, ], add = TRUE, pch = 16, cex = 0.5, col = "red")
mtext("c", side = 3, line = 1, adj = 0, cex = 1)
plot(STATE, col = "grey90")
plot(TornV, add = TRUE, pch = 16, cex = 0.5, col = "red")
mtext("d", side = 3, line = 1, adj = 0, cex = 1)
Figure 1: Kansas tornado reports over the period 1950–2011. The points are the tornado touchdown locations for (a) All, (b) EF-2 to EF-5, © EF-3 to EF-5, and (d) EF-4 to EF-5. The set of all tornado reports includes those without an F-scale category.
Download and overlay city data from www.nws.noaa.gov/geodata/catalog/national/data/
tmp = download.file("http://www.nws.noaa.gov/geodata/catalog/national/data/ci08au12.zip",
"ci08au12.zip", mode = "wb")
unzip("ci08au12.zip")
US = readOGR(dsn = ".", layer = "ci08au12", p4s = ll)
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "ci08au12"
## with 43603 features and 23 fields
## Feature type: wkbPoint with 2 dimensions
UST = spTransform(US, CRS(proj4string(TornA)))
st = over(UST, STATE)
st = !is.na(st$AREA)
Cities = UST[st, ]
There are 871 Kansas cities in the U.S. cities database.
par(mfrow = c(1, 1))
plot(STATE, col = "grey90")
plot(TornA, add = TRUE, pch = 16, cex = 0.5, col = "red")
plot(Cities, pch = 15, col = "black", cex = 0.5, add = TRUE)
Figure 2: Kansas tornado reports and city centers.
Create PPP objects from SPDFs
W = as(STATE, "owin")
## Checking 1 polygon...[Checking polygon with 1329 edges...]
## done.
All.ppp = as.ppp(TornA)
All.ppp = All.ppp[W]
S.ppp = as.ppp(TornS)
S.ppp = S.ppp[W]
C = as.ppp(Cities)
window = as.owin(All.ppp)
C = C[window]
Zc = distmap(C)
rhatA = rhohat(unmark(All.ppp), Zc, bw = 0.25, smoother = "kernel", method = "transform",
adjust = 1)
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
m2km = 0.001
km2 = 1e+10
dist = rhatA$Zc * m2km
rho = rhatA$rho * km2/nYears
hi = rhatA$hi * km2/nYears
lo = rhatA$lo * km2/nYears
rhatA.df = data.frame(dist = dist, rho = rho, hi = hi, lo = lo)
rhatS = rhohat(unmark(S.ppp), Zc, bw = 0.25, smoother = "kernel", method = "transform",
adjust = 1)
dist = rhatS$Zc * m2km
rho = rhatS$rho * km2/nYears
hi = rhatS$hi * km2/nYears
lo = rhatS$lo * km2/nYears
rhatS.df = data.frame(dist = dist, rho = rho, hi = hi, lo = lo)
f = fivenum(Zc$v) * m2km
A statistical description of report clusters near cities and towns is obtained by estimating the spatial report density as a smoothed function of distance from nearest city center. First a 128 \( \times \) 128 grid containing distances from nearest city center is computed. Distances range from a minimum of 0.04 km to a maximum of 33.7 km with a median of 7 km. Fifty percent of all pixels are between 4.6 and 10.3 km from the nearest city.
source("multiplot.txt")
p3a = ggplot(rhatA.df) + geom_ribbon(aes(x = dist, ymin = lo, ymax = hi), alpha = 0.3) +
geom_line(aes(x = dist, y = rho), col = "black") + ylab(expression(paste(hat(rho),
"(Annual Reports/10", {
}^4, " km", {
}^2, ")"))) + xlab("Distance from Nearest City Center (km)") + theme_bw() +
ggtitle("a") + theme(plot.title = element_text(hjust = 0))
p3b = ggplot(rhatS.df) + geom_ribbon(aes(x = dist, ymin = lo, ymax = hi), alpha = 0.3) +
geom_line(aes(x = dist, y = rho), col = "black") + ylab(expression(paste(hat(rho),
"(Annual Reports/10", {
}^4, " km", {
}^2, ")"))) + xlab("Distance from Nearest City Center (km)") + theme_bw() +
ggtitle("b") + theme(plot.title = element_text(hjust = 0))
multiplot(p3a, p3b, cols = 1)
Figure 3: Tornado report density as a function of distance from nearest city center for (a) all tornadoes and (b) strong and violent tornadoes (EF-2 or higher). The gray band is the 95\% CI on the density estimate.
Fit a non-stationary Poisson process model for hazard probability. Or use a non-parametric kernel density smoother.
ArA = sum(TornA$Area)
ArS = sum(TornS$Area)
fitA = rhohat(unmark(All.ppp), Zc, bw = 0.25, smoother = "kernel", method = "transform",
adjust = 1)
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
fitS = rhohat(unmark(S.ppp), Zc, bw = 0.25, smoother = "kernel", method = "transform",
adjust = 1)
beta0A = fitA$rho[1]
beta0S = fitS$rho[1]
beta0LoA = fitA$lo[1]
beta0HiA = fitA$hi[1]
beta0LoS = fitS$lo[1]
beta0HiS = fitS$hi[1]
alphaA = beta0A * stArea/All.ppp$n
alphaS = beta0S * stArea/S.ppp$n
alphaLoA = beta0LoA * stArea/All.ppp$n
alphaHiA = beta0HiA * stArea/All.ppp$n
alphaLoS = beta0LoS * stArea/S.ppp$n
alphaHiS = beta0HiS * stArea/S.ppp$n
AcA = alphaA * ArA
AcS = alphaS * ArS
AcLoA = alphaLoA * ArA
AcHiA = alphaHiA * ArA
AcLoS = alphaLoS * ArS
AcHiS = alphaHiS * ArS
PA = AcA/(stArea * nYears)
PLoA = AcLoA/(stArea * nYears)
PHiA = AcHiA/(stArea * nYears)
PS = AcS/(stArea * nYears)
PLoS = AcLoS/(stArea * nYears)
PHiS = AcHiS/(stArea * nYears)
PLoA * 100
## [1] 0.06331
PA * 100
## [1] 0.06609
PHiA * 100
## [1] 0.06886
PLoS * 100
## [1] 0.05278
PS * 100
## [1] 0.05862
PHiS * 100
## [1] 0.06447
The method results in a statewide hazard probability of 0.0661 (0.0633, 0.0689) \% (95\% CI) per year. The estimate is the annual probability of getting hit by a tornado at any location in Kansas.
Apply method to 15 other states (16 total) using data at http://www.census.gov/geo/www/cob/st2000.html#shp
Kansas (20) Iowa (19) Oklahoma (40) Nebraska (31) Missouri (29) Arkansas (05) Louisiana (22) Illinois (17) Indiana (18) Kentucky (21) Tennessee (47) Mississippi (28) Alabama (01) Ohio (39) Minnesota (27) Wisconsin (55)
Initialize four data frames for storing results. Download state boundary files.
Final0 = numeric()
Final1 = numeric()
ll = "+proj=longlat +ellps=GRS80 +datum=NAD83"
stNo = c("01", "05", "17", "18", "19", "20", "21", "22", "27", "28", "29", "31",
"39", "40", "47", "55")
stName = c("Alabama", "Arkansas", "Illinois", "Indiana", "Iowa", "Kansas", "Kentucky",
"Louisiana", "Minnesota", "Mississippi", "Missouri", "Nebraska", "Ohio",
"Oklahoma", "Tennessee", "Wisconsin")
stAbb = c("AL", "AR", "IL", "IN", "IA", "KS", "KY", "LA", "MN", "MS", "MO",
"NE", "OH", "OK", "TN", "WI")
plotP = list()
beginYear = 1950
endYear = 2011
DT = 1
Loop over other states
for (j in 1:16) {
print(stName[j])
lay = paste("st", stNo[j], "_d00", sep = "")
STATE = readOGR(dsn = ".", layer = lay)
proj4string(STATE) = CRS(ll)
STATE = spTransform(STATE, CRS(proj4string(Torn2)))
stArea = gArea(STATE)
st = over(Torn2, STATE)
st = !is.na(st$AREA)
TornA = Torn2[st, ]
tt = as.vector(table(TornA$FSCALE))
df0 = data.frame(State = stName[j], F0 = tt[2], F1 = tt[3], F2 = tt[4],
F3 = tt[5], F4 = tt[6], F5 = tt[7], M = tt[1])
TornS = subset(TornA, FSCALE >= 2)
TornV = subset(TornA, FSCALE >= 4)
nYears = length(beginYear:endYear)
AnnualDensityA = dim(TornA@data)[1]/stArea * 1e+10/nYears
AnnualDensityS = dim(TornS@data)[1]/stArea * 1e+10/nYears
st = over(UST, STATE)
st = !is.na(st$AREA)
Cities = UST[st, ]
W = as(STATE, "owin")
All.ppp = as.ppp(TornA)
All.ppp = All.ppp[W]
S.ppp = as.ppp(TornS)
S.ppp = S.ppp[W]
C = as.ppp(Cities)
C = C[W]
Zc = distmap(C)
ArA = sum(TornA$Area)
ArS = sum(TornS$Area)
fitA = rhohat(unmark(All.ppp), Zc, bw = 0.25, smoother = "kernel", method = "transform",
adjust = 1)
fitS = rhohat(unmark(S.ppp), Zc, bw = 0.25, smoother = "kernel", method = "transform",
adjust = 1)
beta0A = fitA$rho[1]
beta0S = fitS$rho[1]
beta0LoA = fitA$lo[1]
beta0HiA = fitA$hi[1]
beta0LoS = fitS$lo[1]
beta0HiS = fitS$hi[1]
alphaA = beta0A * stArea/All.ppp$n
alphaS = beta0S * stArea/S.ppp$n
alphaLoA = beta0LoA * stArea/All.ppp$n
alphaHiA = beta0HiA * stArea/All.ppp$n
alphaLoS = beta0LoS * stArea/S.ppp$n
alphaHiS = beta0HiS * stArea/S.ppp$n
AcA = alphaA * ArA
AcS = alphaS * ArS
AcLoA = alphaLoA * ArA
AcHiA = alphaHiA * ArA
AcLoS = alphaLoS * ArS
AcHiS = alphaHiS * ArS
PA = AcA/(stArea * nYears) * 100
PLoA = AcLoA/(stArea * nYears) * 100
PHiA = AcHiA/(stArea * nYears) * 100
PS = AcS/(stArea * nYears) * 100
PLoS = AcLoS/(stArea * nYears) * 100
PHiS = AcHiS/(stArea * nYears) * 100
df1 = data.frame(State = stName[j], PA = PA, PLoA = PLoA, PHiA = PHiA, PS = PS,
PLoS = PLoS, PHiS = PHiS)
Final0 = rbind(Final0, df0)
Final1 = rbind(Final1, df1)
}
## [1] "Alabama"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st01_d00"
## with 2 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 2 polygons...1, [Checking polygon with 1684 edges...]
## 2.
## done.
## Checking for cross-intersection between 2 polygons... 1.
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: 6 out of 45 pixel values were outside the pixel image domain and
## were estimated by convolution
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: 1 out of 9 pixel values were outside the pixel image domain and
## were estimated by convolution
## [1] "Arkansas"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st05_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 2488 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Illinois"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st17_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 2832 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Indiana"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st18_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 1782 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Iowa"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st19_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 2719 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Kansas"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st20_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 1329 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Kentucky"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st21_d00"
## with 2 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 2 polygons...1, [Checking polygon with 3139 edges...]
## 2.
## done.
## Checking for cross-intersection between 2 polygons... 1.
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Louisiana"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st22_d00"
## with 2 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 2 polygons...1, [Checking polygon with 3167 edges...]
## 2.
## done.
## Checking for cross-intersection between 2 polygons... 1.
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Minnesota"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st27_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 4574 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Mississippi"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st28_d00"
## with 5 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 5 polygons...1, [Checking polygon with 2817 edges...]
## 2, 3, 4, 5.
## done.
## Checking for cross-intersection between 5 polygons...1, 2, 3, 4.
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Missouri"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st29_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 3005 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Nebraska"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st31_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 1948 edges...]
## done.
## [1] "Ohio"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st39_d00"
## with 5 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 5 polygons...1, [Checking polygon with 2147 edges...]
## 2, 3, 4, 5.
## done.
## Checking for cross-intersection between 5 polygons...1, 2, 3, 4.
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Oklahoma"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st40_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 2171 edges...]
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
## [1] "Tennessee"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st47_d00"
## with 1 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 1 polygon...[Checking polygon with 2071 edges...]
## done.
## [1] "Wisconsin"
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "st55_d00"
## with 6 features and 10 fields
## Feature type: wkbPolygon with 2 dimensions
## Checking 6 polygons...1, 2, 3, [Checking polygon with 3258 edges...]
## 4, 5, 6.
## done.
## Checking for cross-intersection between 6 polygons...1, 2, 3, 4, 5.
## done.
## Warning: Some tiles with zero area contain points
## Warning: Some weights are zero
Creating resulting tables:
Final0[is.na(Final0)] = 0
Final0$Total = rowSums(Final0[, 2:8])
Final0a = rename(Final0, c(M = "Unknown"))
Final0b = Final0a[order(Final0a$Total, decreasing = TRUE), ]
xtable(Final0b, digits = 0) #Table 2
## % latex table generated in R 3.0.1 by xtable 1.7-1 package
## % Sun Aug 04 15:42:01 2013
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrrrrr}
## \hline
## & State & F0 & F1 & F2 & F3 & F4 & F5 & Unknown & Total \\
## \hline
## 6 & Kansas & 1969 & 862 & 419 & 156 & 35 & 6 & 266 & 3713 \\
## 14 & Oklahoma & 1372 & 1067 & 655 & 184 & 48 & 7 & 73 & 3406 \\
## 12 & Nebraska & 1255 & 728 & 281 & 76 & 23 & 1 & 219 & 2583 \\
## 5 & Iowa & 928 & 673 & 414 & 100 & 39 & 6 & 84 & 2244 \\
## 3 & Illinois & 977 & 659 & 353 & 97 & 23 & 2 & 41 & 2152 \\
## 11 & Missouri & 741 & 787 & 325 & 92 & 33 & 1 & 12 & 1991 \\
## 10 & Mississippi & 524 & 736 & 385 & 129 & 24 & 4 & 59 & 1861 \\
## 1 & Alabama & 576 & 678 & 385 & 123 & 35 & 7 & 7 & 1811 \\
## 8 & Louisiana & 542 & 795 & 293 & 87 & 9 & 1 & 1 & 1728 \\
## 2 & Arkansas & 451 & 602 & 389 & 142 & 27 & 0 & 15 & 1626 \\
## 9 & Minnesota & 819 & 501 & 189 & 50 & 19 & 2 & 10 & 1590 \\
## 4 & Indiana & 393 & 466 & 269 & 82 & 25 & 2 & 40 & 1277 \\
## 16 & Wisconsin & 404 & 483 & 255 & 47 & 17 & 3 & 38 & 1247 \\
## 15 & Tennessee & 434 & 225 & 82 & 25 & 1 & 0 & 282 & 1049 \\
## 13 & Ohio & 293 & 390 & 181 & 41 & 14 & 3 & 14 & 936 \\
## 7 & Kentucky & 321 & 173 & 64 & 16 & 1 & 0 & 210 & 785 \\
## \hline
## \end{tabular}
## \end{table}
Total number of tornadoes in each state: 1811, 1626, 2152, 1277, 2244, 3713, 785, 1728, 1590, 1861, 1991, 2583, 936, 3406, 1049, 1247 Number of strong tornadoes in each state: 550, 558, 475, 378, 559, 616, 81, 390, 260, 542, 451, 381, 239, 894, 108, 322 Number of violent tornadoes in each state: 42, 27, 25, 27, 45, 41, 1, 10, 21, 28, 34, 24, 17, 55, 1, 20
Final1a = Final1[order(Final1$PA, decreasing = TRUE), ] #order for ranking
xtable(Final1a, digits = 3) #Table 3 without added SS2011 stats
## % latex table generated in R 3.0.1 by xtable 1.7-1 package
## % Sun Aug 04 15:42:01 2013
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrrrrr}
## \hline
## & State & PA & PLoA & PHiA & PS & PLoS & PHiS \\
## \hline
## 1 & Alabama & 0.098 & 0.092 & 0.104 & 0.087 & 0.077 & 0.096 \\
## 10 & Mississippi & 0.097 & 0.091 & 0.103 & 0.080 & 0.071 & 0.089 \\
## 2 & Arkansas & 0.093 & 0.087 & 0.099 & 0.080 & 0.072 & 0.089 \\
## 14 & Oklahoma & 0.081 & 0.078 & 0.084 & 0.072 & 0.066 & 0.077 \\
## 6 & Kansas & 0.066 & 0.063 & 0.069 & 0.059 & 0.053 & 0.064 \\
## 4 & Indiana & 0.065 & 0.060 & 0.069 & 0.060 & 0.052 & 0.067 \\
## 5 & Iowa & 0.058 & 0.055 & 0.062 & 0.052 & 0.046 & 0.057 \\
## 8 & Louisiana & 0.055 & 0.052 & 0.058 & 0.042 & 0.037 & 0.047 \\
## 15 & Tennessee & 0.049 & 0.046 & 0.053 & 0.038 & 0.032 & 0.043 \\
## 12 & Nebraska & 0.047 & 0.045 & 0.050 & 0.042 & 0.037 & 0.048 \\
## 3 & Illinois & 0.041 & 0.038 & 0.043 & 0.032 & 0.029 & 0.036 \\
## 11 & Missouri & 0.037 & 0.035 & 0.040 & 0.029 & 0.025 & 0.032 \\
## 16 & Wisconsin & 0.035 & 0.032 & 0.037 & 0.030 & 0.026 & 0.034 \\
## 13 & Ohio & 0.025 & 0.023 & 0.027 & 0.020 & 0.017 & 0.023 \\
## 9 & Minnesota & 0.023 & 0.021 & 0.024 & 0.019 & 0.016 & 0.022 \\
## 7 & Kentucky & 0.017 & 0.015 & 0.019 & 0.013 & 0.010 & 0.015 \\
## \hline
## \end{tabular}
## \end{table}
Adding state population statistics from www.census.gov from 2010 and 1980 to assess tornado exposure (and changes overtime)
tmp = download.file("http://www.census.gov/2010census/csv/pop_change.csv", "pop_change.csv",
mode = "wb")
USPop = read.csv("pop_change.csv", header = TRUE, skip = 2)
StPop = data.frame(USPop$STATE_OR_REGION, USPop$X1980_POPULATION, USPop$X2010_POPULATION)
StPop = rename(StPop, c(USPop.STATE_OR_REGION = "State", USPop.X2010_POPULATION = "Pop2010",
USPop.X1980_POPULATION = "Pop1980"))
Final1b = merge(Final1a, StPop, by = "State", all.y = FALSE, sort = FALSE, keep_order = 1)
Final1b$Exposure2010 = (Final1b$PA/100) * Final1b$Pop2010
Final1b$Exposure1980 = (Final1b$PA/100) * Final1b$Pop1980
Final1b$PercentIncrease = ((Final1b$Pop2010 - Final1b$Pop1980)/Final1b$Pop1980) *
100
Final1c = Final1b[order(Final1b$Exposure2010, decreasing = TRUE), ]
xtable(Final1c[, c(1, 10:12)], digits = 0) #Table 4
## % latex table generated in R 3.0.1 by xtable 1.7-1 package
## % Sun Aug 04 15:42:01 2013
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrrr}
## \hline
## & State & Exposure2010 & Exposure1980 & PercentIncrease \\
## \hline
## 11 & Illinois & 5207 & 4637 & 12 \\
## 1 & Alabama & 4662 & 3798 & 23 \\
## 6 & Indiana & 4196 & 3553 & 18 \\
## 9 & Tennessee & 3140 & 2272 & 38 \\
## 4 & Oklahoma & 3036 & 2449 & 24 \\
## 14 & Ohio & 2927 & 2740 & 7 \\
## 2 & Mississippi & 2878 & 2445 & 18 \\
## 3 & Arkansas & 2706 & 2122 & 28 \\
## 8 & Louisiana & 2490 & 2310 & 8 \\
## 12 & Missouri & 2243 & 1841 & 22 \\
## 13 & Wisconsin & 1981 & 1639 & 21 \\
## 5 & Kansas & 1886 & 1562 & 21 \\
## 7 & Iowa & 1780 & 1702 & 5 \\
## 15 & Minnesota & 1206 & 926 & 30 \\
## 10 & Nebraska & 867 & 745 & 16 \\
## 16 & Kentucky & 736 & 621 & 19 \\
## \hline
## \end{tabular}
## \end{table}
# write.table(Final1b, file='TorExp.txt')
Correlation Analysis with data from Ashley (2007)
NKE = c(1.473, 2.21, 1.926, 1.986, 1.392, 0.822, 0.76, 0.969, 0.442, 1.202,
1.489, 1.194, 0.69, 0.495, 0.83, 0.836) #normalized killer events in order (by state) of Finalb$PA
Final1b$NKE = NKE
cor.test(Final1b$PA, NKE)
##
## Pearson's product-moment correlation
##
## data: Final1b$PA and NKE
## t = 4.11, df = 14, p-value = 0.00106
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3849 0.9039
## sample estimates:
## cor
## 0.7395
cor.test(Final1b$PS, NKE)
##
## Pearson's product-moment correlation
##
## data: Final1b$PS and NKE
## t = 3.901, df = 14, p-value = 0.001597
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3520 0.8966
## sample estimates:
## cor
## 0.7217
cor.test(Final1b$Exposure2010, NKE)
##
## Pearson's product-moment correlation
##
## data: Final1b$Exposure2010 and NKE
## t = 0.9544, df = 14, p-value = 0.3561
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2832 0.6618
## sample estimates:
## cor
## 0.2472
cor.test(Final1b$Exposure1980, NKE)
##
## Pearson's product-moment correlation
##
## data: Final1b$Exposure1980 and NKE
## t = 0.8617, df = 14, p-value = 0.4034
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3052 0.6480
## sample estimates:
## cor
## 0.2244
Correlation between our corrected tornado rates and the normalized killer events (per sq. km X 1000) is 0.6714 (0.2634, 0.8757).
Based on the code provided by Nathan Yau (Flowing Data).
torExp = Final1b
torExp$Exposure1980 = round(torExp$Exposure1980)
torExp$Exposure2010 = round(torExp$Exposure2010)
x0 = c()
y0 = c()
x1 = c()
y1 = c()
startyear = 1980
stopyear = 2010
xoffset = 3.5
yoffset = 0
ystartprev = 0
ystopprev = 0
ythreshold = (max(torExp$Exposure1980) - min(torExp$Exposure1980)) * 0.0015
for (i in length(torExp$State):1) {
vals = torExp[i, ]
ystartdiff = (vals$Exposure1980 + yoffset) - ystartprev
if (abs(ystartdiff) < ythreshold) {
yoffset = yoffset + (ythreshold - ystartdiff)
}
# Calculate slope
slope = (vals$Exposure2010 - vals$Exposure1980)/(stopyear - startyear)
# Intercept
intercept = vals$Exposure1980 + yoffset
# Start and stop coordinates for lines
ystart = intercept
ystop = slope * (stopyear - startyear) + intercept
ystopdiff = ystop - ystopprev
if (abs(ystopdiff) < ythreshold) {
yoffset = yoffset + (ythreshold - ystopdiff)
intercept = vals$Exposure1980 + yoffset
ystart = intercept
ystop = slope * (stopyear - startyear) + intercept
}
# Draw the line for current state
x0 = c(x0, startyear)
y0 = c(y0, ystart)
x1 = c(x1, stopyear)
y1 = c(y1, ystop)
ystartprev = ystart
ystopprev = ystop
}
# Make adjustments for legibility
y0[13] = y0[13] + 25
y1[13] = y1[13] + 25
y0[15] = y0[15] - 25
y1[15] = y1[15] - 25
y0[8] = y0[8] - 5
y1[8] = y1[8] - 5
y0[9] = y0[9] + 25
y1[9] = y1[9] + 25
ymin = min(torExp$Exposure1980)
ymax = max(c(torExp$Exposure1980, torExp$Exposure2010)) + yoffset
par(family = "serif", mar = c(0, 0, 0, 0))
plot(0, 0, type = "n", main = "", xlab = "", ylab = "", xlim = c(1950, 2030),
ylim = c(ymin, ymax * 1.1), bty = "n", las = 1, axes = FALSE)
segments(x0, y0, x1, y1)
text(x0, y0, rev(torExp$Exposure1980), pos = 2, cex = 0.6)
text(x0 - xoffset, y0, rev(torExp$State), pos = 2, cex = 0.6)
text(x1, y1, rev(torExp$Exposure2010), pos = 4, cex = 0.6)
text(x1 + xoffset, y1, rev(torExp$State), pos = 4, cex = 0.6)
# Title
text(1950, ymax * 0.97, "Average Annual\nTornado Exposure\n[Number of People],\n1980 and 2010",
cex = 0.8, pos = 4)
# Year labels
text(startyear, ymax * 1.1, "1980", cex = 0.7, pos = 2, offset = 1)
text(stopyear, ymax * 1.1, "2010", cex = 0.7, pos = 4, offset = 0.5)
