Structural Property Losses From Tornadoes in Florida

James Elsner, Emily Ryan, Georgiana Strode

http://climateandlife.columbia.edu/2017/03/27/thunderstorms-pose-as-much-property-risk-as-hurricanes/ http://www.willisre.com/Media_Room/Press_Releases_(Browse_All)/2017/WillisRe_Impact_of_ENSO_on_US_Tornado_and_Hail_frequencies_Final.pdf

Get directory with files.

if(!file.exists("FloridaPropertyValues")){
download.file(url = "http://myweb.fsu.edu/jelsner/temp/FloridaPropertyValues.zip",
              destfile = "FloridaPropertyValues.zip")
  unzip("FloridaPropertyValues.zip")
}

Tornado reports per 10,000 square km at the state level over the past 30 years (since introduction of Doppler). All tornadoes and EF2+ tornadoes.

#setwd("~/FloridaPropertyValues")
suppressMessages(library(sp))
suppressMessages(library(dplyr))
suppressMessages(library(xtable))
suppressMessages(library(sf))
sfdf = sf::st_read(dsn = "torn", 
                   layer = "torn", 
                   stringsAsFactors = FALSE)

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

df = sfdf %>%
  filter(yr >= 1986) %>%
  group_by(st) %>%
  summarize(nT = n(),
            nST = sum(mag >= 2))
df2 = data.frame(st = state.abb, area = state.area * 2.58999, State = state.name)
df = merge(df, df2, by = 'st')
st_geometry(df) = NULL
df %>%
  mutate(Rate = round(nT/area * 10000),
         Rate2 = round(nST/area * 10000)) %>%
  arrange(desc(Rate)) %>%
  dplyr::select(State:Rate2) %>%
  slice(1:10) %>%
xtable(., digits = 0)

## % latex table generated in R 3.3.3 by xtable 1.8-2 package
## % Mon May  1 10:27:03 2017
## \begin{table}[ht]
## \centering
## \begin{tabular}{rlrr}
##   \hline
##  & State & Rate & Rate2 \\ 
##   \hline
## 1 & Florida & 117 & 5 \\ 
##   2 & Kansas & 117 & 10 \\ 
##   3 & Mississippi & 106 & 18 \\ 
##   4 & Illinois & 100 & 13 \\ 
##   5 & Iowa & 99 & 12 \\ 
##   6 & Oklahoma & 96 & 13 \\ 
##   7 & Maryland & 94 & 6 \\ 
##   8 & Alabama & 92 & 16 \\ 
##   9 & Louisiana & 88 & 10 \\ 
##   10 & Nebraska & 79 & 8 \\ 
##    \hline
## \end{tabular}
## \end{table}

Damage scale. Source: Storm Prediction Center

Rating	Wind Speed (m/s)	Description
EF0	29-37	Minor or no damage: Peels surface off some roofs; some damage to gutters or siding; branches broken off trees; shallow-rooted trees pushed over.
EF1	38-49	Moderate damage: Roofs severely stripped; mobile homes overturned or badly damaged; loss of exterior doors; windows and other glass broken.
EF2	50-60	Considerable damage: Roofs torn off well-constructed houses; foundations of frame homes shifted; mobile homes completely destroyed; large trees snapped or uprooted; light-object missiles generated; cars lifted off ground.
EF3	61-73	Severe damage: Entire stories of well-constructed houses destroyed; severe damage to large buildings such as shopping malls; trains overturned; trees debarked; heavy cars lifted off the ground and thrown; structures with weak foundations are badly damaged.

Other studies: http://climatecenter.fsu.edu/topics/tornadoes http://journals.ametsoc.org/doi/abs/10.1175/JCLI-D-11-00235.1

Florida property value data

http://www.iii.org/article/florida-hurricane-insurance-fact-file

Load the property value dataset.

suppressMessages(library(sf))
PV.sfdf = sf::st_read(dsn = "PropertyValues",
                      layer = "FL_Land_and_Property", 
                      stringsAsFactors = FALSE)

## Reading layer `FL_Land_and_Property' from data source `/Users/jameselsner/Dropbox/FloridaPropertyValues/PropertyValues' using driver `ESRI Shapefile'
## converted into: POLYGON
## Simple feature collection with 363506 features and 10 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -88.05448 ymin: 24.40058 xmax: -79.965 ymax: 31.63519
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

suppressMessages(library(dplyr))
class(PV.sfdf)

## [1] "sf"         "data.frame"

glimpse(PV.sfdf)

## Observations: 363,506
## Variables: 11
## $ LABEL      <chr> "NH 9100", "NH 9200", "NH 9300", "NH 9400", "NH 950...
## $ GOODGRIDID <dbl> 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 3...
## $ Shape_Leng <dbl> 0.03778702, 0.03778693, 0.03778675, 0.03778665, 0.0...
## $ Shape_Area <dbl> 8.906822e-05, 8.906777e-05, 8.906689e-05, 8.906643e...
## $ AV         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ LAND_VAL   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ STRUC_VAL  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ RES_AV     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ RES_LAND_V <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ RES_STRUC_ <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ geometry   <simple_feature> POLYGON((-80.0926045329728 ..., POLYGON(...

There are 363,506 cells. Each cell contains 10 variables.

