A dasymetric method to spatially apportion tornado casualty counts
Tyler Fricker, James B. Elsner, Victor Mesev, and Thomas H. Jagger
Get tornado data and set domain. Subset storms by domain
#download.file("http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
# "tornado.zip", mode = "wb")
#unzip("tornado.zip")
library("rgdal")## Loading required package: sp## rgdal: version: 1.2-7, (SVN revision 660)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/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.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-4TornL = readOGR(dsn = "torn", layer = "torn",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "torn", layer: "torn"
## with 61092 features
## It has 22 fieldslibrary("raster")
xmn = -104
xmx = -82
ymn = 31
ymx = 41
res = .5
#res = 0.25
#res = 0.75
syr = 1955
eyr = 2016
r = raster(xmn = xmn, xmx = xmx,
ymn = ymn, ymx = ymx,
resolution = res)
sp = as(r, 'SpatialPolygons')
spT = spTransform(sp, CRS(proj4string(TornL)))Percent of all tornadoes in the SPC dataset with touchdown locations within the domain.
df = TornL@data
pAll = dim(subset(df, yr >= syr & slat >= ymn & slat <= ymx &
slon >= xmn & slon <= xmx))[1]/dim(subset(df, yr >= syr))[1]
pInt = dim(subset(df, yr >= syr & mag >= 3 & slat >= ymn & slat <= ymx &
slon >= xmn & slon <= xmx))[1]/dim(subset(df, yr >= syr & mag >= 3))[1]
print(pAll); print(pInt)## [1] 0.5233423## [1] 0.6677215The study domain is chosen large enough to capture the most prolific tornado area of the country. As chosen it contains the majority (52%) of all tornadoes in the United States and two-thirds (67%) of all intense (EF3+) tornadoes. A grid of equal lat-lon cells is overlaid at a spatial resolution of 0.5 degrees latitude/longitude resulting in cells north-south and 44 cells east-west for a total of 880 cells.
Create tornado paths by buffering the tracks
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 domain extent and return a vector indicating either NA (tornado is not inside extent) or the cell number. Subset the tornadoes by the raster extent (domain). Remove duplicate paths from the data frame.
library("rgeos")## rgeos version: 0.3-23, (SVN revision 546)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-4
## Polygon checking: TRUEWidth = TornL$wid * .9144
sum(Width[Width == 0])## [1] 0TornP = gBuffer(TornL, byid = TRUE, width = Width/2, capStyle = "FLAT")
tc = over(TornP, spT)
TornP2 = subset(TornP, !is.na(tc) & 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] 355sum(dup)/dim(TornP2@data)[1] * 100## [1] 1.147976TornP3 = subset(TornP2, !dup)
TornP3 = TornP3[!TornP3$mag == -9,]
TornP3$cas = TornP3$fat + TornP3$inj
dim(TornP3@data)## [1] 30546 23l = gLength(TornP3[TornP3$fat > 0, ], byid = TRUE)
p = (4 * l * 111000 * res - l^2)/(pi * (111000 * res)^2)
sum(p > .5)/dim(TornP3[TornP3$fat > 0, ])[1] * 100 ## [1] 61.38728#percentage of all tornadoes with fatalities having a better than even chance of intersecting more than one grid cell
l = gLength(TornP3[TornP3$inj > 0, ], byid = TRUE)
p = (4 * l * 111000 * res - l^2)/(pi * (111000 * res)^2)
sum(p > .5)/dim(TornP3[TornP3$inj > 0, ])[1] * 100 ## [1] 43.10388#percentage of all tornadoes with injuries having a better than even chance of intersecting more than one grid cell
l = gLength(TornP3[TornP3$cas > 0, ], byid = TRUE)
p = (4 * l * 111000 * res - l^2)/(pi * (111000 * res)^2)
sum(p > .5)/dim(TornP3[TornP3$cas > 0, ])[1] * 100 ## [1] 42.97561#percentage of all tornadoes with casualties having a better than even chance of intersecting more than one grid cellNumber of killer tornadoes and their associated fatalities as well as the number of injury-producing tornadoes and their associated injuries.
sum(TornP3$fat > 0); sum(TornP3$fat)## [1] 865## [1] 3553sum(TornP3$inj > 0); sum(TornP3$inj)## [1] 3995## [1] 56424Distribution of tornadoes by casualties and highest damage rating.
library("xtable")
tA = table(TornP3$fat, TornP3$mag)[-1, ]
tA = rbind(tA[1:5,], colSums(tA[6:dim(tA)[1],]))
rownames(tA)[6] = "6+"
xtable(tA, digits = 0)## % latex table generated in R 3.3.3 by xtable 1.8-2 package
## % Thu Oct 5 15:46:44 2017
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrr}
## \hline
## & 0 & 1 & 2 & 3 & 4 & 5 \\
## \hline
## 1 & 5 & 69 & 147 & 155 & 38 & 3 \\
## 2 & 1 & 15 & 37 & 67 & 34 & 1 \\
## 3 & 1 & 3 & 11 & 39 & 29 & 3 \\
## 4 & 0 & 1 & 4 & 18 & 20 & 1 \\
## 5 & 0 & 1 & 1 & 11 & 10 & 1 \\
## 6+ & 0 & 1 & 3 & 28 & 81 & 26 \\
## \hline
## \end{tabular}
## \end{table}tB = table(TornP3$inj, TornP3$mag)[-1, ]
tB = rbind(tB[1:5,], colSums(tB[6:dim(tB)[1],]))
rownames(tB)[6] = "6+"
xtable(tB, digits = 0)## % latex table generated in R 3.3.3 by xtable 1.8-2 package
## % Thu Oct 5 15:46:44 2017
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrrrr}
## \hline
## & 0 & 1 & 2 & 3 & 4 & 5 \\
## \hline
## 1 & 80 & 468 & 440 & 143 & 12 & 0 \\
## 2 & 33 & 259 & 279 & 106 & 15 & 0 \\
## 3 & 14 & 114 & 181 & 75 & 9 & 0 \\
## 4 & 7 & 70 & 129 & 57 & 4 & 0 \\
## 5 & 3 & 48 & 95 & 46 & 11 & 0 \\
## 6+ & 8 & 121 & 403 & 484 & 246 & 35 \\
## \hline
## \end{tabular}
## \end{table}Table 1
There are 30546 tornado paths that intersect the domain. Table 1 shows the distributon of casualties by EF rating.
Map showing the paths of all killer tornadoes within the domain
Use ggmap. This requires the spatial polygons data frame be converted to the geographic projection of the raster.
TornP3T = spTransform(TornP3, CRS(projection(r)))
library("ggmap")## Loading required package: ggplot2library('ggplot2')
TornP3TK = TornP3T[TornP3T$fat > 0,]
#writeOGR(TornP3TK, dsn = "Killer_Torn", layer = "Killer_Torn", overwrite_layer = TRUE,
# driver = "ESRI Shapefile")
BBox = gConvexHull(spT)
df = data.frame(id = getSpPPolygonsIDSlots(BBox))## Warning: use *apply and slot directlyrow.names(df) <- getSpPPolygonsIDSlots(BBox)## Warning: use *apply and slot directlyspdf = SpatialPolygonsDataFrame(BBox, data = df)
#writeOGR(spdf, dsn = "Boundary", layer = "Boundary", overwrite_layer = TRUE,
# driver = "ESRI Shapefile")
mapdata = fortify(TornP3T[TornP3T$fat > 0, ])## Regions defined for each PolygonsMap = get_map(location = c(-((xmx - xmn)/2 - xmx) , (ymx - ymn)/2 + ymn),
source = "stamen",
zoom = 5,
color = "bw",
maptype = "toner")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=36,-93&zoom=5&size=640x640&scale=2&maptype=terrain&sensor=false## Map from URL : http://tile.stamen.com/toner/5/6/11.png## Map from URL : http://tile.stamen.com/toner/5/7/11.png## Map from URL : http://tile.stamen.com/toner/5/8/11.png## Map from URL : http://tile.stamen.com/toner/5/6/12.png## Map from URL : http://tile.stamen.com/toner/5/7/12.png## Map from URL : http://tile.stamen.com/toner/5/8/12.png## Map from URL : http://tile.stamen.com/toner/5/6/13.png## Map from URL : http://tile.stamen.com/toner/5/7/13.png## Map from URL : http://tile.stamen.com/toner/5/8/13.pngggmap(Map, dev = "extent") +
geom_polygon(aes(x = long, y = lat, group = group),
color = "black", fill = "black",
data = mapdata) +
geom_segment(aes(x = xmn, xend = xmx, y = ymn, yend = ymn),
color = "red") +
geom_segment(aes(x = xmn, xend = xmx, y = ymx, yend = ymx),
color = "red") +
geom_segment(aes(x = xmn, xend = xmn, y = ymn, yend = ymx),
color = "red") +
geom_segment(aes(x = xmx, xend = xmx, y = ymn, yend = ymx),
color = "red") +
labs(x = NULL, y = NULL) +
# labs(x = expression(paste("Longitude (", degree, "W)")),
# y = expression(paste("Latitude (", degree, "N)"))) +
theme(legend.position = "none")
Figure 1: Study area
Exploratory analysis on tornado casualties at the level of individual tornadoes. Distribution (monthly, annual), trend, etc.
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:rgeos':
##
## intersect, setdiff, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondf = TornP3@data
dfA = df %>%
group_by(yr) %>%
summarize(nD = sum(fat),
nI = sum(inj),
nkT = sum(fat > 0),
niT = sum(inj > 0))
A = ggplot(dfA, aes(x = yr, y = nD)) +
geom_line() +
geom_smooth() +
scale_x_continuous(breaks = seq(1955, 2015, 10)) +
xlab("Year") + ylab("Number of Tornado Deaths")
B = ggplot(dfA, aes(x = yr, y = nkT)) +
geom_line() +
geom_smooth() +
scale_x_continuous(breaks = seq(1955, 2015, 10)) +
xlab("Year") + ylab("Number of Killer Tornadoes")
C = df %>%
group_by(mo) %>%
summarize(nD = sum(fat)) %>%
mutate(MonthF = factor(month.abb[mo], levels = month.abb)) %>%
ggplot(., aes(x = MonthF, y = nD)) +
geom_bar(stat = "identity") +
scale_y_continuous(limits = c(0, 1500)) +
xlab("") + ylab("Number of Tornado Deaths")
library("VGAM")## Loading required package: stats4## Loading required package: splinesdfD = data.frame(table(df$fat[df$fat > 0]))
dfD$Size = as.integer(dfD$Var1)
#dfD$Size = as.integer(rownames(dfD))
D = ggplot(dfD, aes(x = Size, y = Freq)) +
scale_x_log10() + scale_y_log10() +
geom_point() +
xlab(expression(paste("Number of Deaths"))) +
ylab("Number of Tornadoes")
source("http://peterhaschke.com/Code/multiplot.R")
mat = matrix(c(1, 2, 3, 4), nrow = 2, byrow = TRUE)
A = A + ggtitle("A") +
theme(plot.title=element_text(hjust=0))
B = B + ggtitle("B") +
theme(plot.title=element_text(hjust=0))
C = C + ggtitle("C") +
theme(plot.title=element_text(hjust=0))
D = D + ggtitle("D") +
theme(plot.title=element_text(hjust=0))
multiplot(A, B, C, D, layout = mat)## Loading required package: grid## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
# Export at 10 x 9Figure 2: Tornado deaths.
Repeat for injuries.
A = ggplot(dfA, aes(x = yr, y = nI)) +
geom_line() +
geom_smooth() +
scale_x_continuous(breaks = seq(1955, 2015, 10)) +
xlab("Year") + ylab("Number of Tornado Injuries")
B = ggplot(dfA, aes(x = yr, y = niT)) +
geom_line() +
geom_smooth() +
scale_x_continuous(breaks = seq(1955, 2015, 10)) +
xlab("Year") + ylab("Number of Tornadoes with Injuries")
C = df %>%
group_by(mo) %>%
summarize(nI = sum(inj)) %>%
mutate(MonthF = factor(month.abb[mo], levels = month.abb)) %>%
ggplot(., aes(x = MonthF, y = nI)) +
geom_bar(stat = "identity") +
scale_y_continuous() +
xlab("") + ylab("Number of Tornado Injuries")
dfD = data.frame(table(df$inj[df$inj > 0]))
dfD$Size = as.integer(dfD$Var1)
D = ggplot(dfD, aes(x = Size, y = Freq)) +
scale_x_log10() + scale_y_log10() +
geom_point() +
xlab(expression(paste("Number of Injuries"))) +
ylab("Number of Tornadoes")
source("http://peterhaschke.com/Code/multiplot.R")
mat = matrix(c(1, 2, 3, 4), nrow = 2, byrow = TRUE)
A = A + ggtitle("A") +
theme(plot.title=element_text(hjust=0))
B = B + ggtitle("B") +
theme(plot.title=element_text(hjust=0))
C = C + ggtitle("C") +
theme(plot.title=element_text(hjust=0))
D = D + ggtitle("D") +
theme(plot.title=element_text(hjust=0))
multiplot(A, B, C, D, layout = mat)Proportional allocation of casualties
Determine the number of tornado paths intersecting each cell.
ct = over(spT, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
nt = r
values(nt) = nT
dim(TornP3)[1]## [1] 30546sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)## [1] 37007## [1] 42.05341## [1] 0## [1] 150## [1] 20.13534## [1] 9.640877The average number of tornadoes per cell is 42.1, with a minimum of 0 and a maximum of 150 tornadoes. There are 211 cells that have between 35 and 45 tornado paths. The standard deviation of the tornadoes is 20.1 with a variance to mean ratio (VMR) of 9.64.
Map the number of tornadoes. First get county & state boundaries and convert to the geographic coordinates of the raster.
library("mapproj")## Loading required package: mapslibrary("maptools")## Checking rgeos availability: TRUE##
## Attaching package: 'maptools'## The following object is masked from 'package:xtable':
##
## labelext = as.vector(extent(r))
bndryC = map("county", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(bndryC$names, ":"), function(x) x[1])
bndryCP = map2SpatialPolygons(bndryC, IDs = IDs,
proj4string = CRS(projection(r)))
bndryS = map("state", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(bndryS$names, ":"), function(x) x[1])
bndrySP = map2SpatialPolygons(bndryS, IDs = IDs,
proj4string = CRS(projection(r)))library("rasterVis")
library("wesanderson")
range(values(nt))
rng = seq(0, 160, 20)
cr = wes_palette(8, name = "Zissou", type = "continuous")
levelplot(nt, margin = FALSE,
sub = expression(" Number of Reported Tornadoes"),
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(space = 'bottom'),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1))Casualties in the SPC are given per tornado but where they occurred within the path (or surrounding) is unknown. The likelihood of a casualty obviously increases with the number of people affected so we assign to each cell under the tornado path a fraction of the casualty magnitude based on area of the path within the cell and the proportion of people in the area.
Here we overlay a raster of population density to proportionally allocate causalties geographically. This is done on a regular grid within the domain. Here we choose a grid cell resolution of .5 degree latitude and longitude.
Determine proportion of the tornado path that falls inside each cell. The rows index the tornado number and the columns index the cell number. Call these tornado segments. This takes about 5 minutes.
AreaTor = gArea(TornP3, byid = TRUE)
AreaCell = gArea(spT, byid = TRUE)
inter = gIntersects(spT, TornP3, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
Areas[[i]] = gArea(gIntersection(spT[ids[i, 2], ], TornP3[ids[i, 1], ]))
}
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "Area"
Areas.df$Cell = ids[, 2]
Areas.df$Torn = ids[, 1]Get population density. The raster population data are imported then cropped to match the domain & interpolated to match the resolution of the tornado grid cells. Population density is per square km. Check this. Next add population density values and cell areas to the rows in the Areas.df data frame. Next add per tornado casualties and tornado area in the cell to the columns of Areas.df. Finally create a weight column as a product of cell population and area of tornado within the cell.
library("RColorBrewer")
#download.file(url = "http://myweb.fsu.edu/jelsner/data/usadens.zip",
# destfile = "usadens.zip")
#unzip("usadens.zip")
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, r)
Pop[Pop == 0] = .01
Pop[is.na(Pop)] = .01
#pop = Pop
pop = resample(aggregate(Pop,
fact = c(nrow(Pop)/nrow(r), ncol(Pop)/ncol(r)),
fun = mean), r)
cellStats(pop, stat = "mean"); cellStats(pop, stat = "sd")## [1] 29.27627## [1] 63.82593range(log2(values(Pop)), na.rm = TRUE)## [1] -6.643856 11.716035rng = seq(-8, 12, 2)
cr = brewer.pal(9, "Blues")
cr[10] = "white"
cr = c(cr[10], cr[1:9])
labs = as.character(round(2^rng, 3))
p3a = rasterVis::levelplot(log2(Pop), margin = FALSE,
sub = expression(" 2000 Population Density [people per square km]"),
xlab = NULL, ylab = NULL,
col.regions = cr, at = rng,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15)))
p3a +
latticeExtra::layer(sp.polygons(bndrySP,
col = gray(.15), lwd = 1))
# Export at 8.5 x 6
Areas.df$pop = values(pop)[ids[, 2]]
Areas.df$AreaCell = AreaCell[ids[, 2]]
Areas.df$Tfat = TornP3$fat[ids[, 1]]
Areas.df$Tinj = TornP3$inj[ids[, 1]]
Areas.df$Tloss = TornP3$loss[ids[, 1]]
Areas.df$Year = TornP3$yr[ids[, 1]]
Areas.df$TArea = AreaTor[ids[, 1]]
Areas.df$weight = Areas.df$pop * Areas.df$AreaFigure 3: Population density from the Gridded Population of the World, version 3 For each tornado determine the total area-weighted population. Also determine the total proportion of each tornado inside the domain. Then for each tornado segment adjust the area-weighted population for proportion of tornado inside the domain and compute the proportion of casualties. Note: TArea[1] refers to the first instance of total tornado area as tornadoes often affect more than one cell. AreaCell is in square meters so divide by 1,000,000 before multiplying by population density (/sq km). Divide by number of years so the death rate per 100,000 people is per annum.
df = Areas.df %>%
group_by(Torn) %>%
summarize(tweight = sum(weight),
propin = sum(Area)/TArea[1])
Areas.df$tweight = (df$tweight/df$propin)[ids[, 1]]
Areas.df$propfat = Areas.df$Tfat * Areas.df$weight/Areas.df$tweight
Areas.df$propinj = Areas.df$Tinj * Areas.df$weight/Areas.df$tweight
Areas.df$proploss = Areas.df$Tloss * Areas.df$weight/Areas.df$tweight
# number of fatalities attributed to a given area or cell
df1 = Areas.df %>%
group_by(Cell) %>%
summarize(fat = sum(propfat),
inj = sum(propinj),
loss = sum(proploss),
fatrate = fat/(pop[1] * AreaCell[1]/1000000) * 100000/(eyr - syr + 1),
injrate = inj/(pop[1] * AreaCell[1]/1000000) * 100000/(eyr - syr + 1))Maps
library("rgeos")
EF0m = matrix(c(0, 0, 200, 200, 0, 0, 50, 50, 0, 0),
nrow = 5, ncol = 2)
EF1m = matrix(c(52.9, 52.9, 147.1, 147.1, 52.9, 13.225, 36.775, 36.775, 13.225, 13.225),
nrow = 5, ncol = 2)
EF2m = matrix(c(80, 80, 120, 120, 80, 20, 30, 30, 20, 20),
nrow = 5, ncol = 2)
EF3m = matrix(c(93.3, 93.3, 106.7, 106.7, 93.3, 23.325, 26.675, 26.675, 23.325, 23.325),
nrow = 5, ncol = 2)
trackm = matrix(c(-20, 25, 220, 25), nrow = 2, ncol = 2, byrow = TRUE)
theta = -pi/6
rot = matrix(c(cos(theta), -sin(theta), sin(theta), cos(theta)),
nrow = 2, ncol = 2, byrow = TRUE)
EF0r = EF0m %*% rot
EF1r = EF1m %*% rot
EF2r = EF2m %*% rot
EF3r = EF3m %*% rot
trackr = trackm %*% rot
EF0 = Polygons(list(Polygon(EF0r)), 0)
EF1 = Polygons(list(Polygon(EF1r)), 1)
EF2 = Polygons(list(Polygon(EF2r)), 2)
EF3 = Polygons(list(Polygon(EF3r)), 3)
spsNRC3 = SpatialPolygons(list(EF0, EF1, EF2, EF3))
spsNRC3A.df = fortify(spsNRC3)
spsNRC3A.df$EF = paste("EF", spsNRC3A.df$id, sep = "")
box1 = data.frame(long = c(-50, -50, 0, 0, -50),
lat = c(0, 50, 50, 0, 0))
box2 = data.frame(long = c(0, 0, 50, 50, 0),
lat = c(0, 50, 50, 0, 0))
box3 = data.frame(long = c(50, 50, 100, 100, 50),
lat = c(0, 50, 50, 0, 0))
box4 = data.frame(long = c(100, 100, 150, 150, 100),
lat = c(0, 50, 50, 0, 0))
box5 = data.frame(long = c(150, 150, 200, 200, 150),
lat = c(0, 50, 50, 0, 0))
box6 = data.frame(long = c(-50, -50, 0, 0, -50),
lat = c(50, 100, 100, 50, 50))
box7 = data.frame(long = c(0, 0, 50, 50, 0),
lat = c(50, 100, 100, 50, 50))
box8 = data.frame(long = c(50, 50, 100, 100, 50),
lat = c(50, 100, 100, 50, 50))
box9 = data.frame(long = c(100, 100, 150, 150, 100),
lat = c(50, 100, 100, 50, 50))
box10 = data.frame(long = c(150, 150, 200, 200, 150),
lat = c(50, 100, 100, 50, 50))
box11 = data.frame(long = c(-50, -50, 0, 0, -50),
lat = c(100, 150, 150, 100, 100))
box12 = data.frame(long = c(0, 0, 50, 50, 0),
lat = c(100, 150, 150, 100, 100))
box13 = data.frame(long = c(50, 50, 100, 100, 50),
lat = c(100, 150, 150, 100, 100))
box14 = data.frame(long = c(100, 100, 150, 150, 100),
lat = c(100, 150, 150, 100, 100))
box15 = data.frame(long = c(150, 150, 200, 200, 150),
lat = c(100, 150, 150, 100, 100))
box.df = rbind(box1, box2, box3, box4, box5, box6, box7, box8, box9, box10, box11, box12, box13, box14, box15)
box.df$pop = rep(c("< 5", "5-50", "101-1000", "5-50", "< 5", "< 5", "51-100", "> 1000", "51-100", "< 5", "< 5", "5-50", "101-1000", "5-50", "< 5"), each = 5)
box.df$id = rep(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15), each = 5)
lims = c("< 5", "5-50", "51-100", "101-1000", "> 1000")
brks = c("< 5", "5-50", "51-100", "101-1000", "> 1000")
vals = c("grey90", "grey80", "grey60", "grey40", "grey30")
ggplot(spsNRC3A.df[spsNRC3A.df$EF == "EF0", ], aes(x = long, y = lat)) +
scale_x_continuous(limits = c(-50, 200)) +
scale_y_continuous(limits = c(0, 150)) +
coord_fixed() +
geom_polygon(mapping = aes(long, lat, fill = factor(pop), group = id), data = box.df) +
geom_polygon(mapping = aes(long, lat, fill = factor(pop), group = id), data = box.df, color = "white", show_guide=FALSE) +
geom_polygon(fill = "transparent", color = "black") +
scale_fill_manual(name = expression("Population Density\n(people per square km)"), limits = lims, breaks = brks, values = vals) +
geom_segment(x = trackr[1, 1], xend = trackr[2, 1], y = trackr[1, 2], yend = trackr[2, 2],
arrow = grid::arrow(length = grid::unit(.6, "cm"), type = "closed"), size = 1.0, color = "black") +
xlab("") + ylab("") +
theme(panel.background = element_blank(),
panel.grid.major = element_line(colour = "gray", size = .25, linetype = 'dashed'),
panel.grid.minor = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
strip.background = element_blank()) +
annotate("text", x=-10, y=45, label= "0") +
annotate("text", x=-10, y=95, label= "0") +
annotate("text", x=-10, y=145, label= "0") +
annotate("text", x=40, y=45, label= "1") +
annotate("text", x=40, y=95, label= "2", color = "white") +
annotate("text", x=40, y=145, label= "0") +
annotate("text", x=90, y=45, label= "1", color = "white") +
annotate("text", x=90, y=95, label= "91", color = "white") +
annotate("text", x=90, y=145, label= "1", color = "white") +
annotate("text", x= 140, y=45, label= "0") +
annotate("text", x= 140, y=95, label= "3", color = "white") +
annotate("text", x= 140, y=145, label= "1") +
annotate("text", x= 190, y=45, label= "0") +
annotate("text", x= 190, y=95, label= "0") +
annotate("text", x= 190, y=145, label= "0")## Warning: `show_guide` has been deprecated. Please use `show.legend`
## instead.
Figure 4: Idealized tornado path and population density
Create rasters. fat is the number of deaths within the cell and fatrate is the annual fatality rate per 100,000. inj is the number of injuries within the cell. Subset on cell number ensures a cell without a tornado gets an NA. Total number of deaths in the SPC dataset should match (roughly) the total deaths across the cells.
fat = r
values(fat)[df1$Cell] = df1$fat
sum(TornP3$fat)## [1] 3553sum(df1$fat)## [1] 3525.413sum(TornP3$fat) - sum(df1$fat)## [1] 27.58721(sum(TornP3$fat) - sum(df1$fat))/sum(TornP3$fat) * 100## [1] 0.7764485fatrate = r
values(fatrate)[df1$Cell] = df1$fatrate
inj = r
values(inj)[df1$Cell] = df1$inj
sum(TornP3$inj)## [1] 56424sum(df1$inj)## [1] 55845.9sum(TornP3$inj) - sum(df1$inj)## [1] 578.0962(sum(TornP3$inj) - sum(df1$inj))/sum(TornP3$inj) * 100## [1] 1.024557loss = r
values(loss)[df1$Cell] = df1$loss
sum(TornP3$loss)## [1] 67235980sum(df1$loss)## [1] 62820866sum(TornP3$loss) - sum(df1$loss)## [1] 4415114(sum(TornP3$loss) - sum(df1$loss))/sum(TornP3$loss) * 100## [1] 6.566594Map the estimated number of deaths. The overlay is performed on the transformed spatial polygons data frame containing the tornado paths. It can take a few minutes to render because of the large number of tornado paths.
range(values(fat), na.rm = TRUE)## [1] 0.0000 116.1415range(values(log10(fat)), na.rm = TRUE)## [1] -Inf 2.064988rng = seq(0, 2.5, .5)
labs = as.character(round(10^rng, 0))
cr = brewer.pal(5, "Reds")
rasterVis::levelplot(log10(fat), margin = FALSE,
col.regions = cr, at = rng,
sub = expression(" Estimated Number of Tornado Deaths"),
xlab = NULL, ylab = NULL,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1)) +
latticeExtra::layer(sp.polygons(TornP3T, fill = gray(.4), alpha = .3))
Figure 5: Estimated number of tornado deaths over the period 1955 through 2016
Map the estimated number of injuries.
range(values(inj), na.rm = TRUE)## [1] 0.000 1500.269range(values(log10(inj)), na.rm = TRUE)## [1] -Inf 3.176169rng = seq(0, 3.5, .5)
labs = as.character(round(10^rng, 0))
cr = brewer.pal(7, "Purples")
rasterVis::levelplot(log10(inj), margin = FALSE,
col.regions = cr, at = rng,
sub =expression(" Estimated Number of Tornado Injuries"),
xlab = NULL, ylab = NULL,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1)) +
latticeExtra::layer(sp.polygons(TornP3T, fill = gray(.4), alpha = .3))
Figure 7: Estimated number of tornado injuries over the period 1955 through 2016
Map the estimated property loss.
range(values(loss), na.rm = TRUE)
range(values(log10(loss)), na.rm = TRUE)
rng = seq(0, 3.5, .5)
labs = as.character(round(10^rng, 0))
cr = brewer.pal(7, "Greens")
rasterVis::levelplot(log10(loss), margin = FALSE,
col.regions = cr, at = rng,
sub = expression(" Estimated Property Loss [millions of dollars]"),
xlab = NULL, ylab = NULL,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1)) #+
# latticeExtra::layer(sp.polygons(TornP3T, fill = gray(.4), alpha = .3))Mortality clusters identified using local Moran’s I.
LMfat = MoranLocal(fat)
range(values(LMfat), na.rm = TRUE)## [1] -1.198306 16.526946LMfat[LMfat > 1] = 1
rng = seq(-1, 1, .25)
labs = as.character(rng)
labs[length(labs)] = ">1"
cr = brewer.pal(8, "PiYG")
A = rasterVis::levelplot(LMfat, margin = FALSE,
col.regions = cr, at = rng,
sub = expression(" Local Moran's I [Deaths]"),
xlab = NULL, ylab = NULL,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1))
LMinj = MoranLocal(inj)
range(values(LMinj), na.rm = TRUE)## [1] -0.998788 16.940327LMinj[LMinj > 1] = 1
rng = seq(-1, 1, .25)
labs = as.character(rng)
labs[length(labs)] = ">1"
cr = brewer.pal(8, "PiYG")
B = rasterVis::levelplot(LMinj, margin = FALSE,
col.regions = cr, at = rng,
sub = expression(" Local Moran's I [Injuries]"),
xlab = NULL, ylab = NULL,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1))
A = update(A, main = textGrob("A", x = unit(.05, "npc"),
gp = gpar(fontsize=16)))
B = update(B, main = textGrob("B", x = unit(.05, "npc"),
gp = gpar(fontsize=16)))
library("gridExtra")##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinegrid.arrange(A, B, ncol = 1)
Figure 7: Casualty clusters using values from local Moran’s I for tornado deaths (A) and tornado injuries (B)
Raster stack of 10-yr non-overlapping injury maps
dt = 10
yrs = seq(1955, 2014, dt)
fatS = stack()
injS = stack()
injSc = stack()
for(yr in yrs){
df1 = Areas.df %>%
filter(Year >= yr & Year < (yr + dt)) %>%
group_by(Cell) %>%
summarize(fat = sum(propfat),
inj = sum(propinj))
ra = r
values(ra)[df1$Cell] = df1$fat
fatS = stack(fatS, ra)
values(ra)[df1$Cell] = df1$inj
injc = MoranLocal(ra)
injS = stack(injS, ra)
injSc = stack(injSc, injc)
}
names(fatS) = yrs
names(injS) = yrsrange(values(log10(injS)), na.rm = TRUE)## [1] -Inf 2.969235yrsL = c("1955-1964", "1965-1974", "1975-1984",
"1985-1994", "1995-2004", "2005-2014")
rng = seq(0, 3, .5)
labs = as.character(round(10^rng, 0))
cr = brewer.pal(6, "Purples")
A = rasterVis::levelplot(log10(injS), margin = FALSE,
col.regions = cr, at = rng,
sub = expression(" Estimated Number of Tornado Injuries"),
xlab = NULL, ylab = NULL,
names.attr = yrsL,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1))
injSc[injSc >= 1] = .999
injSc[injSc <= -1] = -.999
rng = seq(-1, 1, .25)
labs = as.character(rng)
labs[length(labs)] = ">1"
cr = brewer.pal(8, "PiYG")
B = rasterVis::levelplot(injSc, margin = FALSE,
col.regions = cr, at = rng,
sub = expression(" Local Moran's I [Injuries]"),
xlab = NULL, ylab = NULL,
names.attr = yrsL,
colorkey = list(space = 'bottom', labels = labs),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1))
A = update(A, main = textGrob("A", x = unit(.05, "npc"),
gp = gpar(fontsize=16)))
B = update(B, main = textGrob("B", x = unit(.05, "npc"),
gp = gpar(fontsize=16)))
grid.arrange(A, B, ncol = 2)
Figure 8: Estimated number of injuries (A) and injury clusters (B) for six non-overlapping decades starting in 1955
Map the mortality rate per 100,000 people. Smooth first.
fatrate3 = focal(fatrate, w = matrix(1/9, nrow = 3, ncol = 3))
fatrate5 = focal(fatrate, w = matrix(1/25, nrow = 5, ncol = 5))
fatrate7 = focal(fatrate, w = matrix(1/49, nrow = 7, ncol = 7))
range(values(log10(fatrate5)), na.rm = TRUE)
range(values(log10(fatrate3)), na.rm = TRUE)
range(values(log10(fatrate7)), na.rm = TRUE)
range(values(fatrate), na.rm = TRUE)
#rng = seq(-6, -3, .5)
rng = seq(-3, .5, .5)
cr = brewer.pal(7, "Reds")
levelplot(log10(fatrate5), margin = FALSE,
col.regions = cr, at = rng,
sub = expression(" Tornado Mortality Rate Per 100,000 People (log10 scale) "),
xlab = NULL, ylab = NULL,
colorkey = list(space = 'bottom'),
par.settings = list(fontsize = list(text = 15))) +
latticeExtra::layer(sp.polygons(bndrySP, col = gray(.15), lwd = 1))
# latticeExtra::layer(sp.polygons(TornP3T, fill = gray(.4), alpha = .3))Validation
Get administrative boundaries. 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.
#download.file("http://www2.census.gov/geo/tiger/GENZ2013/cb_2013_us_county_5m.zip",
# "cb_2013_us_county_5m.zip", mode = "wb")
unzip("cb_2013_us_county_5m.zip")
library("rgdal")## rgdal: version: 1.2-7, (SVN revision 660)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/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.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-4US.sp = readOGR(dsn = ".", layer = "cb_2013_us_county_5m",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "cb_2013_us_county_5m"
## with 3234 features
## It has 9 fields
## Integer64 fields read as strings: ALAND AWATERExtract Tennessee using the state Federal Information Processing Standard (FIPS).
FP = 47
AB = "TN"
ST.sp = US.sp[US.sp$STATEFP == FP, ]
nrow(ST.sp)## [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)))
rgeos::gArea(counties.sp) / 10^6## [1] 108310.3Overlay paths onto the domain extent and return a vector indicating either NA (tornado is not inside extent) or the cell number. Subset the tornadoes by the raster extent (domain). Remove duplicate paths from the data frame.
tc = over(TornP, counties.sp)
TornP2 = subset(TornP, !is.na(tc) & 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.3505697TornP3 = subset(TornP2, !dup)
dim(TornP3@data)## [1] 1137 22Create a latex table
tA = table(TornP3$fat, TornP3$mag)[-1, ]
tA = rbind(tA[1:5,], colSums(tA[6:dim(tA)[1],]))
rownames(tA)[6] = "6+"
xtable(tA, digits = 0)
tB = table(TornP3$inj, TornP3$mag)[-1, ]
tB = rbind(tB[1:5,], colSums(tB[6:dim(tB)[1],]))
rownames(tB)[6] = "6+"
xtable(tB, digits = 0)There are 1137 tornado paths that intersect the state. The table shows the distributon of fatalities by EF rating.
Use ggmap. This requires the spatial polygons data frame be converted to the geographic projection of the raster.
TornP3T = spTransform(TornP3, CRS(projection(r)))
mapdata = fortify(TornP3T[TornP3T$fat > 0, ])## Regions defined for each PolygonsMap = ggmap::get_map(location = "Nashville",
source = "stamen",
zoom = 6,
color = "bw",
maptype = "toner")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Nashville&zoom=6&size=640x640&scale=2&maptype=terrain&sensor=false## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Nashville&sensor=false## Map from URL : http://tile.stamen.com/toner/6/15/23.png## Map from URL : http://tile.stamen.com/toner/6/16/23.png## Map from URL : http://tile.stamen.com/toner/6/17/23.png## Map from URL : http://tile.stamen.com/toner/6/15/24.png## Map from URL : http://tile.stamen.com/toner/6/16/24.png## Map from URL : http://tile.stamen.com/toner/6/17/24.png## Map from URL : http://tile.stamen.com/toner/6/15/25.png## Map from URL : http://tile.stamen.com/toner/6/16/25.png## Map from URL : http://tile.stamen.com/toner/6/17/25.png## Map from URL : http://tile.stamen.com/toner/6/15/26.png## Map from URL : http://tile.stamen.com/toner/6/16/26.png## Map from URL : http://tile.stamen.com/toner/6/17/26.pngggmap(Map, dev = "extent") +
geom_polygon(aes(x = long, y = lat, group = group),
color = "black", fill = "black",
data = mapdata) +
labs(x = NULL, y = NULL) +
# labs(x = expression(paste("Longitude (", degree, "W)")),
# y = expression(paste("Latitude (", degree, "N)"))) +
theme(legend.position = "none")
Determine the number of tornado paths intersecting each county.
ct = over(counties.sp, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
ID = sapply(slot(counties.sp, "polygons"), function(x) slot(x, "ID"))
df = data.frame(nT)
rownames(df) = ID
spdf = SpatialPolygonsDataFrame(counties.sp, df)
sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)## [1] 1375## [1] 14.47368## [1] 1## [1] 54## [1] 9.774959## [1] 6.601625The average number of tornadoes per cell is 14.5, with a minimum of 1 and a maximum of 54 tornadoes. There are 203 cells that have between 35 and 45 tornado paths. The standard deviation of the tornadoes is 9.8 with a variance to mean ratio (VMR) of 6.6. Map number of tornadoes
lcc = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-87 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
spdf = spTransform(spdf, CRS(lcc))
range(nT)## [1] 1 54rng = seq(0, 60, 10)
cr = wesanderson::wes_palette(6, name = "Zissou", type = "continuous")
spplot(spdf, "nT", col = "gray", at = rng,
col.regions = cr,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = expression("Number of Tornadoes (1955-2014)"))
Determine proportion of the tornado path in square meters that falls inside each county. The rows index the tornado number and the columns index the county. Call these tornado segments.
AreaTor = gArea(TornP3, byid = TRUE)
AreaCounty = gArea(counties.sp, byid = TRUE)
inter = gIntersects(counties.sp, TornP3, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
Areas[[i]] = gArea(gIntersection(counties.sp[ids[i, 2], ], TornP3[ids[i, 1], ]))
}
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "Area"
Areas.df$County = ids[, 2]
Areas.df$Torn = ids[, 1]Determine the mean population density per county.
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, countiesun.sp)
PopT = projectRaster(Pop, crs = crs(counties.sp))
pop = raster::extract(PopT, counties.sp, fun = mean, na.rm = TRUE)[, 1]
Areas.df$pop = pop[ids[, 2]]
Areas.df$AreaCounty = AreaCounty[ids[, 2]]
Areas.df$Tfat = TornP3$fat[ids[, 1]]
Areas.df$Tinj = TornP3$inj[ids[, 1]]
Areas.df$TArea = AreaTor[ids[, 1]]
Areas.df$weight = Areas.df$pop * Areas.df$Areadf = Areas.df %>%
group_by(Torn) %>%
summarize(tweight = sum(weight),
propin = sum(Area)/TArea[1])
Areas.df$tweight = (df$tweight/df$propin)[ids[, 1]]
Areas.df$propfat = Areas.df$Tfat * Areas.df$weight/Areas.df$tweight
Areas.df$propinj = Areas.df$Tinj * Areas.df$weight/Areas.df$tweight
# number of fatalities attributed to a given area or county
df1 = Areas.df %>%
group_by(County) %>%
summarize(propfat = sum(propfat),
propinj = sum(propinj),
propfatratio = sum(propfat)/(pop[1] * AreaCounty[1]),
propinjratio = sum(propinj)/(pop[1] * AreaCounty[1]))
sum(TornP3$inj)## [1] 5169sum(df1$propinj)## [1] 3855.013Map estimated fatalities (Tennessee)
mean(df1$propfat)## [1] 2.898484range(df1$propinj)## [1] 0.0000 418.6885range(df1$propfat)## [1] 0.00000 19.77719spdf$propfat = df1$propfat
spdf$propinj = df1$propinj
rng = seq(0, 450, 50)
rng = seq(0, 20, 4)
cr = brewer.pal(5, "Reds")
A = spplot(spdf, "propfat", col = "gray", at = rng,
col.regions = cr,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = expression("Estimated Number of Tornado Deaths"))Compare with county-level casualty data from NCDC Storm Data. http://www.ncdc.noaa.gov/stormevents/ Search Tennessee > All Counties > Sort by Deaths/Inj Clean data in google sheets with add-on to remove cells. Duplicate rows with multiple county names. Read the cleaned data then sum the number of casualties by county.
Casualties = read.csv("storm_data_search_results.csv", header = TRUE)
table(Casualties$CZ_NAME_STR, Casualties$DEATHS_DIRECT)##
## 0 1 2 3 4 5 6 7 8 10 11 13 16
## BEDFORD CO. 4 1 0 0 0 0 0 0 0 0 0 0 0
## BENTON CO. 0 1 0 0 0 0 0 0 0 0 0 0 0
## BLEDSOE CO. 0 0 0 0 1 0 0 0 0 0 0 0 0
## BLOUNT CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## BRADLEY CO. 8 4 1 0 2 0 0 0 0 0 0 0 0
## CANNON CO. 1 2 0 0 0 0 0 0 0 0 0 0 0
## CARROLL CO. 4 2 1 1 0 0 0 0 0 0 0 0 0
## CHEATHAM CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## CHESTER CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## CLAIBORNE CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## COCKE CO. 1 1 0 0 0 0 0 0 0 0 0 0 0
## COFFEE CO. 1 0 1 0 0 0 0 0 0 0 0 0 0
## CROCKETT CO. 3 0 0 1 0 0 0 0 0 0 0 0 0
## CUMBERLAND CO. 3 2 1 0 1 0 0 0 0 0 0 0 0
## DAVIDSON CO. 7 2 0 0 0 0 0 0 0 0 0 0 0
## DEKALB CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## DICKSON CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## DYER CO. 4 1 1 0 0 0 0 0 0 0 0 0 1
## FAYETTE CO. 5 1 0 2 0 0 0 0 0 0 0 0 0
## FENTRESS CO. 3 0 0 0 0 0 0 1 0 0 0 0 0
## FRANKLIN CO. 1 1 0 0 0 1 0 0 0 0 0 0 0
## GIBSON CO. 3 1 1 0 0 0 1 0 0 0 0 0 0
## GILES CO. 4 1 0 1 0 0 0 0 0 0 0 0 0
## GRAINGER CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## GREENE CO. 2 2 0 0 0 0 1 0 0 0 0 0 0
## GRUNDY CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## HAMILTON CO. 4 0 0 0 0 0 0 0 1 0 0 0 0
## HARDEMAN CO. 1 3 0 0 0 0 0 0 0 0 0 0 0
## HARDIN CO. 3 0 0 1 0 0 0 0 0 0 0 0 0
## HAWKINS CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## HAYWOOD CO. 2 1 0 0 0 0 0 0 0 0 0 0 0
## HENDERSON CO. 1 0 0 1 0 0 0 0 0 0 0 0 0
## HENRY CO. 4 0 1 0 0 0 0 0 0 0 0 0 0
## HICKMAN CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## HUMPHREYS CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## JACKSON CO. 1 0 1 0 0 0 0 0 0 0 0 0 0
## JOHNSON CO. 1 0 1 0 0 0 0 0 0 0 0 0 0
## KNOX CO. 3 0 1 0 0 0 0 0 0 0 0 0 0
## LAUDERDALE CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## LAWRENCE CO. 10 0 0 1 0 0 0 0 0 0 0 0 0
## LEWIS CO. 1 1 0 0 0 0 0 0 0 0 0 0 0
## LINCOLN CO. 1 0 1 0 0 0 1 0 0 0 0 0 0
## LOUDON CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## MACON CO. 0 0 0 0 0 0 0 0 0 0 0 1 0
## MADISON CO. 5 1 1 0 0 0 1 0 0 0 1 0 0
## MARION CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## MARSHALL CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## MAURY CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## MCMINN CO. 8 2 0 0 0 0 0 0 0 0 0 0 0
## MCNAIRY CO. 5 0 0 0 1 0 0 0 0 0 0 0 0
## MONROE CO. 7 0 0 0 0 0 0 0 0 0 0 0 0
## MONTGOMERY CO. 5 0 1 0 0 0 0 0 0 0 0 0 0
## MOORE CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## MORGAN CO. 1 0 0 0 0 0 0 1 0 0 0 0 0
## OBION CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## OVERTON CO. 3 0 0 1 0 0 0 0 0 0 0 0 0
## PERRY CO. 0 0 0 1 0 0 0 0 0 0 0 0 0
## PICKETT CO. 1 0 0 0 0 1 0 0 0 0 0 0 0
## POLK CO. 2 1 0 0 0 0 0 0 0 0 0 0 0
## PUTNAM CO. 2 1 0 0 0 0 0 0 0 1 0 0 0
## RHEA CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## ROANE CO. 0 1 0 0 0 0 0 0 0 0 0 0 0
## ROBERTSON CO. 3 0 1 0 0 0 0 0 0 0 0 0 0
## RUTHERFORD CO. 5 0 1 0 0 0 0 0 0 0 0 0 0
## SCOTT CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## SEQUATCHIE CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## SHELBY CO. 12 0 0 2 0 0 0 0 0 0 0 0 0
## SMITH CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## STEWART CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## SULLIVAN CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## SUMNER CO. 7 1 0 0 0 0 0 2 0 0 0 0 0
## TIPTON CO. 3 3 0 0 1 0 0 0 0 0 0 0 0
## TROUSDALE CO. 0 0 1 0 0 0 0 0 0 0 0 0 0
## UNICOI CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## WARREN CO. 3 2 1 0 0 0 0 0 0 0 0 0 0
## WASHINGTON CO. 0 1 0 0 0 0 0 0 0 0 0 0 0
## WAYNE CO. 4 1 0 1 0 0 0 0 0 0 0 0 0
## WEAKLEY CO. 6 0 0 0 0 0 0 0 0 0 0 0 0
## WHITE CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## WILLIAMSON CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## WILSON CO. 4 0 0 0 0 0 0 0 0 0 0 0 0AC = Casualties %>%
group_by(CZ_NAME_STR) %>%
summarize(nD = sum(DEATHS_DIRECT),
nI = sum(INJURIES_DIRECT))
AC[1:10, ]; AC[11:20, ]; AC[21:30, ]## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 BEDFORD CO. 1 9
## 2 BENTON CO. 1 5
## 3 BLEDSOE CO. 4 10
## 4 BLOUNT CO. 1 62
## 5 BRADLEY CO. 14 425
## 6 CANNON CO. 2 7
## 7 CARROLL CO. 7 156
## 8 CHEATHAM CO. 0 1
## 9 CHESTER CO. 0 1
## 10 CLAIBORNE CO. 0 18## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 COCKE CO. 1 11
## 2 COFFEE CO. 2 25
## 3 CROCKETT CO. 3 18
## 4 CUMBERLAND CO. 8 73
## 5 DAVIDSON CO. 2 104
## 6 DEKALB CO. 1 22
## 7 DICKSON CO. 0 5
## 8 DYER CO. 19 101
## 9 FAYETTE CO. 7 76
## 10 FENTRESS CO. 7 157## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 FRANKLIN CO. 6 22
## 2 GIBSON CO. 9 56
## 3 GILES CO. 4 41
## 4 GRAINGER CO. 0 2
## 5 GREENE CO. 8 114
## 6 GRUNDY CO. 0 6
## 7 HAMILTON CO. 8 182
## 8 HARDEMAN CO. 3 11
## 9 HARDIN CO. 3 12
## 10 HAWKINS CO. 0 6Match county names.
NMc = as.character(AC$CZ_NAME_STR)
NMc = unlist(strsplit(NMc, split = " CO."))
AC$NAME = NMc
#AC = AC[order(NMc), ]
nm = data.frame(NAME = toupper(ST.sp$NAME[order(ST.sp$NAME)]), stringsAsFactors=FALSE)
AC2 = left_join(nm, AC[, -1])## Joining, by = "NAME"AC2[is.na(AC2)] = 0Get the list of polygons IDs. Create a data frame of the IDs and the county names with rows alphabetical by county name.
IDs = as.integer(sapply(slot(spdf, "polygons"), function(x) slot(x, "ID")))
Names = toupper(ST.sp$NAME)
IDs = IDs[order(Names)]
df = data.frame(ID = ST.sp$GEOID[order(Names)],
Name = Names[order(Names)],
AC2,
spdf@data[order(Names), ])
cor(df$nD, df$nT, method = "spearman")## [1] 0.5420071cor(df$nI, df$nT, method = "spearman")## [1] 0.7025905cor(df$nD, df$propfat, method = "spearman")## [1] 0.8623928cor(df$nI, df$propinj, method = "spearman")## [1] 0.875547sum(abs(df$nD - df$propfat))/length(df$nD)## [1] 1.183744sum(abs(df$nI - df$propinj))/length(df$nI)## [1] 12.11867Attach the observed casualties to the spatial polygons data frame and make a map.
IDs = as.integer(sapply(slot(spdf, "polygons"), function(x) slot(x, "ID")))
Match = match(IDs, rownames(df))
spdf$nD = df$nD[Match]
spdf$nI = df$nI[Match]
range(spdf$nD)## [1] 0 20rng = seq(0, 20, 4)
spdf$nD[which.max(spdf$nD)] = spdf$nD[which.max(spdf$nD)] - .001
cr = brewer.pal(5, "Reds")
B = spplot(spdf, "nD", col = "gray", at = rng,
col.regions = cr,
colorkey = list(space = "bottom", labels = paste(rng)),
par.settings = list(axis.line = list(col = NA)),
sub = expression("Actual Number of Tornado Deaths"))
A = update(A, main = textGrob("A", x = unit(.05, "npc"), gp=gpar(fontsize=16)))
B = update(B, main = textGrob("B", x = unit(.05, "npc"), gp=gpar(fontsize=16)))
library("gridExtra")
grid.arrange(A, B, ncol = 1)
Figure 9: Estimated deaths (A) and actual deaths (B) from tornadoes by county in the state of Tennessee
Try other states: Kansas, Missouri, Kentucky, Oklahoma, Arkansas
Begin with Kansas
FP = 20
AB = "KS"
ST.sp = US.sp[US.sp$STATEFP == FP, ]
nrow(ST.sp)## [1] 105county = paste(ST.sp$STATEFP, ST.sp$COUNTYFP, sep = "")
countiesun.sp = geometry(spChFIDs(ST.sp, county))
counties.sp = spTransform(countiesun.sp, CRS(proj4string(TornL)))
rgeos::gArea(counties.sp) / 10^6## [1] 210845.3tc = over(TornP, counties.sp)
TornP2 = subset(TornP, !is.na(tc) & 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] 74sum(dup)/dim(TornP2@data)[1] * 100## [1] 1.886792TornP3 = subset(TornP2, !dup)
dim(TornP3@data)## [1] 3848 22TornP3T = spTransform(TornP3, CRS(projection(r)))
mapdata = fortify(TornP3T[TornP3T$fat > 0, ])## Regions defined for each Polygonsct = over(counties.sp, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
ID = sapply(slot(counties.sp, "polygons"), function(x) slot(x, "ID"))
df = data.frame(nT)
rownames(df) = ID
spdf = SpatialPolygonsDataFrame(counties.sp, df)
sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)## [1] 4360## [1] 41.52381## [1] 7## [1] 95## [1] 18.23885## [1] 8.011203AreaTor = gArea(TornP3, byid = TRUE)
AreaCounty = gArea(counties.sp, byid = TRUE)
inter = gIntersects(counties.sp, TornP3, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
Areas[[i]] = gArea(gIntersection(counties.sp[ids[i, 2], ], TornP3[ids[i, 1], ]))
}
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "Area"
Areas.df$County = ids[, 2]
Areas.df$Torn = ids[, 1]
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, countiesun.sp)
PopT = projectRaster(Pop, crs = crs(counties.sp))
pop = raster::extract(PopT, counties.sp, fun = mean, na.rm = TRUE)[, 1]
Areas.df$pop = pop[ids[, 2]]
Areas.df$AreaCounty = AreaCounty[ids[, 2]]
Areas.df$Tfat = TornP3$fat[ids[, 1]]
Areas.df$Tinj = TornP3$inj[ids[, 1]]
Areas.df$TArea = AreaTor[ids[, 1]]
Areas.df$weight = Areas.df$pop * Areas.df$Area
df = Areas.df %>%
group_by(Torn) %>%
summarize(tweight = sum(weight),
propin = sum(Area)/TArea[1])
Areas.df$tweight = (df$tweight/df$propin)[ids[, 1]]
Areas.df$propfat = Areas.df$Tfat * Areas.df$weight/Areas.df$tweight
Areas.df$propinj = Areas.df$Tinj * Areas.df$weight/Areas.df$tweight
# number of fatalities attributed to a given area or county
df1 = Areas.df %>%
group_by(County) %>%
summarize(propfat = sum(propfat),
propinj = sum(propinj),
propfatratio = sum(propfat)/(pop[1] * AreaCounty[1]),
propinjratio = sum(propinj)/(pop[1] * AreaCounty[1]))
sum(TornP3$inj)## [1] 3590sum(df1$propinj)## [1] 2667.525spdf$propfat = df1$propfat
spdf$propinj = df1$propinj
Casualties = read.csv("storm_data_search_results_Kansas.csv", header = TRUE)
table(Casualties$CZ_NAME_STR, Casualties$DEATHS_DIRECT)##
## 0 1 2 3 4 5 6 11 13 15 16 75
## ALLEN CO. 2 0 0 0 0 0 0 0 0 0 0 0
## ANDERSON CO. 3 0 0 1 0 0 0 0 0 0 0 0
## ATCHISON CO. 2 0 0 0 0 0 0 0 0 0 0 0
## BARBER CO. 3 0 0 0 0 0 0 0 0 0 0 0
## BARTON CO. 13 2 0 0 0 0 0 0 0 0 0 0
## BOURBON CO. 3 0 0 0 0 0 0 0 0 0 0 0
## BROWN CO. 1 0 0 0 0 0 0 0 0 0 0 0
## BUTLER CO. 10 0 0 0 0 0 0 0 1 1 0 0
## CHASE CO. 2 0 0 0 0 0 0 0 0 0 0 0
## CHAUTAUQUA CO. 2 0 0 0 0 0 0 0 0 0 0 0
## CHEROKEE CO. 6 1 0 1 0 0 0 0 0 0 0 0
## CHEYENNE CO. 1 0 0 0 0 0 0 0 0 0 0 0
## CLARK CO. 3 0 0 0 0 0 0 0 0 0 0 0
## CLAY CO. 5 0 0 0 0 0 0 0 0 0 0 0
## CLOUD CO. 8 1 0 0 0 0 0 0 0 0 0 0
## COFFEY CO. 3 0 0 0 0 0 0 0 0 0 0 0
## COMANCHE CO. 4 0 0 0 0 0 0 0 0 0 0 0
## COWLEY CO. 10 2 0 0 0 0 0 0 0 0 0 1
## CRAWFORD CO. 6 1 0 1 0 0 0 0 0 0 0 0
## DECATUR CO. 2 0 0 0 0 0 0 0 0 0 0 0
## DICKINSON CO. 3 1 0 0 0 0 0 0 0 0 0 0
## DONIPHAN CO. 3 0 0 0 0 0 0 0 0 0 0 0
## DOUGLAS CO. 6 1 0 0 0 0 0 0 0 0 0 0
## EDWARDS CO. 5 0 0 0 0 0 0 0 0 0 0 0
## ELK CO. 6 2 0 0 0 0 0 0 0 0 0 0
## ELLIS CO. 5 0 0 0 0 0 0 0 0 0 0 0
## ELLSWORTH CO. 3 0 0 0 0 0 0 0 0 0 0 0
## FINNEY CO. 9 1 0 0 0 0 0 0 0 0 0 0
## FORD CO. 10 0 0 0 0 0 0 0 0 0 0 0
## FRANKLIN CO. 6 0 0 1 0 0 0 0 0 0 0 0
## GEARY CO. 2 0 0 0 0 0 0 0 0 0 0 0
## GOVE CO. 2 0 0 0 0 0 0 0 0 0 0 0
## GRAHAM CO. 4 0 0 0 0 0 0 0 0 0 0 0
## GRANT CO. 1 0 0 0 0 0 0 0 0 0 0 0
## GRAY CO. 4 0 0 0 0 0 0 0 0 0 0 0
## GREELEY CO. 1 0 0 0 0 0 0 0 0 0 0 0
## GREENWOOD CO. 7 0 0 0 0 0 0 0 0 0 0 0
## HAMILTON CO. 2 0 0 0 0 0 0 0 0 0 0 0
## HARPER CO. 5 0 0 0 0 0 0 0 0 0 0 0
## HARVEY CO. 7 1 0 0 0 0 0 0 0 0 0 0
## HASKELL CO. 1 0 0 0 0 0 0 0 0 0 0 0
## HODGEMAN CO. 6 0 0 0 0 0 0 0 0 0 0 0
## JACKSON CO. 4 1 0 1 0 0 0 0 0 0 0 0
## JEFFERSON CO. 5 0 0 0 0 0 0 0 0 0 0 0
## JEWELL CO. 3 0 0 0 0 0 0 0 0 0 0 0
## JOHNSON CO. 6 0 0 0 0 0 0 0 0 0 0 0
## KINGMAN CO. 2 0 0 0 0 0 0 0 0 0 0 0
## KIOWA CO. 5 0 0 0 0 0 0 1 0 0 0 0
## LABETTE CO. 4 1 0 0 0 0 0 0 0 0 0 0
## LANE CO. 1 0 0 0 0 0 0 0 0 0 0 0
## LEAVENWORTH CO. 3 2 0 0 0 0 0 0 0 0 0 0
## LINCOLN CO. 2 0 0 0 0 0 0 0 0 0 0 0
## LINN CO. 2 0 0 0 0 0 0 0 0 0 0 0
## LOGAN CO. 3 0 0 0 0 0 0 0 0 0 0 0
## LYON CO. 5 1 0 0 0 0 1 0 0 0 0 0
## MARION CO. 10 1 0 0 0 0 0 0 0 0 0 0
## MARSHALL CO. 3 0 0 0 0 0 0 0 0 0 0 0
## MCPHERSON CO. 8 0 0 0 0 0 0 0 0 0 0 0
## MEADE CO. 8 0 0 0 0 0 0 0 0 0 0 0
## MIAMI CO. 5 0 0 0 1 0 0 0 0 0 0 0
## MITCHELL CO. 6 0 0 0 0 0 0 0 0 0 0 0
## MONTGOMERY CO. 2 1 0 0 0 0 0 0 0 0 0 0
## MORRIS CO. 3 0 0 0 0 0 0 0 0 0 0 0
## MORTON CO. 1 0 0 0 0 0 0 0 0 0 0 0
## NEMAHA CO. 2 0 0 0 0 0 0 0 0 0 0 0
## NEOSHO CO. 3 0 0 0 0 0 0 0 0 0 0 0
## NESS CO. 5 0 0 0 0 0 0 0 0 0 0 0
## NORTON CO. 1 0 0 0 0 0 0 0 0 0 0 0
## OSAGE CO. 2 1 0 0 0 0 0 0 0 0 1 0
## OSBORNE CO. 4 0 0 0 0 0 0 0 0 0 0 0
## OTTAWA CO. 2 2 0 0 0 0 0 0 0 0 0 0
## PAWNEE CO. 2 0 0 0 0 0 0 0 0 0 0 0
## PHILLIPS CO. 3 0 0 0 0 0 0 0 0 0 0 0
## POTTAWATOMIE CO. 3 0 0 0 0 0 0 0 0 0 0 0
## PRATT CO. 5 1 1 0 0 0 0 0 0 0 0 0
## RAWLINS CO. 3 0 0 0 0 0 0 0 0 0 0 0
## RENO CO. 11 0 0 0 0 0 0 0 0 0 0 0
## REPUBLIC CO. 2 0 0 0 0 0 0 0 0 0 0 0
## RICE CO. 10 0 0 0 0 0 0 0 0 0 0 0
## RILEY CO. 2 0 0 0 0 0 0 0 0 0 0 0
## ROOKS CO. 5 0 0 0 0 0 0 0 0 0 0 0
## RUSH CO. 1 0 0 0 0 0 0 0 0 0 0 0
## RUSSELL CO. 4 1 0 0 0 0 0 0 0 0 0 0
## SALINE CO. 4 0 0 0 0 0 0 0 0 0 0 0
## SCOTT CO. 1 1 0 0 0 0 0 0 0 0 0 0
## SEDGWICK CO. 18 0 0 1 1 0 1 0 0 0 0 0
## SEWARD CO. 5 0 0 0 0 0 0 0 0 0 0 0
## SHAWNEE CO. 9 2 0 0 0 0 0 0 0 0 1 0
## SHERIDAN CO. 1 0 0 0 0 0 0 0 0 0 0 0
## SHERMAN CO. 4 0 0 0 0 0 0 0 0 0 0 0
## SMITH CO. 5 0 0 0 0 0 0 0 0 0 0 0
## STAFFORD CO. 16 1 1 0 0 0 0 0 0 0 0 0
## SUMNER CO. 11 0 0 0 0 1 0 0 0 0 0 0
## THOMAS CO. 8 0 0 0 0 0 0 0 0 0 0 0
## TREGO CO. 3 0 0 0 0 0 0 0 0 0 0 0
## WABAUNSEE CO. 2 1 0 0 0 0 0 0 0 0 0 0
## WALLACE CO. 3 0 0 0 0 0 0 0 0 0 0 0
## WASHINGTON CO. 2 0 1 0 0 0 0 0 0 0 0 0
## WICHITA CO. 3 0 0 0 0 0 0 0 0 0 0 0
## WOODSON CO. 3 0 0 0 0 0 0 0 0 0 0 0
## WYANDOTTE CO. 3 0 1 0 0 0 0 0 0 0 0 0AC = Casualties %>%
group_by(CZ_NAME_STR) %>%
summarize(nD = sum(DEATHS_DIRECT),
nI = sum(INJURIES_DIRECT))
AC[1:10, ]; AC[11:20, ]; AC[21:30, ]## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 ALLEN CO. 0 4
## 2 ANDERSON CO. 3 12
## 3 ATCHISON CO. 0 10
## 4 BARBER CO. 0 2
## 5 BARTON CO. 2 42
## 6 BOURBON CO. 0 0
## 7 BROWN CO. 0 5
## 8 BUTLER CO. 28 170
## 9 CHASE CO. 0 0
## 10 CHAUTAUQUA CO. 0 0## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 CHEROKEE CO. 4 56
## 2 CHEYENNE CO. 0 0
## 3 CLARK CO. 0 0
## 4 CLAY CO. 0 31
## 5 CLOUD CO. 1 7
## 6 COFFEY CO. 0 5
## 7 COMANCHE CO. 0 1
## 8 COWLEY CO. 77 290
## 9 CRAWFORD CO. 4 43
## 10 DECATUR CO. 0 5## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 DICKINSON CO. 1 12
## 2 DONIPHAN CO. 0 2
## 3 DOUGLAS CO. 1 45
## 4 EDWARDS CO. 0 7
## 5 ELK CO. 2 8
## 6 ELLIS CO. 0 7
## 7 ELLSWORTH CO. 0 0
## 8 FINNEY CO. 1 41
## 9 FORD CO. 0 0
## 10 FRANKLIN CO. 3 12NMc = as.character(AC$CZ_NAME_STR)
NMc = unlist(strsplit(NMc, split = " CO."))
AC$NAME = NMc
#AC = AC[order(NMc), ]
nm = data.frame(NAME = toupper(ST.sp$NAME[order(ST.sp$NAME)]), stringsAsFactors=FALSE)
AC2 = left_join(nm, AC[, -1])## Joining, by = "NAME"AC2[is.na(AC2)] = 0
IDs = as.integer(sapply(slot(spdf, "polygons"), function(x) slot(x, "ID")))
Names = toupper(ST.sp$NAME)
IDs = IDs[order(Names)]
df = data.frame(ID = ST.sp$GEOID[order(Names)],
Name = Names[order(Names)],
AC2,
spdf@data[order(Names), ])
cor(df$nD, df$nT, method = "spearman")## [1] 0.1884472cor(df$nI, df$nT, method = "spearman")## [1] 0.1389659cor(df$nD, df$propfat, method = "spearman")## [1] 0.7943569cor(df$nI, df$propinj, method = "spearman")## [1] 0.9174576sum(abs(df$nD - df$propfat))/length(df$nD)## [1] 0.9891641sum(abs(df$nI - df$propinj))/length(df$nI)## [1] 9.023746Missouri
FP = 29
AB = "MO"
ST.sp = US.sp[US.sp$STATEFP == FP, ]
nrow(ST.sp)## [1] 115county = paste(ST.sp$STATEFP, ST.sp$COUNTYFP, sep = "")
countiesun.sp = geometry(spChFIDs(ST.sp, county))
counties.sp = spTransform(countiesun.sp, CRS(proj4string(TornL)))
rgeos::gArea(counties.sp) / 10^6## [1] 178687tc = over(TornP, counties.sp)
TornP2 = subset(TornP, !is.na(tc) & 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] 7sum(dup)/dim(TornP2@data)[1] * 100## [1] 0.3237743TornP3 = subset(TornP2, !dup)
dim(TornP3@data)## [1] 2155 22TornP3T = spTransform(TornP3, CRS(projection(r)))
mapdata = fortify(TornP3T[TornP3T$fat > 0, ])## Regions defined for each Polygonsct = over(counties.sp, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
ID = sapply(slot(counties.sp, "polygons"), function(x) slot(x, "ID"))
df = data.frame(nT)
rownames(df) = ID
spdf = SpatialPolygonsDataFrame(counties.sp, df)
sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)## [1] 2497## [1] 21.71304## [1] 7## [1] 48## [1] 9.634337## [1] 4.27487AreaTor = gArea(TornP3, byid = TRUE)
AreaCounty = gArea(counties.sp, byid = TRUE)
inter = gIntersects(counties.sp, TornP3, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
Areas[[i]] = gArea(gIntersection(counties.sp[ids[i, 2], ], TornP3[ids[i, 1], ]))
}
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "Area"
Areas.df$County = ids[, 2]
Areas.df$Torn = ids[, 1]
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, countiesun.sp)
PopT = projectRaster(Pop, crs = crs(counties.sp))
pop = raster::extract(PopT, counties.sp, fun = mean, na.rm = TRUE)[, 1]
Areas.df$pop = pop[ids[, 2]]
Areas.df$AreaCounty = AreaCounty[ids[, 2]]
Areas.df$Tfat = TornP3$fat[ids[, 1]]
Areas.df$Tinj = TornP3$inj[ids[, 1]]
Areas.df$TArea = AreaTor[ids[, 1]]
Areas.df$weight = Areas.df$pop * Areas.df$Area
df = Areas.df %>%
group_by(Torn) %>%
summarize(tweight = sum(weight),
propin = sum(Area)/TArea[1])
Areas.df$tweight = (df$tweight/df$propin)[ids[, 1]]
Areas.df$propfat = Areas.df$Tfat * Areas.df$weight/Areas.df$tweight
Areas.df$propinj = Areas.df$Tinj * Areas.df$weight/Areas.df$tweight
# number of fatalities attributed to a given area or county
df1 = Areas.df %>%
group_by(County) %>%
summarize(propfat = sum(propfat),
propinj = sum(propinj),
propfatratio = sum(propfat)/(pop[1] * AreaCounty[1]),
propinjratio = sum(propinj)/(pop[1] * AreaCounty[1]))
sum(TornP3$inj)## [1] 4866sum(df1$propinj)## [1] 4123.331spdf$propfat = df1$propfat
spdf$propinj = df1$propinj
Casualties = read.csv("storm_data_search_results_Missouri.csv", header = TRUE)
table(Casualties$CZ_NAME_STR, Casualties$DEATHS_DIRECT)##
## 0 1 2 3 4 5 7 8 10 11 14 37 158
## ADAIR CO. 0 0 1 0 0 0 0 0 0 0 0 0 0
## ANDREW CO. 6 1 0 0 0 0 0 0 0 0 0 0 0
## ATCHISON CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## BARRY CO. 6 4 0 0 0 0 0 0 0 0 0 0 0
## BARTON CO. 4 2 0 0 0 0 0 0 0 0 0 0 0
## BATES CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## BENTON CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## BOLLINGER CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## BOONE CO. 6 0 0 0 0 0 0 0 0 0 0 0 0
## BUCHANAN CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## BUTLER CO. 7 1 0 0 0 0 0 0 0 0 0 0 0
## CALDWELL CO. 2 1 0 0 0 0 0 0 0 0 0 0 0
## CALLAWAY CO. 2 1 0 0 0 0 0 0 0 0 0 0 0
## CAMDEN CO. 1 0 0 0 1 0 0 0 0 0 0 0 0
## CAPE GIRARDEAU CO. 6 1 0 0 0 0 0 0 0 0 0 0 0
## CARROLL CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## CARTER CO. 3 1 0 0 0 0 1 0 0 0 0 0 0
## CASS CO. 3 1 1 0 0 0 0 0 0 0 0 0 0
## CEDAR CO. 1 2 0 1 0 0 0 0 0 0 0 0 0
## CHARITON CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## CHRISTIAN CO. 6 1 0 0 0 0 0 0 0 0 0 0 0
## CLARK CO. 1 0 1 0 0 0 0 0 0 0 0 0 0
## CLAY CO. 7 0 0 0 0 0 0 0 0 0 0 0 0
## CLINTON CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## CRAWFORD CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## DADE CO. 6 0 0 0 0 0 0 0 0 0 0 0 0
## DALLAS CO. 2 1 1 0 0 0 0 0 0 0 0 0 0
## DAVIESS CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## DE KALB CO. 2 0 0 1 0 0 0 0 0 0 0 0 0
## DENT CO. 2 0 1 0 0 0 0 0 0 0 0 0 0
## DOUGLAS CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## DUNKLIN CO. 9 0 1 0 0 0 0 0 0 0 0 0 0
## FRANKLIN CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## GASCONADE CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## GENTRY CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## GREENE CO. 8 4 1 0 0 0 0 0 0 0 0 0 0
## GRUNDY CO. 4 0 0 0 0 0 0 0 0 0 0 0 0
## HARRISON CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## HENRY CO. 2 1 0 0 0 0 0 0 0 0 0 0 0
## HICKORY CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## HOWARD CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## HOWELL CO. 14 2 1 0 0 0 0 0 0 0 0 0 0
## IRON CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## JACKSON CO. 5 0 0 0 0 0 0 0 0 0 0 1 0
## JASPER CO. 10 3 1 0 0 0 0 0 0 0 0 0 1
## JEFFERSON CO. 0 1 0 0 0 0 0 0 0 0 0 0 0
## JOHNSON CO. 4 1 0 0 0 0 0 0 0 0 0 0 0
## LACLEDE CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## LAFAYETTE CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## LAWRENCE CO. 2 1 1 0 0 1 0 0 0 0 0 0 0
## LEWIS CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## LINCOLN CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## LINN CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## LIVINGSTON CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## MACON CO. 3 2 0 0 0 0 0 0 0 0 0 0 0
## MADISON CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## MARIES CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## MARION CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## MCDONALD CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## MERCER CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## MILLER CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## MISSISSIPPI CO. 8 1 0 0 0 0 0 0 0 0 0 0 0
## MONITEAU CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## MONROE CO. 2 0 1 0 0 0 0 0 0 0 0 0 0
## MONTGOMERY CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## MORGAN CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## NEW MADRID CO. 4 0 0 0 0 0 0 0 0 0 0 0 0
## NEWTON CO. 7 0 0 1 0 0 0 0 0 0 1 0 0
## NODAWAY CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## OREGON CO. 7 0 0 0 0 0 0 0 0 0 0 0 0
## OSAGE CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## OZARK CO. 6 0 0 0 0 0 0 0 0 0 0 0 0
## PEMISCOT CO. 2 0 1 1 0 0 0 0 0 0 0 0 0
## PERRY CO. 8 0 1 0 0 0 0 0 0 0 0 0 0
## PETTIS CO. 6 2 0 0 0 0 0 0 0 0 0 0 0
## PHELPS CO. 3 0 1 0 0 0 0 0 0 0 0 0 0
## PLATTE CO. 4 0 0 0 0 0 0 0 0 0 0 0 0
## POLK CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## PULASKI CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## PUTNAM CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## RANDOLPH CO. 2 0 0 0 1 0 0 0 0 0 0 0 0
## RAY CO. 3 0 1 0 0 0 0 0 0 0 0 0 0
## REYNOLDS CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## RIPLEY CO. 2 1 0 0 0 0 0 0 0 0 0 0 0
## SALINE CO. 4 0 0 0 0 0 0 0 0 0 0 0 0
## SCHUYLER CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## SCOTLAND CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## SCOTT CO. 7 2 0 2 0 0 0 0 0 0 0 0 0
## SHANNON CO. 0 1 0 1 0 0 0 0 0 0 0 0 0
## SHELBY CO. 1 0 0 0 0 0 0 0 0 0 0 0 0
## ST. CHARLES CO. 13 0 0 0 0 0 0 0 0 0 0 0 0
## ST. CLAIR CO. 3 0 0 0 0 0 0 0 0 0 0 0 0
## ST. FRANCOIS CO. 3 1 0 0 1 0 0 1 0 0 0 0 0
## ST. LOUIS (C) CO. 2 0 0 0 0 0 0 0 0 1 0 0 0
## ST. LOUIS CO. 13 2 0 1 0 0 0 0 1 0 0 0 0
## STE. GENEVIEVE CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## STODDARD CO. 3 1 0 0 0 0 0 0 0 0 0 0 0
## STONE CO. 4 0 0 0 0 0 0 0 0 0 0 0 0
## SULLIVAN CO. 0 1 0 0 0 0 0 0 0 0 0 0 0
## TANEY CO. 2 0 0 0 0 0 0 0 0 0 0 0 0
## TEXAS CO. 5 0 0 0 0 0 0 0 0 0 0 0 0
## VERNON CO. 11 0 0 0 0 0 0 0 0 0 0 0 0
## WARREN CO. 1 1 0 0 0 0 0 0 0 0 0 0 0
## WASHINGTON CO. 5 1 1 1 0 0 0 0 0 0 0 0 0
## WAYNE CO. 4 0 0 0 0 0 0 0 0 0 0 0 0
## WEBSTER CO. 7 0 1 0 0 0 0 0 0 0 0 0 0
## WORTH CO. 5 0 1 0 0 0 0 0 0 0 0 0 0
## WRIGHT CO. 3 0 0 0 0 0 0 0 0 0 0 0 0AC = Casualties %>%
group_by(CZ_NAME_STR) %>%
summarize(nD = sum(DEATHS_DIRECT),
nI = sum(INJURIES_DIRECT))
AC[1:10, ]; AC[11:20, ]; AC[21:30, ]; AC[31:40, ]; AC[41:50, ]; AC[51:60, ]## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 ADAIR CO. 2 6
## 2 ANDREW CO. 1 16
## 3 ATCHISON CO. 0 0
## 4 BARRY CO. 4 44
## 5 BARTON CO. 2 18
## 6 BATES CO. 0 3
## 7 BENTON CO. 0 5
## 8 BOLLINGER CO. 1 23
## 9 BOONE CO. 0 26
## 10 BUCHANAN CO. 0 21## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 BUTLER CO. 1 45
## 2 CALDWELL CO. 1 4
## 3 CALLAWAY CO. 1 5
## 4 CAMDEN CO. 4 30
## 5 CAPE GIRARDEAU CO. 1 8
## 6 CARROLL CO. 0 1
## 7 CARTER CO. 8 85
## 8 CASS CO. 3 26
## 9 CEDAR CO. 5 42
## 10 CHARITON CO. 0 8## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 CHRISTIAN CO. 1 13
## 2 CLARK CO. 2 23
## 3 CLAY CO. 0 30
## 4 CLINTON CO. 0 0
## 5 CRAWFORD CO. 0 5
## 6 DADE CO. 0 6
## 7 DALLAS CO. 3 27
## 8 DAVIESS CO. 0 2
## 9 DE KALB CO. 3 13
## 10 DENT CO. 2 4## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 DOUGLAS CO. 0 2
## 2 DUNKLIN CO. 2 31
## 3 FRANKLIN CO. 0 18
## 4 GASCONADE CO. 0 10
## 5 GENTRY CO. 0 1
## 6 GREENE CO. 6 136
## 7 GRUNDY CO. 0 9
## 8 HARRISON CO. 1 5
## 9 HENRY CO. 1 20
## 10 HICKORY CO. 0 21## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 HOWARD CO. 0 1
## 2 HOWELL CO. 4 63
## 3 IRON CO. 0 9
## 4 JACKSON CO. 37 182
## 5 JASPER CO. 163 1257
## 6 JEFFERSON CO. 1 0
## 7 JOHNSON CO. 1 12
## 8 LACLEDE CO. 0 23
## 9 LAFAYETTE CO. 1 5
## 10 LAWRENCE CO. 8 54## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 LEWIS CO. 0 11
## 2 LINCOLN CO. 0 8
## 3 LINN CO. 0 16
## 4 LIVINGSTON CO. 0 0
## 5 MACON CO. 2 12
## 6 MADISON CO. 0 4
## 7 MARIES CO. 0 1
## 8 MARION CO. 0 1
## 9 MCDONALD CO. 0 8
## 10 MERCER CO. 0 1NMc = as.character(AC$CZ_NAME_STR)
NMc = unlist(strsplit(NMc, split = " CO."))
AC$NAME = NMc
#AC = AC[order(NMc), ]
nm = data.frame(NAME = toupper(ST.sp$NAME[order(ST.sp$NAME)]), stringsAsFactors=FALSE)
AC2 = left_join(nm, AC[, -1])## Joining, by = "NAME"AC2[is.na(AC2)] = 0
IDs = as.integer(sapply(slot(spdf, "polygons"), function(x) slot(x, "ID")))
Names = toupper(ST.sp$NAME)
IDs = IDs[order(Names)]
df = data.frame(ID = ST.sp$GEOID[order(Names)],
Name = Names[order(Names)],
AC2,
spdf@data[order(Names), ])
cor(df$nD, df$nT, method = "spearman")## [1] 0.3189463cor(df$nI, df$nT, method = "spearman")## [1] 0.4001136cor(df$nD, df$propfat, method = "spearman")## [1] 0.8183287cor(df$nI, df$propinj, method = "spearman")## [1] 0.8725444sum(abs(df$nD - df$propfat))/length(df$nD)## [1] 3.479256sum(abs(df$nI - df$propinj))/length(df$nI)## [1] 28.17859Kentucky
FP = 21
AB = "KY"
ST.sp = US.sp[US.sp$STATEFP == FP, ]
nrow(ST.sp)## [1] 120county = paste(ST.sp$STATEFP, ST.sp$COUNTYFP, sep = "")
countiesun.sp = geometry(spChFIDs(ST.sp, county))
counties.sp = spTransform(countiesun.sp, CRS(proj4string(TornL)))
rgeos::gArea(counties.sp) / 10^6## [1] 103608.9tc = over(TornP, counties.sp)
TornP2 = subset(TornP, !is.na(tc) & 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] 1sum(dup)/dim(TornP2@data)[1] * 100## [1] 0.1048218TornP3 = subset(TornP2, !dup)
dim(TornP3@data)## [1] 953 22TornP3T = spTransform(TornP3, CRS(projection(r)))
mapdata = fortify(TornP3T[TornP3T$fat > 0, ])## Regions defined for each Polygonsct = over(counties.sp, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
ID = sapply(slot(counties.sp, "polygons"), function(x) slot(x, "ID"))
df = data.frame(nT)
rownames(df) = ID
spdf = SpatialPolygonsDataFrame(counties.sp, df)
sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)## [1] 1180## [1] 9.833333## [1] 1## [1] 38## [1] 7.082151## [1] 5.100698AreaTor = gArea(TornP3, byid = TRUE)
AreaCounty = gArea(counties.sp, byid = TRUE)
inter = gIntersects(counties.sp, TornP3, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
Areas[[i]] = gArea(gIntersection(counties.sp[ids[i, 2], ], TornP3[ids[i, 1], ]))
}
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "Area"
Areas.df$County = ids[, 2]
Areas.df$Torn = ids[, 1]
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, countiesun.sp)
PopT = projectRaster(Pop, crs = crs(counties.sp))
pop = raster::extract(PopT, counties.sp, fun = mean, na.rm = TRUE)[, 1]
Areas.df$pop = pop[ids[, 2]]
Areas.df$AreaCounty = AreaCounty[ids[, 2]]
Areas.df$Tfat = TornP3$fat[ids[, 1]]
Areas.df$Tinj = TornP3$inj[ids[, 1]]
Areas.df$TArea = AreaTor[ids[, 1]]
Areas.df$weight = Areas.df$pop * Areas.df$Area
df = Areas.df %>%
group_by(Torn) %>%
summarize(tweight = sum(weight),
propin = sum(Area)/TArea[1])
Areas.df$tweight = (df$tweight/df$propin)[ids[, 1]]
Areas.df$propfat = Areas.df$Tfat * Areas.df$weight/Areas.df$tweight
Areas.df$propinj = Areas.df$Tinj * Areas.df$weight/Areas.df$tweight
# number of fatalities attributed to a given area or county
df1 = Areas.df %>%
group_by(County) %>%
summarize(propfat = sum(propfat),
propinj = sum(propinj),
propfatratio = sum(propfat)/(pop[1] * AreaCounty[1]),
propinjratio = sum(propinj)/(pop[1] * AreaCounty[1]))
sum(TornP3$inj)## [1] 3592sum(df1$propinj)## [1] 2917.401spdf$propfat = df1$propfat
spdf$propinj = df1$propinj
Casualties = read.csv("storm_data_search_results_Kentucky.csv", header = TRUE)
table(Casualties$CZ_NAME_STR, Casualties$DEATHS_DIRECT)##
## 0 1 2 3 4 6 7 8 31
## ADAIR CO. 3 0 0 0 0 1 0 0 0
## ALLEN CO. 1 0 0 0 1 0 0 0 0
## ANDERSON CO. 6 0 0 0 0 0 0 0 0
## BALLARD CO. 4 0 0 0 0 0 0 0 0
## BARREN CO. 7 1 1 0 0 0 0 0 0
## BATH CO. 1 0 0 0 0 0 0 0 0
## BELL CO. 0 1 0 0 0 0 0 0 0
## BOONE CO. 12 0 0 0 0 0 0 0 0
## BOURBON CO. 1 0 0 0 0 0 0 0 0
## BOYD CO. 3 0 0 0 0 0 0 0 0
## BOYLE CO. 4 1 0 0 0 0 0 0 0
## BRACKEN CO. 2 1 0 0 0 0 0 0 0
## BREATHITT CO. 0 0 1 0 0 0 0 0 0
## BRECKINRIDGE CO. 3 1 0 0 0 0 0 0 0
## BULLITT CO. 3 0 0 0 0 0 0 0 0
## BUTLER CO. 3 1 0 0 0 0 0 0 0
## CALDWELL CO. 4 0 0 0 0 0 0 0 0
## CALLOWAY CO. 10 1 0 0 0 0 0 0 0
## CAMPBELL CO. 2 0 0 0 0 0 0 0 0
## CARLISLE CO. 1 0 0 0 0 0 0 0 0
## CARROLL CO. 5 1 0 0 0 0 0 0 0
## CARTER CO. 1 0 0 0 0 0 0 0 0
## CASEY CO. 1 0 0 0 0 0 0 0 0
## CHRISTIAN CO. 12 0 0 0 0 0 0 0 0
## CLARK CO. 4 0 0 0 0 0 0 0 0
## CLAY CO. 4 0 0 0 0 0 0 0 0
## CLINTON CO. 2 1 0 0 0 0 0 1 0
## CRITTENDEN CO. 5 0 0 0 0 0 0 0 0
## CUMBERLAND CO. 3 0 0 0 0 0 0 0 0
## DAVIESS CO. 10 0 0 0 0 0 0 0 0
## EDMONSON CO. 2 0 0 0 0 0 0 0 0
## ELLIOTT CO. 1 0 0 0 0 0 0 0 0
## ESTILL CO. 2 0 0 0 0 0 0 0 0
## FAYETTE CO. 7 0 0 0 0 0 0 0 0
## FLEMING CO. 5 0 0 0 0 0 0 0 0
## FLOYD CO. 1 0 0 0 0 0 0 0 0
## FRANKLIN CO. 4 0 0 0 1 0 0 0 0
## FULTON CO. 3 0 0 0 0 0 0 0 0
## GALLATIN CO. 1 0 0 0 0 0 0 0 0
## GARRARD CO. 3 0 0 0 0 0 0 0 0
## GRANT CO. 0 1 0 0 0 0 0 0 0
## GRAVES CO. 5 0 0 0 0 0 0 0 0
## GRAYSON CO. 7 0 0 1 0 0 0 0 0
## GREEN CO. 3 0 0 0 0 0 0 0 0
## GREENUP CO. 2 0 0 0 0 0 0 0 0
## HANCOCK CO. 3 0 0 0 0 0 0 0 0
## HARDIN CO. 7 0 1 0 0 0 0 0 0
## HARRISON CO. 3 0 0 0 0 0 0 0 0
## HART CO. 4 1 0 0 0 0 0 0 0
## HENDERSON CO. 13 0 0 0 0 0 0 0 0
## HENRY CO. 6 1 0 0 0 0 0 0 0
## HICKMAN CO. 2 0 0 0 0 0 0 0 0
## HOPKINS CO. 7 0 0 0 0 0 0 0 0
## JACKSON CO. 1 0 0 0 0 0 0 0 0
## JEFFERSON CO. 6 0 0 1 0 0 0 0 0
## JESSAMINE CO. 2 0 0 0 0 0 0 0 0
## JOHNSON CO. 1 0 1 0 0 0 0 0 0
## KENTON CO. 6 0 0 0 1 0 0 0 0
## LARUE CO. 3 0 0 0 0 0 0 0 0
## LAUREL CO. 10 0 0 0 0 1 0 0 0
## LAWRENCE CO. 1 0 1 0 0 0 0 0 0
## LEE CO. 1 1 0 0 0 0 0 0 0
## LESLIE CO. 1 0 0 0 0 0 0 0 0
## LETCHER CO. 2 0 0 0 0 0 0 0 0
## LEWIS CO. 4 0 0 0 0 0 0 0 0
## LINCOLN CO. 6 0 0 0 0 0 0 0 0
## LIVINGSTON CO. 2 0 0 0 0 0 0 0 0
## LOGAN CO. 8 1 0 0 0 0 0 0 0
## LYON CO. 2 0 0 0 0 0 0 0 0
## MADISON CO. 9 0 1 0 0 0 1 0 0
## MAGOFFIN CO. 1 0 0 0 0 0 0 0 0
## MARION CO. 3 0 0 0 0 0 0 0 0
## MARSHALL CO. 2 2 0 1 0 0 0 0 0
## MASON CO. 4 0 0 0 0 0 0 0 0
## MCCRACKEN CO. 10 0 0 0 0 0 0 0 0
## MCCREARY CO. 6 0 1 0 0 0 0 0 0
## MCLEAN CO. 1 0 0 0 0 0 0 0 0
## MEADE CO. 3 0 0 0 0 0 0 0 1
## MENIFEE CO. 1 0 1 0 0 0 0 0 0
## MERCER CO. 5 1 0 0 0 0 0 0 0
## METCALFE CO. 1 3 0 0 0 0 0 0 0
## MONROE CO. 1 0 0 0 0 0 0 0 0
## MONTGOMERY CO. 3 0 0 0 0 0 0 0 0
## MORGAN CO. 2 0 0 0 0 1 0 0 0
## MUHLENBERG CO. 8 0 0 1 0 0 0 0 0
## NELSON CO. 5 1 0 0 0 0 0 0 0
## NICHOLAS CO. 1 0 0 0 0 0 0 0 0
## OHIO CO. 5 0 0 0 0 0 0 0 0
## OLDHAM CO. 4 0 0 0 0 0 0 0 0
## OWEN CO. 2 0 0 0 0 0 0 0 0
## OWSLEY CO. 1 0 0 0 0 0 0 0 0
## PENDLETON CO. 3 0 0 0 1 0 0 0 0
## PERRY CO. 1 0 0 0 0 0 0 0 0
## POWELL CO. 4 0 0 0 0 0 0 0 0
## PULASKI CO. 9 0 0 0 0 1 0 0 0
## ROBERTSON CO. 1 0 0 0 0 0 0 0 0
## ROCKCASTLE CO. 2 1 0 0 0 0 0 0 0
## ROWAN CO. 1 0 0 0 0 0 0 0 0
## RUSSELL CO. 4 0 1 0 0 0 0 0 0
## SCOTT CO. 9 0 0 0 0 0 0 0 0
## SHELBY CO. 8 0 0 0 0 0 0 0 0
## SIMPSON CO. 7 1 0 0 0 0 0 0 0
## SPENCER CO. 4 0 0 0 0 0 0 0 0
## TAYLOR CO. 4 0 0 0 0 0 0 0 0
## TODD CO. 5 0 0 0 0 0 0 0 0
## TRIGG CO. 3 0 0 0 0 0 0 0 0
## TRIMBLE CO. 3 0 0 0 0 0 0 0 0
## UNION CO. 10 0 0 0 0 0 0 0 0
## WARREN CO. 10 0 1 0 0 0 0 0 0
## WASHINGTON CO. 2 0 0 0 0 0 0 0 0
## WAYNE CO. 4 0 2 0 0 0 0 0 0
## WEBSTER CO. 7 1 0 0 0 0 0 0 0
## WHITLEY CO. 3 0 0 0 0 0 0 0 0
## WOODFORD CO. 3 0 0 0 0 0 0 0 0AC = Casualties %>%
group_by(CZ_NAME_STR) %>%
summarize(nD = sum(DEATHS_DIRECT),
nI = sum(INJURIES_DIRECT))
AC[1:10, ]; AC[11:20, ]; AC[21:30, ]## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 ADAIR CO. 6 65
## 2 ALLEN CO. 4 11
## 3 ANDERSON CO. 0 11
## 4 BALLARD CO. 0 10
## 5 BARREN CO. 3 51
## 6 BATH CO. 0 0
## 7 BELL CO. 1 15
## 8 BOONE CO. 0 41
## 9 BOURBON CO. 0 1
## 10 BOYD CO. 0 0## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 BOYLE CO. 1 43
## 2 BRACKEN CO. 1 8
## 3 BREATHITT CO. 2 7
## 4 BRECKINRIDGE CO. 1 20
## 5 BULLITT CO. 0 24
## 6 BUTLER CO. 1 1
## 7 CALDWELL CO. 0 8
## 8 CALLOWAY CO. 1 35
## 9 CAMPBELL CO. 0 10
## 10 CARLISLE CO. 0 2## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 CARROLL CO. 1 15
## 2 CARTER CO. 0 0
## 3 CASEY CO. 0 0
## 4 CHRISTIAN CO. 0 60
## 5 CLARK CO. 0 13
## 6 CLAY CO. 0 0
## 7 CLINTON CO. 9 65
## 8 CRITTENDEN CO. 0 9
## 9 CUMBERLAND CO. 0 35
## 10 DAVIESS CO. 0 33NMc = as.character(AC$CZ_NAME_STR)
NMc = unlist(strsplit(NMc, split = " CO."))
AC$NAME = NMc
#AC = AC[order(NMc), ]
nm = data.frame(NAME = toupper(ST.sp$NAME[order(ST.sp$NAME)]), stringsAsFactors=FALSE)
AC2 = left_join(nm, AC[, -1])## Joining, by = "NAME"AC2[is.na(AC2)] = 0
IDs = as.integer(sapply(slot(spdf, "polygons"), function(x) slot(x, "ID")))
Names = toupper(ST.sp$NAME)
IDs = IDs[order(Names)]
df = data.frame(ID = ST.sp$GEOID[order(Names)],
Name = Names[order(Names)],
AC2,
spdf@data[order(Names), ])
cor(df$nD, df$nT, method = "spearman")## [1] 0.3358898cor(df$nI, df$nT, method = "spearman")## [1] 0.5613826cor(df$nD, df$propfat, method = "spearman")## [1] 0.7349479cor(df$nI, df$propinj, method = "spearman")## [1] 0.8300582sum(abs(df$nD - df$propfat))/length(df$nD)## [1] 0.8505557sum(abs(df$nI - df$propinj))/length(df$nI)## [1] 12.86175Oklahoma
FP = 40
AB = "OK"
ST.sp = US.sp[US.sp$STATEFP == FP, ]
nrow(ST.sp)## [1] 77county = paste(ST.sp$STATEFP, ST.sp$COUNTYFP, sep = "")
countiesun.sp = geometry(spChFIDs(ST.sp, county))
counties.sp = spTransform(countiesun.sp, CRS(proj4string(TornL)))
rgeos::gArea(counties.sp) / 10^6## [1] 179767.6tc = over(TornP, counties.sp)
TornP2 = subset(TornP, !is.na(tc) & 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] 29sum(dup)/dim(TornP2@data)[1] * 100## [1] 0.8189777TornP3 = subset(TornP2, !dup)
dim(TornP3@data)## [1] 3512 22TornP3T = spTransform(TornP3, CRS(projection(r)))
mapdata = fortify(TornP3T[TornP3T$fat > 0, ])## Regions defined for each Polygonsct = over(counties.sp, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
ID = sapply(slot(counties.sp, "polygons"), function(x) slot(x, "ID"))
df = data.frame(nT)
rownames(df) = ID
spdf = SpatialPolygonsDataFrame(counties.sp, df)
sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)## [1] 4009## [1] 52.06494## [1] 21## [1] 110## [1] 19.10383## [1] 7.009636AreaTor = gArea(TornP3, byid = TRUE)
AreaCounty = gArea(counties.sp, byid = TRUE)
inter = gIntersects(counties.sp, TornP3, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
Areas[[i]] = gArea(gIntersection(counties.sp[ids[i, 2], ], TornP3[ids[i, 1], ]))
}
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "Area"
Areas.df$County = ids[, 2]
Areas.df$Torn = ids[, 1]
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, countiesun.sp)
PopT = projectRaster(Pop, crs = crs(counties.sp))
pop = raster::extract(PopT, counties.sp, fun = mean, na.rm = TRUE)[, 1]
Areas.df$pop = pop[ids[, 2]]
Areas.df$AreaCounty = AreaCounty[ids[, 2]]
Areas.df$Tfat = TornP3$fat[ids[, 1]]
Areas.df$Tinj = TornP3$inj[ids[, 1]]
Areas.df$TArea = AreaTor[ids[, 1]]
Areas.df$weight = Areas.df$pop * Areas.df$Area
df = Areas.df %>%
group_by(Torn) %>%
summarize(tweight = sum(weight),
propin = sum(Area)/TArea[1])
Areas.df$tweight = (df$tweight/df$propin)[ids[, 1]]
Areas.df$propfat = Areas.df$Tfat * Areas.df$weight/Areas.df$tweight
Areas.df$propinj = Areas.df$Tinj * Areas.df$weight/Areas.df$tweight
# number of fatalities attributed to a given area or county
df1 = Areas.df %>%
group_by(County) %>%
summarize(propfat = sum(propfat),
propinj = sum(propinj),
propfatratio = sum(propfat)/(pop[1] * AreaCounty[1]),
propinjratio = sum(propinj)/(pop[1] * AreaCounty[1]))
sum(TornP3$inj)## [1] 7554sum(df1$propinj)## [1] 5371.122spdf$propfat = df1$propfat
spdf$propinj = df1$propinj
Casualties = read.csv("storm_data_search_results_Oklahoma.csv", header = TRUE)
table(Casualties$CZ_NAME_STR, Casualties$DEATHS_DIRECT)##
## 0 1 2 3 4 5 6 7 8 10 11 12 13 16 20 24
## ADAIR CO. 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## ALFALFA CO. 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## ATOKA CO. 3 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## BEAVER CO. 4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## BECKHAM CO. 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## BLAINE CO. 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## BRYAN CO. 8 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## CADDO CO. 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## CANADIAN CO. 3 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0
## CARTER CO. 10 0 0 2 0 0 0 0 1 0 0 0 0 0 0 0
## CHEROKEE CO. 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## CHOCTAW CO. 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## CIMARRON CO. 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## CLEVELAND CO. 9 1 0 2 0 0 0 0 0 0 1 0 0 0 0 1
## COAL CO. 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## COMANCHE CO. 4 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## COTTON CO. 5 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## CRAIG CO. 10 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## CREEK CO. 4 0 1 0 0 2 0 0 0 0 0 0 1 0 0 0
## CUSTER CO. 4 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## DELAWARE CO. 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## DEWEY CO. 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## ELLIS CO. 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## GARFIELD CO. 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## GARVIN CO. 8 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## GRADY CO. 11 2 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## GRANT CO. 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## GREER CO. 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## HARMON CO. 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## HARPER CO. 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## HASKELL CO. 2 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## HUGHES CO. 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## JACKSON CO. 8 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0
## JEFFERSON CO. 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## JOHNSTON CO. 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## KAY CO. 11 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## KINGFISHER CO. 9 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## KIOWA CO. 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## LATIMER CO. 4 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## LE FLORE CO. 8 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0
## LINCOLN CO. 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## LOGAN CO. 7 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## MAJOR CO. 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## MARSHALL CO. 3 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## MAYES CO. 10 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## MCCLAIN CO. 6 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## MCCURTAIN CO. 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## MCINTOSH CO. 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## MURRAY CO. 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## MUSKOGEE CO. 3 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## NOBLE CO. 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## NOWATA CO. 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## OKFUSKEE CO. 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## OKLAHOMA CO. 23 2 1 0 0 0 0 0 0 0 0 1 0 0 0 0
## OKMULGEE CO. 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0
## OSAGE CO. 6 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## OTTAWA CO. 3 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## PAWNEE CO. 5 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## PAYNE CO. 5 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## PITTSBURG CO. 8 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## PONTOTOC CO. 4 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## POTTAWATOMIE CO. 10 1 2 0 1 0 0 0 0 0 0 0 0 0 0 0
## PUSHMATAHA CO. 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## ROGER MILLS CO. 3 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## ROGERS CO. 10 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## SEMINOLE CO. 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## SEQUOYAH CO. 4 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0
## STEPHENS CO. 13 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## TEXAS CO. 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## TILLMAN CO. 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## TULSA CO. 16 1 1 0 0 1 0 1 0 0 0 0 0 0 0 0
## WAGONER CO. 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## WASHINGTON CO. 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## WASHITA CO. 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## WOODS CO. 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## WOODWARD CO. 2 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0AC = Casualties %>%
group_by(CZ_NAME_STR) %>%
summarize(nD = sum(DEATHS_DIRECT),
nI = sum(INJURIES_DIRECT))
AC[1:10, ]; AC[11:20, ]; AC[21:30, ]## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 ADAIR CO. 0 1
## 2 ALFALFA CO. 0 6
## 3 ATOKA CO. 2 52
## 4 BEAVER CO. 2 15
## 5 BECKHAM CO. 0 4
## 6 BLAINE CO. 2 10
## 7 BRYAN CO. 3 12
## 8 CADDO CO. 0 26
## 9 CANADIAN CO. 17 151
## 10 CARTER CO. 14 84## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 CHEROKEE CO. 2 5
## 2 CHOCTAW CO. 0 28
## 3 CIMARRON CO. 0 3
## 4 CLEVELAND CO. 42 642
## 5 COAL CO. 1 2
## 6 COMANCHE CO. 3 110
## 7 COTTON CO. 4 14
## 8 CRAIG CO. 1 28
## 9 CREEK CO. 25 282
## 10 CUSTER CO. 4 19## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 DELAWARE CO. 0 28
## 2 DEWEY CO. 0 6
## 3 ELLIS CO. 0 6
## 4 GARFIELD CO. 0 22
## 5 GARVIN CO. 5 24
## 6 GRADY CO. 14 108
## 7 GRANT CO. 0 8
## 8 GREER CO. 0 1
## 9 HARMON CO. 0 6
## 10 HARPER CO. 0 7NMc = as.character(AC$CZ_NAME_STR)
NMc = unlist(strsplit(NMc, split = " CO."))
AC$NAME = NMc
#AC = AC[order(NMc), ]
nm = data.frame(NAME = toupper(ST.sp$NAME[order(ST.sp$NAME)]), stringsAsFactors=FALSE)
AC2 = left_join(nm, AC[, -1])## Joining, by = "NAME"AC2[is.na(AC2)] = 0
IDs = as.integer(sapply(slot(spdf, "polygons"), function(x) slot(x, "ID")))
Names = toupper(ST.sp$NAME)
IDs = IDs[order(Names)]
df = data.frame(ID = ST.sp$GEOID[order(Names)],
Name = Names[order(Names)],
AC2,
spdf@data[order(Names), ])
cor(df$nD, df$nT, method = "spearman")## [1] 0.2539735cor(df$nI, df$nT, method = "spearman")## [1] 0.4804905cor(df$nD, df$propfat, method = "spearman")## [1] 0.7892082cor(df$nI, df$propinj, method = "spearman")## [1] 0.8627194sum(abs(df$nD - df$propfat))/length(df$nD)## [1] 1.786353sum(abs(df$nI - df$propinj))/length(df$nI)## [1] 21.69797Arkansas
FP = "05"
AB = "AR"
ST.sp = US.sp[US.sp$STATEFP == FP, ]
nrow(ST.sp)## [1] 75county = paste(ST.sp$STATEFP, ST.sp$COUNTYFP, sep = "")
countiesun.sp = geometry(spChFIDs(ST.sp, county))
counties.sp = spTransform(countiesun.sp, CRS(proj4string(TornL)))
rgeos::gArea(counties.sp) / 10^6## [1] 136986.2tc = over(TornP, counties.sp)
TornP2 = subset(TornP, !is.na(tc) & 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.238379TornP3 = subset(TornP2, !dup)
dim(TornP3@data)## [1] 1674 22TornP3T = spTransform(TornP3, CRS(projection(r)))
mapdata = fortify(TornP3T[TornP3T$fat > 0, ])## Regions defined for each Polygonsct = over(counties.sp, TornP3, returnList = TRUE)
nT = sapply(ct, function(x) nrow(x))
ID = sapply(slot(counties.sp, "polygons"), function(x) slot(x, "ID"))
df = data.frame(nT)
rownames(df) = ID
spdf = SpatialPolygonsDataFrame(counties.sp, df)
sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)## [1] 2070## [1] 27.6## [1] 8## [1] 82## [1] 13.73297## [1] 6.833137AreaTor = gArea(TornP3, byid = TRUE)
AreaCounty = gArea(counties.sp, byid = TRUE)
inter = gIntersects(counties.sp, TornP3, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
Areas[[i]] = gArea(gIntersection(counties.sp[ids[i, 2], ], TornP3[ids[i, 1], ]))
}
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "Area"
Areas.df$County = ids[, 2]
Areas.df$Torn = ids[, 1]
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, countiesun.sp)
PopT = projectRaster(Pop, crs = crs(counties.sp))
pop = raster::extract(PopT, counties.sp, fun = mean, na.rm = TRUE)[, 1]
Areas.df$pop = pop[ids[, 2]]
Areas.df$AreaCounty = AreaCounty[ids[, 2]]
Areas.df$Tfat = TornP3$fat[ids[, 1]]
Areas.df$Tinj = TornP3$inj[ids[, 1]]
Areas.df$TArea = AreaTor[ids[, 1]]
Areas.df$weight = Areas.df$pop * Areas.df$Area
df = Areas.df %>%
group_by(Torn) %>%
summarize(tweight = sum(weight),
propin = sum(Area)/TArea[1])
Areas.df$tweight = (df$tweight/df$propin)[ids[, 1]]
Areas.df$propfat = Areas.df$Tfat * Areas.df$weight/Areas.df$tweight
Areas.df$propinj = Areas.df$Tinj * Areas.df$weight/Areas.df$tweight
# number of fatalities attributed to a given area or county
df1 = Areas.df %>%
group_by(County) %>%
summarize(propfat = sum(propfat),
propinj = sum(propinj),
propfatratio = sum(propfat)/(pop[1] * AreaCounty[1]),
propinjratio = sum(propinj)/(pop[1] * AreaCounty[1]))
sum(TornP3$inj)## [1] 4735sum(df1$propinj)## [1] 4481.686spdf$propfat = df1$propfat
spdf$propinj = df1$propinj
Casualties = read.csv("storm_data_search_results_Arkansas.csv", header = TRUE)
table(Casualties$CZ_NAME_STR, Casualties$DEATHS_DIRECT)##
## 0 1 2 3 4 5 6 7 10 12 14 34
## ARKANSAS CO. 7 0 0 0 0 0 0 0 0 0 0 0
## ASHLEY CO. 5 0 0 2 0 0 0 0 0 0 0 0
## BAXTER CO. 4 1 0 1 0 0 0 0 0 0 0 0
## BENTON CO. 7 0 0 0 0 0 0 0 0 0 0 0
## BOONE CO. 4 1 0 0 0 0 0 0 0 0 0 0
## BRADLEY CO. 4 0 0 0 0 0 0 1 0 0 0 0
## CALHOUN CO. 2 0 0 0 0 0 0 0 0 0 0 0
## CARROLL CO. 3 0 0 0 0 0 0 0 0 0 0 0
## CHICOT CO. 5 1 0 0 0 0 0 0 0 0 0 0
## CLARK CO. 6 0 0 0 0 0 1 0 0 0 0 0
## CLAY CO. 0 1 0 0 0 0 0 0 0 0 0 0
## CLEBURNE CO. 7 2 1 0 0 0 0 0 0 0 0 0
## CLEVELAND CO. 1 0 0 0 0 0 0 0 0 0 0 0
## COLUMBIA CO. 6 1 1 0 0 0 0 0 0 0 0 0
## CONWAY CO. 6 4 2 0 0 0 0 0 0 0 0 0
## CRAIGHEAD CO. 6 0 0 1 0 0 0 0 0 0 0 1
## CRAWFORD CO. 6 0 0 0 0 0 0 0 0 0 0 0
## CRITTENDEN CO. 5 0 0 0 0 0 1 0 0 0 0 0
## CROSS CO. 3 1 0 0 0 0 0 0 0 0 0 0
## DALLAS CO. 0 1 0 0 0 0 0 0 0 0 0 0
## DESHA CO. 5 0 0 0 0 0 0 0 0 0 0 0
## DREW CO. 4 1 0 0 0 0 0 0 0 0 0 0
## FAULKNER CO. 15 3 1 0 1 0 1 0 0 1 0 0
## FRANKLIN CO. 4 1 1 0 0 0 0 0 0 0 0 0
## FULTON CO. 3 2 1 0 0 0 0 0 0 0 0 0
## GARLAND CO. 11 1 0 0 0 0 0 0 0 0 0 0
## GRANT CO. 3 1 0 0 0 0 0 0 0 0 0 0
## GREENE CO. 5 2 0 0 0 0 0 0 0 0 0 0
## HEMPSTEAD CO. 6 0 0 0 0 1 0 0 0 0 0 0
## HOT SPRING CO. 10 0 0 0 0 0 0 0 0 0 0 0
## HOWARD CO. 6 0 0 1 0 0 0 0 0 0 0 0
## INDEPENDENCE CO. 11 1 0 0 0 0 0 1 0 0 0 0
## IZARD CO. 4 0 2 0 0 0 0 0 0 0 0 0
## JACKSON CO. 9 2 1 0 0 0 0 0 0 0 0 0
## JEFFERSON CO. 5 1 0 0 0 0 0 0 0 0 0 0
## JOHNSON CO. 8 4 0 0 0 0 0 0 0 0 0 0
## LAWRENCE CO. 3 1 0 0 0 0 0 0 0 0 0 0
## LEE CO. 2 0 0 0 0 0 0 0 0 0 0 0
## LINCOLN CO. 3 0 0 0 0 1 0 0 0 0 0 0
## LITTLE RIVER CO. 3 2 0 0 0 0 0 0 0 0 0 0
## LOGAN CO. 8 1 0 0 0 0 0 0 0 0 0 0
## LONOKE CO. 12 0 1 0 0 1 0 0 0 0 0 0
## MADISON CO. 2 0 1 0 0 0 0 0 0 0 0 0
## MARION CO. 7 0 0 2 0 0 0 0 0 0 0 0
## MILLER CO. 3 0 0 0 0 0 0 0 0 0 0 0
## MISSISSIPPI CO. 10 0 1 0 0 0 0 0 0 0 0 0
## MONROE CO. 2 0 0 0 0 0 0 0 0 0 0 0
## MONTGOMERY CO. 2 0 0 0 0 0 0 0 0 0 0 0
## NEVADA CO. 4 0 0 0 0 0 0 0 0 0 0 0
## NEWTON CO. 2 0 0 0 0 0 0 0 0 0 0 0
## OUACHITA CO. 3 0 0 0 0 0 0 0 0 0 0 0
## PERRY CO. 4 0 0 0 0 0 0 0 0 0 0 0
## PHILLIPS CO. 5 0 0 0 0 0 0 0 0 0 0 0
## PIKE CO. 3 0 0 0 0 0 0 0 0 0 0 0
## POINSETT CO. 8 0 0 0 0 1 0 0 0 0 0 0
## POLK CO. 5 0 0 1 0 0 0 0 0 0 0 0
## POPE CO. 7 2 0 0 0 1 0 0 0 0 0 0
## PULASKI CO. 15 3 1 3 0 1 0 0 0 0 0 0
## RANDOLPH CO. 5 0 0 0 0 0 0 0 0 0 0 0
## SALINE CO. 9 2 0 0 0 0 0 0 1 0 0 0
## SCOTT CO. 2 0 0 0 0 0 0 0 0 0 0 0
## SEARCY CO. 4 0 0 0 0 0 0 0 0 0 0 0
## SEBASTIAN CO. 6 0 1 0 0 0 0 0 0 0 1 0
## SEVIER CO. 3 0 0 0 0 0 0 0 0 0 0 0
## SHARP CO. 6 0 0 0 0 0 0 0 0 0 0 0
## ST. FRANCIS CO. 3 2 0 0 1 0 0 0 0 0 0 0
## STONE CO. 2 1 0 0 0 1 0 0 0 0 0 0
## UNION CO. 6 1 1 0 0 0 0 0 0 0 0 0
## VAN BUREN CO. 7 2 0 2 0 0 0 0 0 0 0 0
## WASHINGTON CO. 6 1 0 0 1 0 0 0 0 0 0 0
## WHITE CO. 17 4 3 0 0 0 0 0 0 0 0 0
## WOODRUFF CO. 7 1 1 0 0 0 0 0 0 0 0 0
## YELL CO. 1 0 0 0 0 0 0 0 0 0 0 0AC = Casualties %>%
group_by(CZ_NAME_STR) %>%
summarize(nD = sum(DEATHS_DIRECT),
nI = sum(INJURIES_DIRECT))
AC[1:10, ]; AC[11:20, ]; AC[21:30, ]## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 ARKANSAS CO. 0 23
## 2 ASHLEY CO. 6 50
## 3 BAXTER CO. 4 70
## 4 BENTON CO. 0 24
## 5 BOONE CO. 1 66
## 6 BRADLEY CO. 7 59
## 7 CALHOUN CO. 0 1
## 8 CARROLL CO. 0 16
## 9 CHICOT CO. 1 43
## 10 CLARK CO. 6 118## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 CLAY CO. 1 1
## 2 CLEBURNE CO. 4 35
## 3 CLEVELAND CO. 0 0
## 4 COLUMBIA CO. 3 19
## 5 CONWAY CO. 8 52
## 6 CRAIGHEAD CO. 37 616
## 7 CRAWFORD CO. 0 82
## 8 CRITTENDEN CO. 6 126
## 9 CROSS CO. 1 19
## 10 DALLAS CO. 1 6## # A tibble: 10 x 3
## CZ_NAME_STR nD nI
## <fctr> <int> <int>
## 1 DESHA CO. 0 33
## 2 DREW CO. 1 4
## 3 FAULKNER CO. 27 516
## 4 FRANKLIN CO. 3 26
## 5 FULTON CO. 4 28
## 6 GARLAND CO. 1 56
## 7 GRANT CO. 1 3
## 8 GREENE CO. 2 62
## 9 HEMPSTEAD CO. 5 13
## 10 HOT SPRING CO. 0 62NMc = as.character(AC$CZ_NAME_STR)
NMc = unlist(strsplit(NMc, split = " CO."))
AC$NAME = NMc
#AC = AC[order(NMc), ]
nm = data.frame(NAME = toupper(ST.sp$NAME[order(ST.sp$NAME)]), stringsAsFactors=FALSE)
AC2 = left_join(nm, AC[, -1])## Joining, by = "NAME"AC2[is.na(AC2)] = 0
IDs = as.integer(sapply(slot(spdf, "polygons"), function(x) slot(x, "ID")))
Names = toupper(ST.sp$NAME)
IDs = IDs[order(Names)]
df = data.frame(ID = ST.sp$GEOID[order(Names)],
Name = Names[order(Names)],
AC2,
spdf@data[order(Names), ])
cor(df$nD, df$nT, method = "spearman")## [1] 0.5364371cor(df$nI, df$nT, method = "spearman")## [1] 0.5609784cor(df$nD, df$propfat, method = "spearman")## [1] 0.8844517cor(df$nI, df$propinj, method = "spearman")## [1] 0.9378985sum(abs(df$nD - df$propfat))/length(df$nD)## [1] 1.513527sum(abs(df$nI - df$propinj))/length(df$nI)## [1] 16.53203*** Table of county level data ***
Case studies (Xenia, OH tornado (1974) and Piedmont, AL tornado (1994))
Case Study 1
#Map = get_map(location = "Xenia Ohio", zoom = 9)
#str(Map)
xmn = -84.1
xmx = -83.5
ymn = 39.5
ymx = 40
res = .09
r = raster(xmn = xmn, xmx = xmx,
ymn = ymn, ymx = ymx,
resolution = res)
ext = as.vector(extent(r))
Xe = subset(TornP, yr == 1974 & mo == 4 & mag == 5 & st == "OH")
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, extent(r))
#pop = resample(aggregate(Pop,
# fact = c(nrow(Pop)/nrow(r), ncol(Pop)/ncol(r)),
# fun = mean), r)
pop = Pop
cellStats(pop, stat = "mean"); cellStats(pop, stat = "sd")## [1] 93.46892## [1] 169.7215XeT = spTransform(Xe, CRS(projection(r)))
sp = as(pop, 'SpatialPolygons')
spT = spTransform(sp, CRS(proj4string(TornL)))
AT = gArea(Xe, byid = TRUE)
AC = gArea(spT, byid = TRUE)
inter = gIntersects(spT, Xe, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
A = list()
for(i in 1:nrow(ids)){
A[[i]] = gArea(gIntersection(spT[ids[i, 2], ], Xe[ids[i, 1], ]))
}
A.df = data.frame(matrix(unlist(A), ncol = 1, byrow = TRUE))
names(A.df) = "Area"
A.df$Cell = ids[, 2]
A.df$Torn = ids[, 1]
A.df$pop = values(pop)[ids[, 2]]
A.df$AC = AC[ids[, 2]]
A.df$Tfat = Xe$fat[ids[, 1]]
A.df$Tinj = Xe$inj[ids[, 1]]
A.df$TArea = AT[ids[, 1]]
A.df$weight = A.df$pop * A.df$Area
A.df$probfat = A.df$weight/sum(A.df$weight) * 100
A.df$propinj = A.df$Tinj * A.df$weight/sum(A.df$weight)
pfat = pop
values(pfat) = NA
values(pfat)[A.df$Cell] = A.df$probfat
Map = get_map(location = c(-((xmx - xmn)/2 - xmx) , (ymx - ymn)/2 + ymn),
source = "stamen",
zoom = 10,
color = "bw",
maptype = "toner")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.75,-83.8&zoom=10&size=640x640&scale=2&maptype=terrain&sensor=false## Map from URL : http://tile.stamen.com/toner/10/272/387.png## Map from URL : http://tile.stamen.com/toner/10/273/387.png## Map from URL : http://tile.stamen.com/toner/10/274/387.png## Map from URL : http://tile.stamen.com/toner/10/272/388.png## Map from URL : http://tile.stamen.com/toner/10/273/388.png## Map from URL : http://tile.stamen.com/toner/10/274/388.png## Map from URL : http://tile.stamen.com/toner/10/272/389.png## Map from URL : http://tile.stamen.com/toner/10/273/389.png## Map from URL : http://tile.stamen.com/toner/10/274/389.png#ggmap(Map, dev = "extent")
cr = brewer.pal(8, "Reds")
pfatP = rasterToPoints(pfat)
pfatP.df = data.frame(pfatP)
names(pfatP.df)[3] = "Probability"
XeT = spTransform(Xe, CRS(projection(r)))
path = fortify(XeT)## Regions defined for each Polygonstrack = fortify(XeT)## Regions defined for each Polygonsxs = track$long[which.min(track$long)]
xe = track$long[which.max(track$long)]
ys = track$lat[which.min(track$lat)]
ye = track$lat[which.max(track$lat)]
deaths = data.frame(x = c(-83.929, -83.973, -83.878), y = c(39.685, 39.672, 39.716), deaths = c(25, 9, 2))
CS1 = ggmap(Map, extent="device") +
geom_tile(aes(x, y, fill = Probability), data = pfatP.df, alpha = .75) +
geom_segment(aes(x = xs, xend = xe, y = ys, yend = ye),
color = "white", size = 3, arrow = arrow(length = unit(.03, "npc"))) +
geom_segment(aes(x = xs, xend = xe, y = ys, yend = ye),
color = "black", size = .5, arrow = arrow(length = unit(.03, "npc"))) +
geom_point(aes(x = x, y = y, color = factor(deaths)), data = deaths, size = 4) +
# geom_polygon(aes(x = long, y = lat), data = path) +
scale_fill_gradient(low = cr[1], high = cr[8], guide = guide_colorbar(title = "Probability of at Least One Death (%)")) +
# scale_color_discrete(guide = guide_legend(title = "Number of Deaths")) +
scale_color_manual(guide = guide_legend(title = "Number of Deaths", order = 1), values = c("gray", "darkgray", "black")) +
coord_equal() +
ylab(expression(paste("Latitude [", {}^o,"N]"))) +
xlab(expression(paste("Longitude [", {}^o,"E]"))) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
legend.position = "bottom") +
ggtitle("A") +
theme(plot.title=element_text(hjust=0))## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## insteadCase Study 2
#Map = get_map(location = "Jacksonville Alabama", zoom = 9)
#str(Map)
xmn = -86.6
xmx = -84.9
ymn = 33.3
ymx = 34.5
res = .09
r = raster(xmn = xmn, xmx = xmx,
ymn = ymn, ymx = ymx,
resolution = res)
ext = as.vector(extent(r))
STL = subset(TornP, yr == 1994 & mo == 3 & time == "10:55:00")
Pop = raster("usadens/usads00g/w001001.adf")
Pop = crop(Pop, extent(r))
#pop = resample(aggregate(Pop,
# fact = c(nrow(Pop)/nrow(r), ncol(Pop)/ncol(r)),
# fun = mean), r)
pop = Pop
cellStats(pop, stat = "mean"); cellStats(pop, stat = "sd")## [1] 42.9535## [1] 50.74915STLT = spTransform(STL, CRS(projection(r)))
sp = as(pop, 'SpatialPolygons')
spT = spTransform(sp, CRS(proj4string(TornL)))
AT = gArea(STL, byid = TRUE)
AC = gArea(spT, byid = TRUE)
inter = gIntersects(spT, STL, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
A = list()
for(i in 1:nrow(ids)){
A[[i]] = gArea(gIntersection(spT[ids[i, 2], ], STL[ids[i, 1], ]))
}
A.df = data.frame(matrix(unlist(A), ncol = 1, byrow = TRUE))
names(A.df) = "Area"
A.df$Cell = ids[, 2]
A.df$Torn = ids[, 1]
A.df$pop = values(pop)[ids[, 2]]
A.df$AC = AC[ids[, 2]]
A.df$Tfat = STL$fat[ids[, 1]]
A.df$Tinj = STL$inj[ids[, 1]]
A.df$TArea = AT[ids[, 1]]
A.df$weight = A.df$pop * A.df$Area
A.df$probfat = A.df$weight/sum(A.df$weight) * 100
A.df$propinj = A.df$Tinj * A.df$weight/sum(A.df$weight)
pfat = pop
values(pfat) = NA
values(pfat)[A.df$Cell] = A.df$probfat
Map = get_map(location = c(-((xmx - xmn)/2 - xmx) , (ymx - ymn)/2 + ymn),
source = "stamen",
zoom = 10,
color = "bw",
maptype = "toner")## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=33.9,-85.75&zoom=10&size=640x640&scale=2&maptype=terrain&sensor=false## Map from URL : http://tile.stamen.com/toner/10/266/408.png## Map from URL : http://tile.stamen.com/toner/10/267/408.png## Map from URL : http://tile.stamen.com/toner/10/268/408.png## Map from URL : http://tile.stamen.com/toner/10/269/408.png## Map from URL : http://tile.stamen.com/toner/10/266/409.png## Map from URL : http://tile.stamen.com/toner/10/267/409.png## Map from URL : http://tile.stamen.com/toner/10/268/409.png## Map from URL : http://tile.stamen.com/toner/10/269/409.png## Map from URL : http://tile.stamen.com/toner/10/266/410.png## Map from URL : http://tile.stamen.com/toner/10/267/410.png## Map from URL : http://tile.stamen.com/toner/10/268/410.png## Map from URL : http://tile.stamen.com/toner/10/269/410.png#ggmap(Map, dev = "extent")
cr = brewer.pal(8, "Reds")
pfatP = rasterToPoints(pfat)
pfatP.df = data.frame(pfatP)
names(pfatP.df)[3] = "Probability"
STLT = spTransform(STL, CRS(projection(r)))
path = fortify(STLT)## Regions defined for each Polygonstrack = fortify(STLT)## Regions defined for each Polygonsxs1 = track$long[which.min(track$long)]
xe1 = track$long[which.max(track$long)]
ys1 = track$lat[which.min(track$lat)]
ye1 = track$lat[which.max(track$lat)]
deaths = data.frame(x = c(-86.043, -85.909, -85.617102), y = c(33.764, 33.820, 33.976087), deaths = c(1, 1, 20))
CS2 = ggmap(Map, extent="device") +
geom_tile(aes(x, y, fill = Probability), data = pfatP.df, alpha = .75) +
geom_segment(aes(x = xs1, xend = xe1, y = ys1, yend = ye1),
color = "white", size = 2, arrow = arrow(length = unit(.03, "npc"))) +
geom_segment(aes(x = xs1, xend = xe1, y = ys1, yend = ye1),
color = "black", size = .25, arrow = arrow(length = unit(.03, "npc"))) +
geom_point(aes(x = x, y = y, color = factor(deaths)), data = deaths, size = 4) +
# geom_polygon(aes(x = long, y = lat), data = path) +
scale_fill_gradient(low = cr[1], high = cr[8], guide = guide_colorbar(title = "Probability of at Least One Death (%)")) +
# scale_color_discrete(guide = guide_legend(title = "Number of Deaths")) +
scale_color_manual(guide = guide_legend(title = "Number of Deaths", order = 1), values = c("darkgray", "black")) +
coord_equal() +
ylab(expression(paste("Latitude [", {}^o,"N]"))) +
xlab(expression(paste("Longitude [", {}^o,"E]"))) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank(),
legend.position = "bottom") +
ggtitle("B") +
theme(plot.title=element_text(hjust=0))## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## insteadPlot Case Studies Together
mat = matrix(c(1, 2), nrow = 1, byrow = TRUE)
multiplot(CS1, CS2, layout = mat)
#Plot Removes geom_segment for one of the two tornadoesFigure 11: The probability of at least one death by raster cell in th (A) 1974 Xenia, OH and (B) 1994 Piedmont, AL tornadoes