County-Level Tornado Rates Adjusted for Population
James B. Elsner and Thomas H. Jagger
County-level analysis by individual states. The framework uses spatial polygons for organizing the data.
FP = 47
AB = "TN"Get tornado data
The model is built using data from the SPC tornado dataset. Data are stored as a spatial lines data frame. Transform the native coordinates to a Mercator projection.
if(!file.exists("torn")){
download.file(url = "http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
destfile = "tornado.zip")
unzip("tornado.zip")
}
library("rgdal")## Loading required package: sp## rgdal: version: 1.0-4, (SVN revision 548)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.2/Resources/library/rgdal/proj
## Linking to sp version: 1.1-1TornL = readOGR(dsn = "torn", layer = "torn",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "torn", layer: "torn"
## with 60114 features
## It has 22 fieldsTornL = spTransform(TornL, CRS("+init=epsg:3857"))Get County Boundaries
Boundaries are from https://www.census.gov/geo/maps-data/data/cbf/cbf_counties.html. The 20m = 1:20,000,000 low resolution boundaries. Here we use the 5m = 1:5,000,000.
if(!file.exists("cb_2014_us_county_5m.shp")){
download.file("http://www2.census.gov/geo/tiger/GENZ2014/shp/cb_2014_us_county_5m.zip",
"cb_2014_us_county_5m.zip", mode = "wb")
unzip("cb_2014_us_county_5m.zip")
}
US.sp = readOGR(dsn = ".", layer = "cb_2014_us_county_5m",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "cb_2014_us_county_5m"
## with 3233 features
## It has 9 fields## Warning in readOGR(dsn = ".", layer = "cb_2014_us_county_5m",
## stringsAsFactors = FALSE): Z-dimension discardedExtract the state using the state Federal Information Processing Standard (FIPS). Print out the number of counties and the area of state in square km. Check with Wolfram alpha.
library("rgeos")## rgeos version: 0.3-11, (SVN revision 479)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.1-0
## Polygon checking: TRUEST.sp = US.sp[US.sp$STATEFP == FP, ]
nrow(ST.sp) #number of counties## [1] 95county = paste(ST.sp$STATEFP, ST.sp$COUNTYFP, sep = "")
countiesun.sp = geometry(spChFIDs(ST.sp, county))
counties.sp = spTransform(countiesun.sp, CRS(proj4string(TornL)))
gArea(counties.sp) / 10^6 #total area of the state in square km## [1] 166482Organize county-level population data
Use the population estimates from the Census Bureau from 1970-2012 see http://www.nber.org/data/census-intercensal-county-population.html. Use the 2012 value for 2013 & 2014. Census U.S. Intercensal County Population Data, 1970-2014 from http://www.nber.org/data/census-intercensal-county-population.html.
library("reshape2")
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionL = "http://myweb.fsu.edu/jelsner/data/county_population.csv"
#L = "county_population.csv"
#pop = read.csv(L) %>%
# filter(county_fips < 1000 & county_fips > 0)
pop = read.csv(L) %>%
filter(state_fips %in% FP, county_fips != 0) %>%
mutate(pop2013 = pop2012,
pop2014 = pop2012,
pop2015 = pop2012,
county = as.character(fips)) %>%
select(county, pop1970:pop2012, pop2013, pop2014, pop2015)
Pop.df = melt(pop, id.vars = "county") %>%
mutate(yr = as.numeric(substring(variable, first = 4, last = 7))) %>%
rename(pop = value) %>%
mutate(lpop = log10(pop)) %>%
select(county, pop, yr, lpop)
Pop.df$ID = match(Pop.df$county, county) #generate spatial ID
Pop.df = Pop.df[order(Pop.df$ID), ] #order by spatial IDPC.df = Pop.df %>%
group_by(ID) %>%
summarize(Change = (pop[yr == max(yr)] - pop[yr == min(yr)])/pop[yr == max(yr)] * 100) %>%
data.frame()
row.names(PC.df) = county
spdf = SpatialPolygonsDataFrame(counties.sp, PC.df)
spdf$Name = ST.sp$NAME
spdf$area = rgeos::gArea(counties.sp, byid = TRUE)Latest population values by county are needed for correlation and choropleth. The log of the population density is the common log.
poplatest = subset(Pop.df, yr == max(yr))
spdf$pop = poplatest$pop
spdf$Lpop = log10(spdf$pop)
spdf$density = spdf$pop/spdf$area * 10^6
spdf$Ldensity = log10(spdf$density)Buffer tornado tracks to make tornado paths
Add a buffer to the tracks based on the width specified in the attribute table to create a spatial polygons data frame of tornado paths. Overlay paths onto the extent and return a vector indicating either NA (is not inside extent) or the county number. Remove duplicate paths. Start with 1970 for a 46-year climatology.
syr = 1970
eyr = 2015
Width = TornL$wid * .9144
sum(Width[Width == 0])## [1] 0TornP = gBuffer(TornL, byid = TRUE, width = Width/2, capStyle = "FLAT")
tc = over(TornP, counties.sp)
TornP2 = subset(TornP, !is.na(tc))
TornP2 = subset(TornP2, yr >= syr)
df = data.frame(SLAT = TornP2$slat,
SLON = TornP2$slon,
ELAT = TornP2$elat,
ELON = TornP2$elon,
DATE = TornP2$date)
dup = duplicated(df)
sum(dup)## [1] 4sum(dup)/dim(TornP2@data)[1] * 100## [1] 0.3960396TornP3 = subset(TornP2, !dup)
dim(TornP3@data)## [1] 1006 22dim(TornP3@data)[1] / (gArea(counties.sp) / 10^6) * 10000 / (2015 - 1970 + 1) # 1.31: Avg annual state rate per 100 sq. km## [1] 1.313629All EF0+ tornadoes since 1970. Overlay tornado paths on counties. Each element of the list is a county subset of the tornado track file. The data frame nTor.df contains the per-county tornado count, area, and rate.
The PathsAll.df data frame contains the tornado path attributes listed by county. Information is repeated for each county in cases where the path intersects more than one county.
TornP3 = subset(TornP3, mag >= 0) #EF0+
ct = over(counties.sp, TornP3, returnList = TRUE)
area = gArea(counties.sp, byid = TRUE)
names(ct) = county
nT = sapply(ct, function(x) nrow(x))
nTd = sapply(ct, function(x) length(unique(x$date)))
Nyears = diff(range(TornP3$yr)) + 1
nTor.df = data.frame(state = AB,
county = county,
Nyears = Nyears,
nT = nT,
area = area/1000000,
nTd = nTd,
rate = nT/(area / 10^6) * 10000 / Nyears,
stringsAsFactors = FALSE, ID = 1:length(county)
)
spdf$state = AB
spdf$county = county
spdf$Nyears = nTor.df$Nyears
spdf$nT = nTor.df$nT
spdf$nTd = nTor.df$nTd
spdf$area = nTor.df$area
spdf$rate = nTor.df$rate
PathsAll.df = plyr::ldply(ct, data.frame)
colnames(PathsAll.df)[1] = "county"
sort(spdf$area)## [1] 471.2718 496.6895 697.3847 702.3864 740.9926 776.2061 829.2528
## [8] 845.3091 902.6165 977.0308 986.1245 1033.8463 1046.0849 1050.0793
## [15] 1053.1146 1086.0437 1108.3788 1116.1085 1205.4636 1220.4810 1225.5208
## [22] 1235.5835 1246.6954 1278.6924 1283.8842 1295.9205 1296.7618 1318.3275
## [29] 1319.5332 1349.8143 1367.7236 1383.7417 1411.9075 1475.6909 1497.6494
## [36] 1557.7475 1598.5121 1598.6117 1663.2101 1696.5629 1702.9063 1712.6330
## [43] 1716.8347 1725.5238 1733.2496 1740.4123 1755.7491 1773.5320 1854.4404
## [50] 1864.3619 1918.0113 1950.8371 1980.4611 1992.2595 1992.6467 1996.2897
## [57] 1998.6578 2001.0987 2069.2755 2083.5847 2084.9052 2087.5182 2087.7313
## [64] 2094.9583 2135.9904 2174.0580 2181.1142 2191.8166 2194.0038 2207.3125
## [71] 2220.6882 2220.8744 2232.7677 2235.2461 2241.8724 2303.9633 2318.6417
## [78] 2322.2051 2330.6950 2356.2049 2377.2276 2378.8104 2379.4917 2390.0942
## [85] 2390.1021 2419.9674 2420.8936 2457.4910 2493.8274 2550.0284 2607.1324
## [92] 2715.1687 2750.6281 2871.9707 3046.9015spdf$Name[order(spdf$area)]## [1] "Trousdale" "Moore" "Hamblen" "Pickett" "Unicoi"
## [6] "Lake" "Houston" "Meigs" "Hancock" "Loudon"
## [11] "Union" "Sequatchie" "Clay" "Cannon" "Crockett"
## [16] "Van Buren" "Lewis" "Chester" "Grainger" "Johnson"
## [21] "Cheatham" "Macon" "Jefferson" "Jackson" "Bradley"
## [26] "DeKalb" "Smith" "Washington" "Rhea" "Decatur"
## [31] "Anderson" "Carter" "Grundy" "Marshall" "White"
## [36] "Roane" "Putnam" "Bledsoe" "Perry" "McMinn"
## [41] "Coffee" "Warren" "Polk" "Sullivan" "Benton"
## [46] "Overton" "Cocke" "Claiborne" "Tipton" "Bedford"
## [51] "Robertson" "Dickson" "Stewart" "Marion" "Fentress"
## [56] "Campbell" "Hawkins" "Lauderdale" "Henderson" "Morgan"
## [61] "Knox" "Dyer" "Davidson" "Haywood" "Scott"
## [66] "Sumner" "Montgomery" "McNairy" "Madison" "Humphreys"
## [71] "Obion" "Lincoln" "Blount" "Franklin" "Hamilton"
## [76] "Williamson" "Hardin" "Weakley" "Wilson" "Sevier"
## [81] "Carroll" "Giles" "Henry" "Lawrence" "Gibson"
## [86] "Hickman" "Maury" "Rutherford" "Greene" "Monroe"
## [91] "Hardeman" "Cumberland" "Fayette" "Wayne" "Shelby"sort(spdf$nT)## [1] 1 1 2 2 2 2 2 3 3 3 4 4 4 4 4 4 5 5 5 5 6 6 6
## [24] 7 7 7 7 7 7 8 8 8 8 8 8 9 9 9 9 9 9 9 10 10 11 11
## [47] 12 12 12 12 12 12 13 13 13 14 14 14 15 15 15 15 15 15 15 15 16 16 16
## [70] 17 17 17 18 18 18 19 19 19 20 20 21 22 23 23 24 26 26 28 29 30 30 30
## [93] 32 36 38mean(spdf$nT); median(spdf$nT); sd(spdf$nT)## [1] 12.83158## [1] 12## [1] 8.409281spdf$Name[order(spdf$nT)]## [1] "Hamblen" "Campbell" "Hancock" "Clay" "Hawkins"
## [6] "Johnson" "Unicoi" "Grainger" "Sevier" "Union"
## [11] "Morgan" "Carter" "Meigs" "Rhea" "Sullivan"
## [16] "Van Buren" "Anderson" "Grundy" "Cocke" "Roane"
## [21] "Trousdale" "Houston" "Jackson" "Pickett" "Perry"
## [26] "Loudon" "Decatur" "Lake" "Bledsoe" "Chester"
## [31] "Obion" "Jefferson" "Lewis" "Sequatchie" "Washington"
## [36] "Cannon" "Macon" "Crockett" "Scott" "Smith"
## [41] "Greene" "Claiborne" "Hardin" "DeKalb" "Polk"
## [46] "Cheatham" "White" "Haywood" "Moore" "Warren"
## [51] "Overton" "Fentress" "Hickman" "Stewart" "Knox"
## [56] "Dickson" "Marshall" "Maury" "Henry" "Giles"
## [61] "McNairy" "Lauderdale" "Wayne" "Blount" "Bedford"
## [66] "Putnam" "Fayette" "Henderson" "Carroll" "Monroe"
## [71] "Marion" "Weakley" "Humphreys" "Tipton" "Cumberland"
## [76] "Hardeman" "Robertson" "McMinn" "Benton" "Williamson"
## [81] "Bradley" "Coffee" "Gibson" "Dyer" "Franklin"
## [86] "Lawrence" "Montgomery" "Madison" "Hamilton" "Davidson"
## [91] "Wilson" "Lincoln" "Rutherford" "Sumner" "Shelby"cor.test(spdf$nT, spdf$area, conf.level = .9)##
## Pearson's product-moment correlation
##
## data: spdf$nT and spdf$area
## t = 7.5399, df = 93, p-value = 3.063e-11
## alternative hypothesis: true correlation is not equal to 0
## 90 percent confidence interval:
## 0.4982281 0.7113552
## sample estimates:
## cor
## 0.6159385cor.test(spdf$nT, spdf$pop, conf.level = .9)##
## Pearson's product-moment correlation
##
## data: spdf$nT and spdf$pop
## t = 5.9727, df = 93, p-value = 4.238e-08
## alternative hypothesis: true correlation is not equal to 0
## 90 percent confidence interval:
## 0.3917367 0.6392034
## sample estimates:
## cor
## 0.5265342cor.test(spdf$area, spdf$pop, conf.level = .9)##
## Pearson's product-moment correlation
##
## data: spdf$area and spdf$pop
## t = 3.5683, df = 93, p-value = 0.0005704
## alternative hypothesis: true correlation is not equal to 0
## 90 percent confidence interval:
## 0.1882945 0.4880855
## sample estimates:
## cor
## 0.3470239Area is now given in square km. Rate is tornadoes per sq m per 10,000 years. Choropleth map of tornado counts & tornado days.
library("wesanderson")
range(spdf$nT)## [1] 1 38rng = seq(0, 40, 5)
crq = wes_palette(8, name = "Zissou", type = "continuous")
spplot(spdf, "nT", col = "white", at = rng,
col.regions = crq,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = "Number of EF0+ Tornadoes (1970-2015)")
library("RColorBrewer")
range(spdf$nTd)## [1] 1 32rng = seq(0, 40, 5)
crq = brewer.pal(8, "GnBu")
spplot(spdf, "nTd", col = "white", at = rng,
col.regions = crq,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = "Number of EF0+ Tornado Days (1970-2014)")
range(spdf$rate)## [1] 0.1088977 5.2521655rng = seq(0, 6, 1)
crq = wes_palette(6, name = "Zissou", type = "continuous")
spplot(spdf, "rate", col = "white", at = rng,
col.regions = crq,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = "Annual Tornado Rate Per 100 sq km (1970-2015)") 
mean(spdf$rate) # 1.57## [1] 1.572299Annual tornado counties
The data frame nTor.year is the number of tornado counties by year.
nTor.year = PathsAll.df %>%
group_by(yr) %>%
summarize(nT = n()) %>%
data.frame()
nTor.year$yr2 = nTor.year$yr
sum(nTor.year$nT)## [1] 1219The data frame nTor.day0 contains the tornado county and tornado day county by year and county.
nTor.day0 = PathsAll.df %>%
group_by(county, yr) %>%
summarize(nT = length(county),
nTd = length(unique(date)))
nTor.day0$county = as.factor(nTor.day0$county)
nTor.day0$yr = as.factor(nTor.day0$yr)
nTor.day0 = left_join(expand.grid(county = levels(as.factor(PathsAll.df$county)),
yr = levels(as.factor(PathsAll.df$yr))),
nTor.day0)## Joining by: c("county", "yr")nTor.day0[is.na(nTor.day0)] = 0
nTor.day0$county = as.character(nTor.day0$county)
nTor.day0$yr = as.integer(as.character(nTor.day0$yr))
nTor.day = base::merge(nTor.df[c("state", "county", "area", "ID")],
nTor.day0) %>%
arrange(ID, yr)Merge tornado and population data.
popNtor = base::merge(nTor.day, Pop.df, all = FALSE) %>%
mutate(density = pop/area,
Ldensity = log2(density),
yr2 = yr) %>% #needed for the INLA models
arrange(ID, yr)Trends in statewide tornado counts
Regress tornado-county counts on year. First with Year as a fixed covariate then with Year as a second order random walk term of the form
\[x[i] - 2 x[i-1] + x[i-2] \sim \mathcal{N}(0,\tau^{-1})\]
Some global controls for INLA. Use as needed.
control = list(
predictor = list(compute = TRUE),
inla = list(strategy = "laplace",
fast = FALSE,
stencil = 7,
npoints = 198,
int.strategy = "grid",
dz = .5),
results = list(return.marginals.random = TRUE),
compute = list(config = TRUE, mlik = TRUE, cpo = TRUE, dic = TRUE, po = TRUE),
family = list(variant = 1, hyper = list(theta = list(prior = "loggamma", param = c(1, 1)))))Linear Trend (%/yr)
Torn.df = as.data.frame(TornP3) %>%
group_by(yr) %>%
summarize(nT = n())
Torn.df$yr2 = Torn.df$yr
#source("http://www.math.ntnu.no/inla/givemeINLA.R")
library("INLA")## Loading required package: Matrix## Loading required package: splinesformula = nT ~ yr
modelt0 = inla(formula = formula,
family = "nbinomial",
quantiles = c(.05, .5, .95),
data = Torn.df,
control.compute = control$compute,
control.predictor = control$predictor)
TM = modelt0$summary.fixed$mean[2]
TL = modelt0$summary.fixed$'0.05quant'[2]
TH = modelt0$summary.fixed$'0.95quant'[2]
(exp(TM) - 1) * 100## [1] 2.47951(exp(TL) - 1) * 100## [1] 1.32612(exp(TH) - 1) * 100## [1] 3.644021Nonlinear trend
formula = nT ~ yr + f(yr2, model = "rw2")
modelt1 = inla(formula = formula,
family = "nbinomial",
quantiles = c(.05, .5, .95),
data = Torn.df,
control.compute = control$compute,
control.predictor = control$predictor)
out = data.frame(Year = nTor.year$yr,
QL = modelt1$summary.fitted.values$"0.05quant",
QM = modelt1$summary.fitted.values$"0.5quant",
QH = modelt1$summary.fitted.values$"0.95quant",
nT = nTor.year$nT,
ST = "Tennessee")
library("ggplot2")
ggplot(out, aes(x = Year, y = QM)) +
geom_line(color = 'blue') +
geom_ribbon(aes(x = Year, ymax = QH, ymin = QL), color = "gray", alpha = .25) +
geom_point(aes(x = Year, y = nT)) +
ylab("Number of EF0+ Tornadoes") +
scale_y_continuous(limits = c(0, 150))
stateRate = modelt1$summary.fitted.values$mean[Nyears] / (gArea(counties.sp) / 10^6) * 10000Models for county-level tornado counts
Spatial neighborhood definition as an inla graph
library("spdep")
nb = poly2nb(spdf)
nb2INLA("tornb.inla", nb)
tornb.inla = inla.read.graph("tornb.inla")The area multiplied by Nyears indicates the square-meter years exposed to tornadoes (E). Scale the exposure to have a mean close to unity. The model has a spatial id. E is scaled to have a mean of one. constr = constant rate. The DIC measures how well the model fits the data; the larger this is, the worse the fit. Year is centered on 1991.
Model with spatial term. Determine the model that results in the lowest DIC.
meanArea = mean(popNtor$area)
formula = nT ~ f(ID, model = "besag", graph = tornb.inla, constr = TRUE) +
Ldensity + I(yr - 1991) +
f(yr2, model = "iid")
model = inla(formula = formula, family = "nbinomial", E = area/meanArea,
data = popNtor,
quantiles = c(.05, .5, .95),
control.compute = control$compute,
control.predictor = control$predictor,
control.results = control$results,
control.family = control$family
)
summary(model)##
## Call:
## c("inla(formula = formula, family = \"nbinomial\", data = popNtor, ", " quantiles = c(0.05, 0.5, 0.95), E = area/meanArea, control.compute = control$compute, ", " control.predictor = control$predictor, control.family = control$family, ", " control.results = control$results)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.9213 15.1902 1.2137 17.3251
##
## Fixed effects:
## mean sd 0.05quant 0.5quant 0.95quant mode kld
## (Intercept) -2.4576 0.1825 -2.7605 -2.4562 -2.1599 -2.4532 0
## Ldensity 0.1821 0.0349 0.1250 0.1819 0.2399 0.1815 0
## I(yr - 1991) 0.0242 0.0075 0.0119 0.0242 0.0366 0.0241 0
##
## Random effects:
## Name Model
## ID Besags ICAR model
## yr2 IID model
##
## Model hyperparameters:
## mean sd
## size for the nbinomial observations (overdispersion) 0.9279 0.1125
## Precision for ID 7.2531 2.7643
## Precision for yr2 2.6845 0.6878
## 0.05quant 0.5quant
## size for the nbinomial observations (overdispersion) 0.7582 0.9192
## Precision for ID 3.7596 6.7469
## Precision for yr2 1.7007 2.6087
## 0.95quant mode
## size for the nbinomial observations (overdispersion) 1.127 0.9005
## Precision for ID 12.472 5.8509
## Precision for yr2 3.932 2.4671
##
## Expected number of effective parameters(std dev): 63.59(5.285)
## Number of equivalent replicates : 68.73
##
## Deviance Information Criterion (DIC) ...: 5162.63
## Effective number of parameters .........: 64.67
##
## Watanabe-Akaike information criterion (WAIC) ...: 5163.53
## Effective number of parameters .................: 62.63
##
## Marginal log-Likelihood: -2711.66
## CPO and PIT are computed
##
## Posterior marginals for linear predictor and fitted values computedModel diagnostics. The predictive quality of the model is assessed by the cross-validated log score. A smaller value of the score indicates better prediction quality (Gneiting and Raftery 2007).
modPIT = model$cpo$pit - runif(length(model$cpo$cpo)) * model$cpo$cpo
goftest::ad.test(modPIT)##
## Anderson-Darling test of goodness-of-fit
## Null hypothesis: uniform distribution
##
## data: modPIT
## An = 0.84585, p-value = 0.4491-mean(log(model$cpo$cpo))## [1] 0.5908513cor.test(popNtor$nT, model$summary.fitted.values$mean, conf.level = .9)##
## Pearson's product-moment correlation
##
## data: popNtor$nT and model$summary.fitted.values$mean
## t = 27.493, df = 4368, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 90 percent confidence interval:
## 0.3626662 0.4050988
## sample estimates:
## cor
## 0.3840853Brier score
brier.score <- function(x, m){
with(m, {mean(x^2) - 2 * mean(x * mean) + mean(mean^2 + sd^2)})
}
brier.score(popNtor[["nT"]], model[["summary.fitted.values"]])## [1] 0.4610103ggplot() +
geom_histogram(aes(x = modPIT), color = "white",
fill = "grey", binwidth = .05) +
xlab("Probability") +
ylab("Frequency") +
ggtitle("Distribution of modified PIT") #looks good!
The spatial random term by county is given in the data frame model\(summary.random\)ID. Normalize the stateRate by the number of counties. Multiply the exponentiated spatial random term by the statewide mean rate per county.
df = model$summary.random$ID
names(df) = c("ID", "mean", "sd", "QL", "QM", "QH", "mode", "kld")
df = df %>%
mutate(QL = exp(QL) - 1,
QH = exp(QH) - 1,
Sig = sign(QL) == sign(QH),
sd = sd,
ctyRate = stateRate * exp(mean),
ctyPerState = (exp(mean) - 1) * 100,
ID = county)
sum(df$Sig)## [1] 50spdfR = spdf
spdfR@data = df
range(spdfR$ctyPerState)## [1] -65.27702 77.21458spdf$Name[order(spdfR$ctyPerState)]## [1] "Johnson" "Sullivan" "Carter" "Hawkins" "Unicoi"
## [6] "Washington" "Greene" "Hamblen" "Cocke" "Sevier"
## [11] "Hancock" "Grainger" "Jefferson" "Campbell" "Union"
## [16] "Claiborne" "Knox" "Anderson" "Blount" "Morgan"
## [21] "Roane" "Scott" "Loudon" "Rhea" "Monroe"
## [26] "Meigs" "Cumberland" "Bledsoe" "Fentress" "McMinn"
## [31] "Hamilton" "Clay" "Polk" "Pickett" "Van Buren"
## [36] "Overton" "Jackson" "Sequatchie" "Putnam" "Bradley"
## [41] "White" "Warren" "Hardin" "Grundy" "Obion"
## [46] "Maury" "Shelby" "Smith" "Hickman" "DeKalb"
## [51] "Macon" "Wayne" "Marion" "Williamson" "Decatur"
## [56] "Fayette" "Perry" "Dickson" "Davidson" "McNairy"
## [61] "Lewis" "Cheatham" "Chester" "Weakley" "Henderson"
## [66] "Henry" "Carroll" "Cannon" "Haywood" "Giles"
## [71] "Rutherford" "Tipton" "Hardeman" "Robertson" "Houston"
## [76] "Montgomery" "Marshall" "Trousdale" "Humphreys" "Stewart"
## [81] "Gibson" "Bedford" "Wilson" "Lauderdale" "Coffee"
## [86] "Lake" "Crockett" "Lawrence" "Madison" "Benton"
## [91] "Dyer" "Franklin" "Sumner" "Lincoln" "Moore"rng = seq(-80, 80, 20)
library("RColorBrewer")
crq = rev(brewer.pal(8, "RdBu"))
spplot(spdfR, "ctyPerState", col = "white", at = rng,
col.regions = crq,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = "Tornado (EF0+) Occurrence Rate [% Difference from Statewide Avg Rate]")
range(spdfR$ctyRate)## [1] 0.8581178 4.3795493rng = seq(0, 5, 1)
crq = brewer.pal(5, "Oranges")
spplot(spdfR, "ctyRate", col = "white", at = rng,
col.regions = crq,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = "Expected Annual Tornado (EF0+) Rate Per 100 sq. km")
range(spdfR$sd)## [1] 0.138139 0.311528rng = seq(.1, .35, .05)
crq = brewer.pal(5, "Greens")
spplot(spdfR, "sd", col = "white", at = rng,
col.regions = crq,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = "Standard Error")
Write out the shapefiles.
writeOGR(spdfR, dsn = AB, layer = "Rates", overwrite_layer = TRUE,
driver = "ESRI Shapefile")## Warning in writeOGR(spdfR, dsn = AB, layer = "Rates", overwrite_layer =
## TRUE, : Field names abbreviated for ESRI Shapefile driver