Attribute	Description
LABEL	This label is second part of the U.S. National Grid identifier for 1-km individual grid cells. Sometimes the first part of the USNG identifier is omitted locally if the prefix is understood. The first part of the USNG identifier has been removed. It is either 16R for UTM zsone 16 or 17R for UTM zone 17.
GOODGRIDID	Positive real numbers that are automatically generated.
AV	Assessed value of properties, including land and structural value for all land use categories. All property values are aggregated to the USNG 1-km grid cell. Properties spanning multiple grids are disaggregated proportionally through interpolation.
LAND_VAL	Property land values are aggregated to 1-km grid cell. Land value is determined by the county property appraiser. This field contains all property parcels regardless of land use code.
STRUC_VAL	Structural value is a calculated field aggregated to 1-km grid cell. Structural value is calculated by subtracting land value from assessed value. If land value is unreasonable (e.g. zero or is greater than the assessed value), structural value is 90% of assessed value. All land use categories are considered.
RES_AV	Assessed values for only Residential categories are aggregated to a 1-km grid cell. Assessed value is determined by the county property appraiser. Only property parcels with a land use code of Residential are considered.
RES_LAND_VAL	Land values for residential properties are aggregated to a 1-km grid cell. Land value is determined by the county property appraiser. This field only aggregates property parcels with a Residential land use code.
RES_STRUC_VAL	Structural value is a calculated field for residential properties aggregated to 1-km grid cell. Structural value is calculated by subtracting land value from assessed value. If land value is unreasonable (e.g. zero or is greater than the assessed value), structural value is 90% of assessed value. Only land use categories of Residential are considered.

Interest is on structural values. Remove cells with structural values equal to zero.

df = PV.sfdf %>%
  dplyr::select(id = GOODGRIDID,
         sv = STRUC_VAL, 
         rsv = RES_STRUC_,
         geometry = geometry) %>%
  filter(sv > 0) %>%
  mutate(nrsv = sv - rsv) %>%
  sf::st_sf()
class(df)

## [1] "sf"         "data.frame"

dim(df)[1]

## [1] 91180

There are 91,180 1-km cells with non-zero structural property.

Create a choropleth map of the structural property values. The spplot method takes too long so here we use the plot method on the simple feature geometry.

suppressMessages(library(RColorBrewer))
cr = brewer.pal(8, "PuRd")
par(mar = c(0,0,1,0))
plot(df$geometry[df$sv > 1000000], col = cr[1], border = cr[1])
plot(df$geometry[df$sv > 10000000], col = cr[2], border = cr[2], add = TRUE)
plot(df$geometry[df$sv > 25000000], col = cr[3], border = cr[3], add = TRUE)
plot(df$geometry[df$sv > 50000000], col = cr[4], border = cr[4], add = TRUE)
plot(df$geometry[df$sv > 75000000], col = cr[5], border = cr[5], add = TRUE)
plot(df$geometry[df$sv > 100000000], col = cr[6], border = cr[6], add = TRUE)
plot(df$geometry[df$sv > 200000000], col = cr[7], border = cr[7], add = TRUE)
plot(df$geometry[df$sv > 250000000], col = cr[8], border = cr[8], add = TRUE)

suppressMessages(library(USAboundaries))
#suppressMessages(library(USAboundariesData)) #installed from source. 
FL = us_counties(states = "Florida", resolution = "high")
plot(FL, border = "grey", lwd = .5, add = TRUE)

leg.txt = c("$1-10M", "$10-25M", "$25-50M", "$50-75M", "$75-100M",
            "$100-200M", "$200-250M", "> $250M")
legend(x = -87, y = 28.7, leg.txt, pch = 15, 
       col = cr,
       title = "Structural Property Values\n (2014)",
       box.col = "transparent", box.lwd = .5)

#text(x = -85, y = 25, labels = "Source: Florida Resources and Environmental Analysis Center", cex = .7)

Summarize the statewide structural values.

df %>%
  summarize(n(), 
            sum(sv), min(sv), max(sv), mean(sv), median(sv),
            sum(rsv), min(rsv), max(rsv), mean(rsv), median(rsv),
            sum(nrsv), min(nrsv), max(nrsv), mean(nrsv), median(nrsv))

##     n()      sum(sv)      min(sv)    max(sv) mean(sv) median(sv)
## 1 91180 941592518065 2.237679e-05 8767929967 10326744   484880.8
##      sum(rsv) min(rsv)  max(rsv) mean(rsv) median(rsv)    sum(nrsv)
## 1 6.18948e+11        0 410798523   6788199    184876.5 322644500962
##   min(nrsv)  max(nrsv) mean(nrsv) median(nrsv)
## 1         0 8748864378    3538545     87495.02

The statewide total structural value of all property is $942 Billion and the total residential structural value is $619 Billion. At the per cell level, the statewide average total value is $10 Million with a median of $484,881. Since the mean >> median the values are strongly right skewed (long right tail).

The cell with the highest total structural value ($8.7 B) includes Allen D. Nease Senior High School in Ponte Vedra, in St. John’s County. The cell with the highest residential structural value ($410 M) includes the area around Lake Eola Park in downtown Orlando.

Plot the paths on a leaflet map.

suppressMessages(library(leaflet))
leaflet(df[df$sv == max(df$sv), ]) %>%
  addPolygons(stroke = FALSE,
              fillOpacity = .25,
              smoothFactor = 0) %>%
  addTiles()

leaflet(df[df$rsv == max(df$rsv), ]) %>%
  addPolygons(stroke = FALSE,
              fillOpacity = .25,
              smoothFactor = 0) %>%
  addTiles()

Distribution of property exposures by grid cell.

dfL = reshape2::melt(df, measure.vars = c("rsv", "nrsv"))
suppressMessages(library(ggplot2))
suppressMessages(library(scales))
suppressMessages(library(ggthemes))
ggplot(dfL, aes(log10(value), fill = variable)) +
  geom_histogram(breaks = seq(0, 10, .5), color = "white", position = "stack") +
  scale_x_continuous(breaks = seq(0, 8, 2), 
                     labels = c("$1", "$100", "$10,000", "$1,000,000", "$100,000,000")) +
  ylab("Number of One Kilometer Grid Cells") + xlab("") +
  theme_minimal() +
  theme(legend.title = element_blank(), legend.position = "bottom") +
  scale_fill_discrete(labels = c("Residential", "Non Residential")) +
  ggtitle("Florida Structural Values") 

## Warning: Removed 36630 rows containing non-finite values (stat_bin).

#  labs(caption = "Department of Revenue [2014]")

On average, residential property values exceed non-residential values. The most common residential exposures at the cell level are between $100,000 and $1,000,000 while the most common non-residential exposures are between $50,000 and $500,000.

sum(df$rsv > 50000000)/length(df$rsv) * 100

## [1] 3.290195

sum(df$nrsv > 50000000)/length(df$nrsv) * 100

## [1] 1.191051

The percentage of cells in the state with values exceeding $50 M is 3.3% for residential values and 1.2% for non-residential values.

Create a spatial polygons data frame of the cells having at least some structural property value. The CRS is geographic.

svPg = as(df, "Spatial")

Load and summarize the tornado data

Tornado data. Data from the Storm Prediction Center (SPC) http://www.spc.noaa.gov/gis/svrgis/. Track = center line (straight line in the tornado dataset). Path = width and length (area). Buffer the track to create a path.

download.file("http://www.spc.noaa.gov/gis/svrgis/zipped/tornado.zip",
              "tornado.zip", mode = "wb")
unzip("tornado.zip")

Create a point map for genesis locations with size, color, and alpha level a function of F/EF rating.

suppressMessages(library(sp))
sfdfL = sf::st_read(dsn = "torn", 
                   layer = "torn", 
                   stringsAsFactors = FALSE)

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

sfdfL = sfdfL %>%
  filter(st == "FL", yr >= 1950) %>%
  sf::st_sf()
spdfP = as.data.frame(sfdfL)
coordinates(spdfP) = ~ slon + slat
proj4string(spdfP) = "+init=epsg:4326"
sfdfP = sf::st_as_sf(spdfP)

## Warning in st_as_sf.Spatial(spdfP): column "geometry" will be overwritten
## by geometry column

suppressMessages(library(scales))
par(mar = c(0,0,1,0))
cr = c(alpha("black", .2), alpha("blue", .3), alpha("green", .4),
       alpha("orange", .5), col = alpha("red", .6))
plot(sfdfP$geometry[sfdfP$mag == 0], cex = .1, pch = 19, col = cr[1])
plot(sfdfP$geometry[sfdfP$mag == 1], cex = .25, pch = 19, col = cr[2], add = TRUE)
plot(sfdfP$geometry[sfdfP$mag == 2], cex = .75, pch = 19, col = cr[3], add = TRUE)
plot(sfdfP$geometry[sfdfP$mag == 3], cex = 1, pch = 19, col = cr[4], add = TRUE)
plot(sfdfP$geometry[sfdfP$mag == 4], cex = 1.25, pch = 19, col = cr[5], add = TRUE)

plot(FL, border = "grey", lwd = .5, add = TRUE)

leg.txt = c("F/EF0", "F/EF1", "F/EF2", "F/EF3", "F/EF4")
legend(x = -87, y = 28.7, leg.txt, pch = 19, 
       col = cr,
       pt.cex = c(.1, .25, .75, 1, 1.25),
       title = "Tornado Genesis\n (1950-2015)",
       box.col = "transparent", box.lwd = .5)

#text(x = -85, y = 25, labels = "Source: Storm Prediction Center", cex = .7)

There are 3,233 tornado reports with genesis locations attributed to the state in the 66-year period 1950–2015 for an annual average of 49 reports per year. Sixty-one percent of the tornadoes are rated F/EF0 and 28% are rated F/EF1. There are only 37 F/EF3 and two F/EF4 reports since 1950.

sfdfL %>%
  group_by(mag) %>%
  summarize(nT = n(),
            perc = nT/dim(sfdfL)[1] * 100)

## Simple feature collection with 5 features and 3 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: 815196.8 ymin: -1494784 xmax: 1616488 ymax: -728318.2
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
## # A tibble: 5 × 4
##     mag    nT        perc          geometry
##   <dbl> <int>       <dbl>  <simple_feature>
## 1     0  1959 60.59387566 <MULTILINESTR...>
## 2     1   917 28.36374884 <MULTILINESTR...>
## 3     2   318  9.83606557 <MULTILINESTR...>
## 4     3    37  1.14444788 <MULTILINESTR...>
## 5     4     2  0.06186205 <MULTILINESTR...>

Though our main concern is the potential for wind damage losses from tornadoes in Florida, we use event attributes from tornadoes that occurred in Florida and across southern portions of the neighboring states of Georgia and Alabama.

Load the data. Remove tornadoes with genesis locations outside Florida, Georgia, or Alabama. Convert to a SpatialPolygonsDataframe.

sfdfL2 = sf::st_read(dsn = "torn", 
                     layer = "torn", 
                     stringsAsFactors = FALSE)

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

sfdfL2 = sfdfL2 %>%
  filter(st == "FL" | st == "GA" | st == "AL", yr >= 1950) %>%
  sf::st_sf()
sfdfL2$wid = ifelse(sfdfL2$yr >= 1995, sfdfL2$wid * pi/4, sfdfL2$wid)
TornL = as(sfdfL2, "Spatial")

Annual average path width (meters).

sfdfL2 %>%
  group_by(yr) %>%
  summarize(mW = mean(wid) * .9144) %>%
ggplot(., aes(x = yr, y = mW)) +
  geom_point() +
  geom_smooth()

## `geom_smooth()` using method = 'loess'

Add a buffer to the track. Make the track a path. The SPC provides average path widths before 1995 and maximum widths thereafter. Since we are using rectangles representing the mean path we decrease the later widths by a factor of pi/4 to account for this change following Grieser and Terenzi (2016). See above.

TornP = rgeos::gBuffer(TornL, 
                       byid = TRUE, 
                       width = TornL$wid * .9144/2, 
                       capStyle = "ROUND")
TornP$AreaPath = rgeos::gArea(TornP, 
                              byid = TRUE)

Four percent of tornado paths have an area that exceeds one square km. The total area of all paths with genesis in Florida is 713 km\(^2\), which is .5% of the total area of the state.

TornPFL = TornP[TornP$st == "FL", ]
sum(TornPFL$AreaPath/10^6 > 1)/length(TornPFL$AreaPath) * 100

## [1] 3.742138

sum(TornPFL$AreaPath)

## [1] 713071272

sum(TornPFL$AreaPath)/10^6/df2$area[df2$st == "FL"] * 100

## [1] 0.4701471

Create a spatial polygons object of tornadoes with geographic CRS. Create a spatial polygons object of cells with a planar CRS.

TornPg = spTransform(TornP, CRS = CRS(proj4string(svPg)))
TornPFLg = spTransform(TornPFL, CRS = CRS(proj4string(svPg)))
svP = spTransform(svPg, CRS = CRS(proj4string(TornP)))

We now have two copies of the tornado paths and two copies of the structural property value grids. We also have a point file for Florida genesis tornadoes.

Plot the paths for tornadoes with genesis in Florida.

sfdfPFL = st_as_sf(TornPFLg)
plot(sfdfPFL$geometry)
plot(FL, border = "grey", lwd = .5, add = TRUE)

Plot the total path area by year.

as.data.frame(TornPFL) %>%
  filter(mag >= 0, yr >= 1988) %>%
  group_by(yr) %>%
  summarize(nT = n(),
            AreaPath = sum(AreaPath)) %>%
ggplot(., aes(yr, AreaPath/10^6)) +
  geom_line() +
  geom_point() +
  scale_y_log10(limits = c(.01, 1000), breaks = c(.01, .1, 1, 10, 100, 1000),
                labels = c(.01, .1, 1., 10, 100, 1000)) +
  scale_x_continuous(breaks = seq(1990, 2015, 5)) +
  ylab(expression(paste("Approximate Total Tornado Path Area [km", {}^2,"]"))) + 
  xlab("Year") +
  theme_minimal()

Tornadoes from tropical cyclones

The data were given by Todd Moore to Emily Ryan as a csv file. Read in the .csv file. Convert to spatial points using start long and lat. Add a CRS and convert to a simple features data frame.

FLtcTors = read.csv(file = "FLtctors.csv", 
                    header = TRUE, 
                    sep = ",")
coordinates(FLtcTors) = ~ StartLong + StartLat
proj4string(FLtcTors) = CRS("+init=epsg:4326")
FLtcTors.sfdf = st_as_sf(FLtcTors)
#plot(FLtcTors.sfdf$geometry)

Label the TC tornadoes

id = paste(FLtcTors$Year, FLtcTors$YearlyTracking_FromONETOR, sep = "")
id2 = paste(TornP$yr, TornP$om, sep = "")
TC = id2 %in% id
TornP$TC = TC

Overlay all tornadoes on the grid

Overlay an approximate tornado path on top of the structural value grid. Example: the 1956 Pembroke Pines to Ft Lauderdale F3 tornado.

strg = sfdfPFL %>%
  filter(mag >= 3)
strg[2, ]

## Simple feature collection with 1 feature and 23 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -80.2309 ymin: 25.99919 xmax: -80.1291 ymax: 26.07081
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##    om   yr mo dy       date     time tz st stf stn mag inj fat loss closs
## 2 139 1956  4 10 1956-04-10 17:48:00  3 FL  12   2   3  20   0    6     0
##   slat   slon  elat   elon len wid fc AreaPath
## 2   26 -80.23 26.07 -80.13 7.8 200  0  2386135
##                         geometry
## 2 POLYGON((-80.1305488154248 ...

leaflet(strg[2, ]) %>%
  addPolygons(stroke = FALSE,
              fillOpacity = .25,
              smoothFactor = 0) %>%
  addTiles()

Create figure showing a single tornado path on top of the property value grids.

cr = brewer.pal(8, "PuRd")
plot(x = c(-80.25, -80.13), y = c(25.97, 26.1), bty = "n", col = "transparent")
plot(strg$geometry[2], add = TRUE)
plot(df$geometry[df$sv > 1000000], col = adjustcolor(cr[1], alpha.f = .5), border = "white", add = TRUE)
plot(df$geometry[df$sv > 10000000], col = adjustcolor(cr[2], alpha.f = .5), border = "white", add = TRUE)
plot(df$geometry[df$sv > 25000000], col = adjustcolor(cr[3], alpha.f = .5), border = "white", add = TRUE)
plot(df$geometry[df$sv > 50000000], col = adjustcolor(cr[4], alpha.f = .5),  border = "white", add = TRUE)
plot(df$geometry[df$sv > 75000000], col = adjustcolor(cr[5], alpha.f = .5),  border = "white", add = TRUE)
plot(df$geometry[df$sv > 100000000], col = adjustcolor(cr[6], alpha.f = .5), border = "white", add = TRUE)
plot(df$geometry[df$sv > 200000000], col = adjustcolor(cr[7], alpha.f = .5), border = "white", add = TRUE)
plot(df$geometry[df$sv > 250000000], col = adjustcolor(cr[8], alpha.f = .5), border = "white", add = TRUE)
plot(strg$geometry[2], add = TRUE, lwd = 2)

Determine the number of tornado paths intersecting each cell.

ct = over(svP, TornP, returnList = TRUE)  
nT = sapply(ct, function(x) nrow(x))
ID = sapply(slot(svP, "polygons"), function(x) slot(x, "ID"))
df = data.frame(nT)
rownames(df) = ID
spdfnT = SpatialPolygonsDataFrame(svP, df)
sum(nT); mean(nT); min(nT); max(nT); sd(nT); var(nT)/mean(nT)

## [1] 10936

## [1] 0.1199386

## [1] 0

## [1] 20

## [1] 0.4353907

## [1] 1.580518

Map the number of tornadoes by grid cell. First convert sp to sf. Use a non-linear scale.

spdfnTg = spTransform(spdfnT, CRS("+init=epsg:4326"))
sfdfnTg = st_as_sf(spdfnTg)
suppressMessages(library(wesanderson))
cr = wes_palette(6, name = "Zissou", type = "continuous")
par(mar = c(0, 0, 1, 0))
plot(sfdfnTg$geometry[sfdfnTg$nT >= 1], col = cr[1], border = cr[1])
plot(sfdfnTg$geometry[sfdfnTg$nT >= 2], col = cr[2], border = cr[2], add = TRUE)
plot(sfdfnTg$geometry[sfdfnTg$nT >= 4], col = cr[3], border = cr[3], add = TRUE)
plot(sfdfnTg$geometry[sfdfnTg$nT >= 8], col = cr[4], border = cr[4], add = TRUE)
plot(sfdfnTg$geometry[sfdfnTg$nT >= 16], col = cr[5], border = cr[5], add = TRUE)
plot(sfdfnTg$geometry[sfdfnTg$nT >= 20], col = cr[6], border = cr[6], add = TRUE)
plot(FL, border = "grey", lwd = .5, add = TRUE)
leg.txt = c("1", "2-3", "4-7", "8-15", "16-19", ">19")
legend(x = -86.578, y = 29.174, leg.txt, pch = 15, 
       col = cr,
       title = "Number of Tornadoes\n (1950-2015)",
       box.col = "transparent", box.lwd = .5)

#text(x = -85, y = 25, labels = "Source: Storm Prediction Center", cex = .7)

Plot the cell with greatest tornado occurrence.

table(sfdfnTg$nT)

## 
##     0     1     2     3     4     5     6     7     8     9    10    11 
## 82020  8080   791   152    52    34    16     8     9     5     1     4 
##    12    14    16    17    20 
##     3     2     1     1     1

strg = sfdfnTg %>%
  filter(nT == 20)
strg[1, ]

## Simple feature collection with 1 feature and 1 field
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -82.67449 ymin: 27.76235 xmax: -82.66421 ymax: 27.7715
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   nT                       geometry
## 1 20 POLYGON((-82.6642112463519 ...

leaflet(strg[1, ]) %>%
  addPolygons(stroke = FALSE,
              fillOpacity = .25,
              smoothFactor = 0) %>%
  addTiles()

The total number of cells affected by tornadoes is 10,936 with 8,080 of the cells having been hit only once and 791 cells having been hit twice. The average number of tornadoes per cell with property is .12 with a minimum of zero and a maximum of 20 (once every 3 or 4 years on average) in the Palmetto Park area of St. Petersburg.

Determine the 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.

inter = rgeos::gIntersects(svP, TornP, byid = TRUE)
ids = which(inter, arr.ind = TRUE)
Areas = list()
for(i in 1:nrow(ids)){
    Areas[[i]] = rgeos::gArea(rgeos::gIntersection(svP[ids[i, 2], ], TornP[ids[i, 1], ]))
    }
Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
names(Areas.df) = "IntersectArea"
Areas.df$Cell = ids[, 2]
Areas.df$Torn = ids[, 1]

Use the proportion of path area within each cell to apportion the structural property values. This is different than the TorMC model of Strader et al. (2016).

AreaTor = rgeos::gArea(TornP, byid = TRUE)
AreaCell = rgeos::gArea(svP, byid = TRUE)
Areas.df$CellArea = AreaCell[Areas.df$Cell]
Areas.df$TornArea = AreaTor[Areas.df$Torn]
Areas.df$sv = svP$sv[Areas.df$Cell]
Areas.df$rsv = svP$rsv[Areas.df$Cell]
Areas.df$nrsv = svP$nrsv[Areas.df$Cell]
Areas.df$om = TornP$om[Areas.df$Torn]

df2 = Areas.df %>%
  mutate(perc = IntersectArea/CellArea,
         sv = sv * perc,
         rsv = rsv * perc,
         nrsv = nrsv * perc) %>%
  group_by(Torn) %>%
  summarize(sv = sum(sv),
            rsv = sum(rsv),
            nrsv = sum(nrsv),
            om = first(om))

Subset TornP by the tornadoes in df2 and save in TornP2 since not all tornadoes occurred over areas with structures (e.g., water) or with structural value grids (e.g., Georgia). Add the per tornado total and residential structural values as new columns in TornP2.

TornP2 = TornP[df2$Torn, ]
TornP2$sv = df2$sv
TornP2$rsv = df2$rsv
TornP2$nrsv = df2$nrsv

Compute the annual average exposure.

AvgSV = sum(TornP2$sv)/(max(TornP2$yr) - min(TornP2$yr) + 1)
AvgRSV = sum(TornP2$rsv)/(max(TornP2$yr) - min(TornP2$yr) + 1)
AvgNRSV = sum(TornP2$nrsv)/(max(TornP2$yr) - min(TornP2$yr) + 1)
AvgSV; AvgRSV; AvgNRSV

## [1] 156556581

## [1] 94417389

## [1] 62139192

Total number of tornadoes by month.

suppressMessages(library(scales))
bymonth = as.data.frame(TornP2) %>%
  filter(mag >= 0) %>%
  group_by(mo) %>%
  summarize(nT = n()) %>%
  mutate(moF = factor(month.abb, levels = month.abb))

ggplot(bymonth, aes(moF, nT)) +
  geom_histogram(stat = "identity") +
  scale_y_continuous(label = comma) +
  xlab("") +
  ylab("Total Number of Tornadoes")

## Warning: Ignoring unknown parameters: binwidth, bins, pad

Average property value exposure by month.

exbymonth = as.data.frame(TornP2) %>%
  filter(mag >= 0) %>%
  group_by(mo) %>%
  summarize(mean = mean(nrsv)) %>%
  mutate(moF = factor(month.abb, levels = month.abb))
  
ggplot(exbymonth, aes(moF, mean)) +
  geom_histogram(stat = "identity") +
  scale_y_continuous(labels = dollar) +
  xlab("") +
  ylab("Average Exposure")

## Warning: Ignoring unknown parameters: binwidth, bins, pad

Total exposure by EF rating.

as.data.frame(TornP2) %>%
  group_by(mag) %>%
  summarize(nT = n()/dim(TornP2)[1] * 100,
    totalrsv = sum(rsv)/sum(TornP2$rsv) * 100,
    totalnrsv = sum(nrsv)/sum(TornP2$nrsv) * 100)

## # A tibble: 5 × 4
##     mag          nT  totalrsv totalnrsv
##   <dbl>       <dbl>     <dbl>     <dbl>
## 1     0 58.86925795  8.996752  9.500809
## 2     1 29.54063604 26.648205 21.482755
## 3     2 10.28268551 30.931940 39.288039
## 4     3  1.27208481 25.187894 20.912429
## 5     4  0.03533569  8.235209  8.815967

While EF0 tornadoes account for 59% of all Florida tornadoes they represent only 9% of the property exposure. In contrast, tornadoes rated EF2 or higher account for less than 12% of all tornadoes but 64% of residential and 69% of non-residential property exposure.

Percentage of exposure for TC and non-TC tornadoes.

as.data.frame(TornP2) %>%
  filter(yr >= 1995) %>%
  group_by(TC) %>%
  summarize(totalrsv = sum(rsv),
            totalnrsv = sum(nrsv))

## # A tibble: 2 × 3
##      TC   totalrsv  totalnrsv
##   <lgl>      <dbl>      <dbl>
## 1 FALSE 1683673072 1054537758
## 2  TRUE  447295735  219839667

1683673072/(1683673072+447295735) * 100

## [1] 79.00975

1054537758/(1054537758+219839667) * 100

## [1] 82.74925

Model for losses using actual losses since 2007

Relationship between actual losses and loss exposure. Before 1996 this was the breakdown: 0 - Unknown, 1 - <$50, 2 - $50 - $500, 3 - $500 - $5,000, 4 - $5,000 - $50,000, 5 - $50,000 - $500,000, 6 - $500,000 - $5,000,000, 7 - $5,000,000 - $50,000,000, 8 - $ 50,000,000 - $500,000,000, 9 - >$500,000,000.

Losses are in Millions.

Loss versus exposure.

TornP2.df = as.data.frame(TornP2)
TornP2.df = TornP2.df %>%
  filter(sv > 0 & yr >= 2007 & loss > 0) %>%
  mutate(loss = loss * 10^6)

ggplot(TornP2.df, aes(x = log10(sv), y = log10(loss))) +
  geom_point() +
#  geom_abline(slope = 1, intercept = 0) +
  scale_x_continuous(breaks = seq(0, 8, 2), 
                     labels = c("$1", "$100", "$10,000", "$1,000,000",
                                "$100,000,000")) +
  scale_y_continuous(breaks = seq(0, 8, 2), 
                     labels = c("$1", "$100", "$10,000", "$1,000,000",
                                "$100,000,000")) +
  theme_minimal() +
  ylab("Losses") + xlab("Total Structural Property Value")

Models (rsv & nrsv) for actual losses based on 2014 exposures.

modelRSV = lm(log(loss) ~ log(rsv) + yr, data = TornP2.df[TornP2.df$rsv > 0, ])
summary(modelRSV)

## 
## Call:
## lm(formula = log(loss) ~ log(rsv) + yr, data = TornP2.df[TornP2.df$rsv > 
##     0, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9195 -1.1890 -0.2964  1.1362  4.9691 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 344.39323  113.88797   3.024  0.00294 ** 
## log(rsv)      0.33023    0.04761   6.936 1.12e-10 ***
## yr           -0.16776    0.05665  -2.962  0.00356 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.714 on 150 degrees of freedom
## Multiple R-squared:  0.2765, Adjusted R-squared:  0.2669 
## F-statistic: 28.66 on 2 and 150 DF,  p-value: 2.869e-11

wtRSV = coef(modelRSV)[2]
(2^coef(modelRSV)[2] - 1) * 100

## log(rsv) 
## 25.72148

modelNRSV = lm(log(loss) ~ log(nrsv) + yr, data = TornP2.df[TornP2.df$nrsv > 0, ])
summary(modelNRSV)

## 
## Call:
## lm(formula = log(loss) ~ log(nrsv) + yr, data = TornP2.df[TornP2.df$nrsv > 
##     0, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5874 -1.2369 -0.1844  1.3764  5.2840 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 296.14096  115.48080   2.564   0.0113 *  
## log(nrsv)     0.26851    0.04117   6.523 9.26e-10 ***
## yr           -0.14319    0.05745  -2.492   0.0137 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.75 on 155 degrees of freedom
## Multiple R-squared:  0.2372, Adjusted R-squared:  0.2274 
## F-statistic:  24.1 on 2 and 155 DF,  p-value: 7.701e-10

wtNRSV = coef(modelNRSV)[2]
(2^coef(modelNRSV)[2] - 1) * 100

## log(nrsv) 
##  20.45596

The coefficient is not 1 because: exposure is not the only component of a structure’s vulnerability, differences in structural durability, variation in the wind speed within the path,

The estimated exposure is multiplied by this percentage to obtain an estimate of loss.

AvgSL = AvgRSV * wtRSV + AvgNRSV * wtNRSV
AvgSL

## log(rsv) 
## 47864292

Determine the probable maximum loss using a Monte Carlo approach

We estimate the annual probable maximum loss (PML) as the value of the largest loss that could occur from tornadoes. This helps insurance underwriting decisions including the amount to cede to a reinsurer. The estimate is done using a Monte Carlo procedure. The approach is similar to Strader et al. (2016). Here we area weight the exposure by the percentage of grid cell affected by the path and we retain the spatial clustering by using a bootstrap resampling.

Convert to a psp object. This requires the maptools package.

suppressMessages(library(maptools))
suppressMessages(library(spatstat))
sfdf3 = sfdfL2 %>%
  mutate(EF = as.factor(mag),
         W = wid * .9144) %>%
  sf::st_sf()
TornL3 = as(sfdf3, "Spatial")
TornF.psp = as(TornL3["EF"], "psp")
TornW.psp = as(TornL3["W"], "psp")

Create owin objects from the state border of Florida and from the boundary box (edges).

suppressMessages(library(mapproj))
bndryS = map("state", 
             region = 'florida', 
             fill = TRUE,
             plot = FALSE)
bndrySP = map2SpatialPolygons(bndryS, 
                              IDs = bndryS$names,
                              proj4string = CRS("+proj=longlat +ellps=GRS80"))
bndrySP = spTransform(bndrySP, CRS(sf::st_crs(sfdf)$proj4string))
FL = as(bndrySP, 'owin')
BB = as.owin(edges(FL))
plot(BB)
plot(FL, add = TRUE)

Subset the psp object by the window.

TornF.psp = TornF.psp[BB]
TornW.psp = TornW.psp[BB]
plot(TornF.psp)

Create vectors of lengths, angles, and locations of tracks.

lengths = lengths.psp(TornF.psp)
angles = angles.psp(TornF.psp)
mpts = coords(midpoints.psp(TornF.psp))
widths = marks(TornW.psp)
maxEF = marks(TornF.psp)

Repeat M times. This takes a long time.

#M = 10000
rate = length(angles)/(2015 - 1950 + 1)
M = 1
rate = length(angles)
tsv = numeric()
trsv = numeric()
tnrsv = numeric()
for(i in 1:M){
  print(i)
  n = 1:length(angles)
  size = rpois(1, lambda = rate)
  idx = sample(n, size = size, replace = TRUE)
  aa = sample(angles, size = size, replace = TRUE)
  ww = widths[idx]
  ll = lengths[idx]
  ef = maxEF[idx]
  pt = mpts[sample(n, size, replace = TRUE), ]
  df = data.frame(xmid = pt[, 1], ymid = pt[, 2], 
                length = ll, angle = aa)
  random.psp = as.psp(df, window = BB)
  randomL = as(random.psp, "SpatialLines")
  rsldf = SpatialLinesDataFrame(randomL, 
                              data = data.frame(ww, ef), 
                              match.ID = FALSE)
  proj4string(rsldf) = CRS(proj4string(TornL))
  rspdf = rgeos::gBuffer(rsldf, 
                       byid = TRUE, 
                       width = rsldf$ww / 2, 
                       capStyle = "ROUND")
  inter = rgeos::gIntersects(svP, rspdf, byid = TRUE)
  ids = which(inter, arr.ind = TRUE)
  Areas = list()
    for(i in 1:nrow(ids)){
      Areas[[i]] = rgeos::gArea(rgeos::gIntersection(svP[ids[i, 2], ], rspdf[ids[i, 1], ]))
      }
  Areas.df = data.frame(matrix(unlist(Areas), ncol = 1, byrow = TRUE))
  names(Areas.df) = "IntersectArea"
  Areas.df$Cell = ids[, 2]
  Areas.df$Torn = ids[, 1]
  AreaTor = rgeos::gArea(rspdf, byid = TRUE)
  AreaCell = rgeos::gArea(svP, byid = TRUE)
  Areas.df$CellArea = AreaCell[Areas.df$Cell]
  Areas.df$TornArea = AreaTor[Areas.df$Torn]
  Areas.df$sv = svP$sv[Areas.df$Cell]
  Areas.df$rsv = svP$rsv[Areas.df$Cell]
  Areas.df$nrsv = svP$nrsv[Areas.df$Cell]
  Areas.df = Areas.df %>%
    mutate(perc = IntersectArea/CellArea,
           tsv = sv * perc,
           trsv = rsv * perc,
           tnrsv = nrsv * perc)
tsv = c(tsv, sum(Areas.df$tsv))
trsv = c(trsv, sum(Areas.df$trsv))
tnrsv = c(tnrsv, sum(Areas.df$tnrsv))
}

## [1] 1

Plot the paths for fake tornadoes in Florida. Use only when M = 1 above.

rsfdf = st_as_sf(rspdf)
plot(rsfdf$geometry)
plot(FL, border = "grey", lwd = .5, add = TRUE)

Write out results when M = 10000.

unc = data.frame(index = 1:M, tsv, trsv, tnrsv)
write.table(unc, file = "unc.txt")

Load the data from a previous run of the simulation.

L = "unc.txt"
unc = read.table(L, header = TRUE)

Quantifying risk https://understandinguncertainty.org/node/622

Upper quantiles of the MC samples.

exRloss5 = quantile(unc$trsv, prob = .95) * wtRSV # 5%
exRloss1 = quantile(unc$trsv, prob = .99) * wtRSV # 1%
exRloss01 = quantile(unc$trsv, prob = .999) * wtRSV # .1%
exNRloss5 = quantile(unc$tnrsv, prob = .95) * wtNRSV # 5%
exNRloss1 = quantile(unc$tnrsv, prob = .99) * wtNRSV # 1%
exNRloss01 = quantile(unc$tnrsv, prob = .999) * wtNRSV # .1%

exRloss5; exRloss1; exRloss01; exNRloss5; exNRloss1; exNRloss01

##       95% 
## 122686432

##       99% 
## 231167653

##     99.9% 
## 442318422

##      95% 
## 78481301

##       99% 
## 195769481

##     99.9% 
## 606211564

exRloss5 + exNRloss5; exRloss1 + exNRloss1; exRloss01 + exNRloss01

##       95% 
## 201167733

##       99% 
## 426937134

##      99.9% 
## 1048529987

There is a 5% chance that the annual exposure will exceed $203 million, a 1% chance that it will exceed $430 million, and a .1% chance that it will exceed $1 billion.

Plot a histogram.

unc$trsl = unc$trsv * wtRSV
unc$tnrsl = unc$tnrsv * wtNRSV
unc$ttl = unc$trsl + unc$tnrsl
max(unc$ttl)

## [1] 3135297732

ggplot(unc, aes(ttl)) +
  geom_density(fill = "grey") +
  scale_x_log10(label = scales::dollar) +
  geom_vline(xintercept = AvgSL, color = "red") +
  geom_vline(xintercept = 30750000, color = "blue") +
  xlab("Annual Florida Property Loss from Tornadoes\n [2014 Dollars]") +
  theme_minimal() 

as.data.frame(TornP2) %>%
  filter(yr >= 2007) %>%
  summarise(aveL = sum(loss)/9)

##    aveL
## 1 30.75

Chance that the annual loss will exceed X dollars.

probs = seq(.98, 1, .0001)
qts = quantile(unc$ttl, probs)
qts.df = data.frame(qts, pr = 1 - probs)
ggplot(qts.df, aes(y = qts, x = pr)) +
  geom_line(size = 2) +
  scale_x_reverse(label = scales::percent) +
  scale_y_log10(label = scales::dollar, 
                breaks = c(100000000, 200000000, 500000000, 1000000000, 2000000000, 5000000000), 
                limits = c(100000000, 5000000000)) +
  ylab("Property Loss") + xlab("Annual Exceedance Probability") +
  theme_minimal()