Hurricane Data Provisions, R Code Demonstration
Description
This document demonstrates how to use R to determine HIP-WI county level triggers.
Important Information
If you would like to run this code for a specific storm, search for
Dorian, 2019 within the code and replace it with the
desired storm.
If you are a Windows user, you may need to install Rtools in addition to R in order to run this example.
The Code
The following sections break up this demonstration into parts. The complete code without breaks is available at the end of the document.
Loading Packages
The following code downloads the necessary R packages to a temporary
directory, and loads for the current R session. If you run into issues
running this section, you can always do a traditional install
install.packages("package_name_here") and loading
library(package_name_here).
If you use this traditional method, remove
Sys.setenv(TMPDIR=tempdir()) and
if (!require("pacman")) install.packages("pacman") from the
code, and restart your R session.
Sys.setenv(TMPDIR=tempdir())
if (!require("pacman")) install.packages("pacman")
pacman::p_install_version("tidyverse", "1.3.1")
pacman::p_load(tidyverse)
pacman::p_load(lubridate)
pacman::p_install_version("pins", "1.0.1")
pacman::p_load(pins)
pacman::p_install_version("sf", "1.0-7")
pacman::p_load(sf)
pacman::p_install_version("tigris", "1.6.1")
pacman::p_load(tigris)
pacman::p_install_version("geosphere", "1.5-14")
pacman::p_load(geosphere)
pacman::p_install_version("units", "0.8-0")
pacman::p_load(units)
pacman::p_install_version("janitor", "2.1.0")
pacman::p_load(janitor)
pacman::p_install_version("curl", "4.3.2")
pacman::p_load(curl)
pacman::p_install_version("rnaturalearth", "0.1.0")
pacman::p_load(rnaturalearth)
pacman::p_install_version("kableExtra", "1.3.4")
pacman::p_load(kableExtra)
pacman::p_install_version("lutz", "0.3.1")
pacman::p_load(lutz)
theme_set(theme_bw())Prepare Shape Files
The following code prepares shape files for mapping and determining county triggers.
counties_t <- counties(
state = c(
'01',
'05',
'09',
'10',
'12',
'13',
'15',
'22',
'23',
'24',
'25',
'28',
'33',
'34',
'36',
'37',
'42',
'44',
'45',
'48',
'50',
'51'
),
cb = FALSE,
year = 2020,
progress_bar = FALSE,
class = "sf"
) %>%
mutate(
center_point = st_centroid(geometry),
county_timezone = tz_lookup(center_point, method = "accurate")
) %>%
janitor::clean_names() %>%
select(
namelsad,
geoid,
geometry,
county_timezone
) %>%
rename(
county_name = namelsad
)
counties_4326 <- st_transform(counties_t, crs=4326) %>%
st_sf(crs=4326)
remove(counties_t)
## Altering the Monroe County Florida Shapefile
## "Any part of Monroe County, Florida south of Latitude 25˚ North, as identified by the United States Census Bureau, will not be used for the purpose of determining an adjacent county loss trigger under the Hurricane Insurance Protection-Wind Index Endorsement."
counties_4326_12087 <- counties_4326 %>%
filter(
geoid == "12087"
)
bbox_original_12087 <-st_bbox(counties_4326_12087)
bbox_crop_12087 <- bbox_original_12087
bbox_crop_12087[[2]] <- 25
counties_4326_12087_crop <- st_crop(
counties_4326_12087 ,
bbox_crop_12087
) %>%
mutate(
geometry = geometry %>% st_sfc(crs=4326)
) %>%
st_sf(crs=4326)
counties_4326_12087_psudo <- st_difference(
counties_4326_12087 %>% st_make_valid(),
counties_4326_12087_crop %>% st_make_valid()
) %>%
select(
county_name,
geoid,
geometry,
county_timezone
) %>%
mutate(
county_name = "Monroe Keys County",
geoid = "12087.2",
geoid_neighbor = "12087"
)
## Altering the Honolulu County Shapefile
## "The trailing islands in the state of Hawaii, defined by the United States Census Bureau as Census Tract 114.98 of Honolulu County, shall not constitute a triggered county or county adjacent to a triggered county for purposes of determining a county loss trigger under the Hurricane Insurance Protection-Wind Index Endorsement."
counties_4326_15003 <- counties_4326 %>%
filter(
geoid == "15003"
)
bbox_original_15003 <-st_bbox(counties_4326_15003)
bbox_crop_15003 <- bbox_original_15003
bbox_crop_15003[[4]] <- 22
counties_4326_15003_crop <- st_crop(
counties_4326_15003,
bbox_crop_15003
)
## Updating the County File
counties_4326_update <- counties_4326 %>%
mutate(
geometry = ifelse(
geoid == "12087",
counties_4326_12087_crop$geometry,
ifelse(
geoid == "15003",
counties_4326_15003_crop$geometry,
geometry
)
) %>%
st_sfc(crs=4326)
) %>%
st_sf(crs=4326)
## US Census Bureau County Adjacency File
county_adjacency_pin <- pin("https://www2.census.gov/geo/docs/reference/county_adjacency.txt")
county_adjacency <- read_delim(
county_adjacency_pin,
"\t",
escape_double = FALSE,
trim_ws = TRUE,
col_names = FALSE
) %>%
select(
X2,
X4
) %>%
rename(
geoid = X2,
geoid_neighbor = X4
) %>%
fill(geoid, .direction="down") %>%
group_by(geoid) %>%
nest() %>%
ungroup()
counties_4326_update <- left_join(
counties_4326_update %>%
select(
county_name,
geoid,
county_timezone
),
county_adjacency
) %>%
full_join(
counties_4326_12087_psudo %>%
tibble() %>%
nest(data=c("geoid_neighbor"))
) %>%
drop_na(
geometry,
data
)
counties_5070_update <- counties_4326_update %>% st_transform(crs=5070)Download Hurricane Data
The following temporarily downloads International Best Track Archive for Climate Stewardship (IBTrACS) data from NOAA. If you are going to run this program more than once, consider downloading the data locally to speed up processing and reduce burden on NOAA servers.
# IBTrACS
ibtracs_shape_zip <- tempfile()
curl_download(
'https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/shapefile/IBTrACS.since1980.list.v04r00.points.zip',
ibtracs_shape_zip
)
ibtracs_shape_unzip <- tempfile()
unzip(
zipfile = ibtracs_shape_zip,
exdir = ibtracs_shape_unzip
)
remove(ibtracs_shape_zip)
ibtracs_shape <- list.files(
ibtracs_shape_unzip,
pattern = ".shp$",
full.names = TRUE
)
remove(ibtracs_shape_unzip)
ibtracs_4326 <- read_sf(ibtracs_shape) %>%
clean_names() %>%
select(
usa_atcfid,
usa_wind,
usa_pres,
name,
season,
basin,
iso_time,
starts_with("usa_r64"),
track_type
) %>%
filter(
basin %in% c("NA", "EP")
) %>%
mutate(
name = name %>% str_to_title(),
storm_name = str_c(
name,
usa_atcfid %>% str_trunc(4, side = "left", ellipsis = ""),
sep = ", "
) %>%
str_trim(),
season = season %>% as.numeric(),
iso_time = iso_time %>% as_datetime()
) %>%
select(-name) %>%
mutate(
across(
starts_with("usa_r"),
as.numeric
)
) %>%
filter(
storm_name=="Dorian, 2019"
)
remove(ibtracs_shape)Construct Hurricane Corridors
The following code creates hurricane corridors using the methods described in the HDP.
# Constructing Hurricane Corridors
bearing_map <- function(geometry, lon_lead, lat_lead) {
bearing(
st_coordinates(geometry),
c(lon_lead, lat_lead)
) %>%
set_units(degrees)
}
## Another storm can be selected here
points_4326 <- ibtracs_4326 %>%
arrange(iso_time) %>%
group_by(usa_atcfid) %>%
mutate(
lat = map(geometry, 2) %>% unlist(),
lon = map(geometry, 1) %>% unlist(),
lat_lead = lead(lat),
lon_lead = lead(lon),
lat_lag = lag(lat),
lon_lag = lag(lon),
usa_wind_lead = lead(usa_wind),
usa_wind_lag = lag(usa_wind),
point_prior_hu = if_else(
usa_wind_lead >= 64 & usa_wind < 64,
1,
0
),
point_after_hu = if_else(
usa_wind_lag >= 64 & usa_wind < 64,
1,
0
),
hu = if_else(usa_wind >= 64, 1, 0) %>% as_factor,
# Nautical Miles to Meters
wind_extent_hu_m = pmax(
usa_r64_ne,
usa_r64_se,
usa_r64_nw,
usa_r64_sw,
na.rm = T
)*1852,
bearing_accurate = pmap(
.l = list(
geometry=geometry,
lon_lead=lon_lead,
lat_lead=lat_lead
),
.f = bearing_map
)
) %>%
ungroup() %>%
rowwise() %>%
mutate(
distance_to_next = distGeo(
c(lon, lat),
c(lon_lead, lat_lead)
),
wind_extent_hu_m = wind_extent_hu_m %>% replace_na(0)
) %>%
ungroup()
buffers_5070 <- points_4326 %>%
st_transform(crs=5070) %>%
ungroup() %>%
mutate(
geometry = st_buffer(
st_geometry(.),
dist = wind_extent_hu_m,
nQuadSegs = 180
)
) %>%
ungroup()
## Interpolating Hurricane Status Transition Points
### Hurricane Dissolution (end points)
points_4326_e1 <- points_4326 %>%
arrange(iso_time) %>%
group_by(usa_atcfid) %>%
mutate(
lag_wind_extent_hu_m = lag(wind_extent_hu_m),
lag_distance_to_next = lag(distance_to_next),
bearing_accurate = lag(bearing_accurate),
lag_iso_time = lag(iso_time)
) %>%
ungroup() %>%
filter(point_after_hu==1) %>%
mutate(
mid_distance_to_lead = lag_distance_to_next*((usa_wind_lag - 64)/(usa_wind_lag - usa_wind)),
wind_extent_hu_m_guess_last =
pmax(
lag_wind_extent_hu_m/2,
lag_wind_extent_hu_m*(1 -((usa_wind_lag - 64)/(usa_wind_lag - usa_wind)))
),
iso_time = iso_time - minutes(
round(
difftime(
iso_time,
lag_iso_time,
units = "mins"
)*(mid_distance_to_lead/lag_distance_to_next)
)
),
wind_extent_hu_m = wind_extent_hu_m_guess_last,
interpolated_hu = T,
hu = 1 %>% as_factor,
distance_to_next = mid_distance_to_lead,
usa_wind = 64
) %>%
arrange(iso_time) %>%
select(
usa_atcfid,
storm_name,
season,
iso_time,
wind_extent_hu_m,
bearing_accurate,
lon,
lon_lag,
lat,
lat_lag,
usa_wind,
interpolated_hu,
hu,
distance_to_next,
basin
) %>%
ungroup()
points_4326_e2 <- points_4326_e1
for (i in 1:dim(points_4326_e2)[1]) {
points_4326_e2$geometry[i] <- destPoint(
c(
points_4326_e2$lon_lag[i],
points_4326_e2$lat_lag[i]
),
points_4326_e2$bearing_accurate[i],
points_4326_e2$distance_to_next[i]
) %>% st_point() %>% st_sfc()
}
interpolated_end_hu_5070 <- points_4326_e2 %>%
select(
-lat,
-lat_lag,
-lon,
-lon_lag
) %>%
st_transform(crs=5070) %>%
ungroup() %>%
mutate(
geometry = st_buffer(
st_geometry(.),
dist = wind_extent_hu_m,
nQuadSegs = 180
)
) %>%
ungroup()
### Hurricane Inception (start points)
points_4326_s1 <- points_4326 %>%
arrange(iso_time) %>%
group_by(usa_atcfid) %>%
mutate(
lead_wind_extent_hu_m = lead(wind_extent_hu_m),
lead_iso_time = lead(iso_time)
) %>%
filter(point_prior_hu==1) %>%
ungroup() %>%
mutate(
mid_distance_to_lead = distance_to_next*(1-(usa_wind_lead - 64)/(usa_wind_lead - usa_wind)),
wind_extent_hu_m_guess_first =
pmax(
lead_wind_extent_hu_m/2,
lead_wind_extent_hu_m*(1 -((usa_wind_lead - 64)/(usa_wind_lead - usa_wind)))
),
iso_time = iso_time + minutes(
round(
difftime(
lead_iso_time,
iso_time,
units = "mins"
) *(mid_distance_to_lead/distance_to_next)
)
),
wind_extent_hu_m = wind_extent_hu_m_guess_first,
interpolated_hu = T,
hu = 1 %>% as_factor,
distance_to_next = mid_distance_to_lead,
usa_wind = 64,
) %>%
arrange(iso_time) %>%
select(
usa_atcfid,
storm_name,
season,
iso_time,
wind_extent_hu_m,
bearing_accurate,
lon,
lat,
distance_to_next,
usa_wind,
interpolated_hu,
hu,
basin
)
points_4326_s2 <- points_4326_s1
for (i in 1:dim(points_4326_s2)[1]) {
points_4326_s2$geometry[i] <- destPoint(
c(
points_4326_s2$lon[i],
points_4326_s2$lat[i]
),
points_4326_s2$bearing_accurate[i],
points_4326_s2$distance_to_next[i]
) %>% st_point() %>% st_sfc()
}
interpolated_start_hu_5070 <- points_4326_s2 %>%
select(
-lat,
-lon
) %>%
st_transform(crs=5070) %>%
ungroup() %>%
mutate(
geometry = st_buffer(
st_geometry(.),
dist = wind_extent_hu_m,
nQuadSegs = 180
)
) %>%
ungroup()
remove(points_4326)
### Combine Interpolated Data
buffers_5070_new <- buffers_5070 %>%
as_tibble() %>%
mutate(
interpolated_hu = F
) %>%
full_join(
interpolated_start_hu_5070 %>%
as_tibble()
) %>%
full_join(
interpolated_end_hu_5070 %>%
as_tibble()
) %>%
st_sf(crs=5070)
remove(buffers_5070)
### Creating Polygon Pairs
status_change_function <- function(x) {
x <- rle(x)$lengths
rep(seq_along(x), times=x)
}
poly_5070 <- buffers_5070_new %>%
ungroup() %>%
arrange(usa_atcfid, iso_time) %>%
mutate(
hu_status_change = status_change_function(hu %>% as.numeric())
) %>%
group_by(usa_atcfid, hu, hu_status_change) %>%
mutate(
geometry_lag = lag(geometry) %>% st_sfc()
) %>%
slice(2:n()) %>%
ungroup()
geo_union_function <- function(x, y, ...) st_union(x, y)
poly_5070$geometry <- map2(
poly_5070$geometry,
poly_5070$geometry_lag,
geo_union_function
)
poly_5070 <- poly_5070 %>%
select(-geometry_lag) %>%
mutate(
geometry = st_sfc(geometry, crs = 5070)
)
### Convex Hull
hull_5070 <- poly_5070 %>%
group_by(usa_atcfid, iso_time) %>%
mutate(
geometry = st_convex_hull(geometry)
) %>%
ungroup()
remove(poly_5070)Find Triggered Counties
The following code finds the intersection of the hurricane corridor and the county shape files prepared in previous steps.
# Finding Triggered Counties
triggered_direct <- st_intersection(
hull_5070,
counties_5070_update %>%
st_make_valid()
) %>%
tibble() %>%
select(-geometry) %>%
arrange(usa_atcfid, iso_time) %>%
group_by(usa_atcfid, geoid) %>%
ungroup() %>%
select(
usa_atcfid,
geoid,
county_name,
data,
storm_name,
season,
iso_time,
interpolated_hu
)
triggered_direct_4326_t <- triggered_direct %>%
select(
usa_atcfid,
geoid,
iso_time,
county_name
) %>%
group_by(usa_atcfid, geoid) %>%
slice_min(iso_time) %>%
ungroup() %>%
left_join(
counties_4326_update
) %>%
st_sf(crs=4326)
triggered <- triggered_direct %>%
unnest("data")%>%
mutate(
trigger = if_else(
geoid == geoid_neighbor,
"direct",
"indirect"
),
county_name = ifelse(
trigger == "indirect",
NA,
county_name
),
geoid = if_else(
geoid == geoid_neighbor,
geoid,
geoid_neighbor
)
) %>%
select(
-geoid_neighbor,
-county_name
) %>%
left_join(
counties_4326_update %>%
tibble() %>%
select(
county_name,
geoid,
county_timezone
)
) %>%
mutate(
local_time = iso_time %>% with_tz(county_timezone)
) %>%
select(
-iso_time
) %>%
group_by(
usa_atcfid,
geoid,
storm_name,
county_name,
county_timezone
) %>%
summarize(
start_local_time = min(local_time),
end_local_time = max(local_time),
) %>%
ungroup() %>%
mutate(
trigger_interval_local_time = interval(start_local_time, end_local_time)
)
triggered_4326 <- triggered %>%
left_join(
counties_4326 %>%
tibble() %>%
select(
geoid,
geometry
)
) %>%
filter(
geoid != 12087.2
) %>%
st_sf(crs=4326)Triggered Counties Map
This code creates as basic map of the hurricane corridor and the triggered counties.
# Map
north_america_4326 <- ne_countries(
continent = 'north america',
returnclass = "sf"
) %>%
st_transform(crs=4326)
ggplot() +
geom_sf(
data = north_america_4326,
fill = "#f7f7f7"
) +
geom_sf(
data = counties_4326,
fill = "#f7f7f7"
) +
geom_sf(
data = triggered_4326,
fill = "#67a9cf"
) +
geom_sf(
data = hull_5070 %>%
st_transform(crs=4326) %>%
st_union() %>%
st_crop(
xmin = triggered_4326 %>%
st_bbox() %>%
pluck(1) - .3,
ymin = triggered_4326 %>%
st_bbox() %>%
pluck(2) - .3,
xmax = triggered_4326 %>%
st_bbox() %>%
pluck(3) + .3,
ymax = triggered_4326 %>%
st_bbox() %>%
pluck(4) + .3,
),
color = "#ef8a62",
alpha = .5,
fill = "#ef8a62",
size = 1,
show.legend = FALSE
) +
geom_sf(
data = triggered_direct_4326_t,
fill = NA
) +
geom_sf(
data = hull_5070 %>%
st_transform(crs=4326) %>%
st_union() %>%
st_crop(
xmin = triggered_4326 %>%
st_bbox() %>%
pluck(1) - .3,
ymin = triggered_4326 %>%
st_bbox() %>%
pluck(2) - .3,
xmax = triggered_4326 %>%
st_bbox() %>%
pluck(3) + .3,
ymax = triggered_4326 %>%
st_bbox() %>%
pluck(4) + .3,
),
color = "#ef8a62",
fill = NA,
size = 1,
show.legend = FALSE
) +
coord_sf(
xlim =
c(
triggered_4326 %>%
st_bbox() %>%
pluck(1),
triggered_4326 %>%
st_bbox() %>%
pluck(3)
),
ylim =
c(
triggered_4326 %>%
st_bbox() %>%
pluck(2),
triggered_4326 %>%
st_bbox() %>%
pluck(4)
)
) +
theme(
panel.background = element_rect(fill = "aliceblue")
) +
ggtitle(
paste0(
"Hurricane ",
hull_5070$storm_name %>% first()
),
subtitle = paste0("HIP-WI County Triggers")
) 
Triggered Counties Table
This section produces a table that lists all of the triggered counties, along with the county’s trigger period. For a HIP-WI policy to trigger, the policy must be active during some part of the trigger interval.
## Triggered Counties
triggered_4326 %>%
tibble() %>%
mutate(
county_name = county_name %>% str_to_title()
) %>%
select(
county_name,
geoid,
storm_name,
trigger_interval_local_time
) %>%
kbl(
col.names = c("County", "FIP", "Storm Name", "Trigger Interval Local Time")
) %>%
kable_styling(
bootstrap_options = c(
"condensed",
"responsive"
)
)| County | FIP | Storm Name | Trigger Interval Local Time |
|---|---|---|---|
| Beaufort County | 37013 | Dorian, 2019 | 2019-09-06 02:00:00 EDT–2019-09-06 14:00:00 EDT |
| Bertie County | 37015 | Dorian, 2019 | 2019-09-06 08:00:00 EDT–2019-09-06 11:00:00 EDT |
| Bladen County | 37017 | Dorian, 2019 | 2019-09-05 20:00:00 EDT–2019-09-06 08:00:00 EDT |
| Brunswick County | 37019 | Dorian, 2019 | 2019-09-05 17:00:00 EDT–2019-09-06 08:00:00 EDT |
| Camden County | 37029 | Dorian, 2019 | 2019-09-06 08:00:00 EDT–2019-09-06 14:00:00 EDT |
| Carteret County | 37031 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 14:00:00 EDT |
| Chowan County | 37041 | Dorian, 2019 | 2019-09-06 08:00:00 EDT–2019-09-06 11:00:00 EDT |
| Columbus County | 37047 | Dorian, 2019 | 2019-09-05 17:00:00 EDT–2019-09-06 08:00:00 EDT |
| Craven County | 37049 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 11:00:00 EDT |
| Currituck County | 37053 | Dorian, 2019 | 2019-09-06 05:00:00 EDT–2019-09-06 14:00:00 EDT |
| Dare County | 37055 | Dorian, 2019 | 2019-09-06 05:00:00 EDT–2019-09-06 14:00:00 EDT |
| Duplin County | 37061 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 08:00:00 EDT |
| Gates County | 37073 | Dorian, 2019 | 2019-09-06 08:30:00 EDT–2019-09-06 11:00:00 EDT |
| Hyde County | 37095 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 14:00:00 EDT |
| Jones County | 37103 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 11:00:00 EDT |
| Lenoir County | 37107 | Dorian, 2019 | 2019-09-06 02:00:00 EDT–2019-09-06 11:00:00 EDT |
| Martin County | 37117 | Dorian, 2019 | 2019-09-06 05:00:00 EDT–2019-09-06 11:00:00 EDT |
| New Hanover County | 37129 | Dorian, 2019 | 2019-09-05 17:00:00 EDT–2019-09-06 08:00:00 EDT |
| Onslow County | 37133 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 11:00:00 EDT |
| Pamlico County | 37137 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 14:00:00 EDT |
| Pasquotank County | 37139 | Dorian, 2019 | 2019-09-06 08:00:00 EDT–2019-09-06 14:00:00 EDT |
| Pender County | 37141 | Dorian, 2019 | 2019-09-05 17:00:00 EDT–2019-09-06 08:00:00 EDT |
| Perquimans County | 37143 | Dorian, 2019 | 2019-09-06 08:00:00 EDT–2019-09-06 14:00:00 EDT |
| Pitt County | 37147 | Dorian, 2019 | 2019-09-06 02:00:00 EDT–2019-09-06 11:00:00 EDT |
| Robeson County | 37155 | Dorian, 2019 | 2019-09-05 17:00:00 EDT–2019-09-06 02:00:00 EDT |
| Sampson County | 37163 | Dorian, 2019 | 2019-09-05 23:00:00 EDT–2019-09-06 08:00:00 EDT |
| Tyrrell County | 37177 | Dorian, 2019 | 2019-09-06 05:00:00 EDT–2019-09-06 14:00:00 EDT |
| Washington County | 37187 | Dorian, 2019 | 2019-09-06 05:00:00 EDT–2019-09-06 14:00:00 EDT |
| Wayne County | 37191 | Dorian, 2019 | 2019-09-06 02:00:00 EDT–2019-09-06 05:00:00 EDT |
| Berkeley County | 45015 | Dorian, 2019 | 2019-09-05 08:00:00 EDT–2019-09-05 20:00:00 EDT |
| Charleston County | 45019 | Dorian, 2019 | 2019-09-05 08:00:00 EDT–2019-09-05 20:00:00 EDT |
| Clarendon County | 45027 | Dorian, 2019 | 2019-09-05 11:00:00 EDT–2019-09-05 17:00:00 EDT |
| Colleton County | 45029 | Dorian, 2019 | 2019-09-05 08:00:00 EDT–2019-09-05 20:00:00 EDT |
| Dillon County | 45033 | Dorian, 2019 | 2019-09-05 17:00:00 EDT–2019-09-06 02:00:00 EDT |
| Dorchester County | 45035 | Dorian, 2019 | 2019-09-05 08:00:00 EDT–2019-09-05 20:00:00 EDT |
| Georgetown County | 45043 | Dorian, 2019 | 2019-09-05 08:00:00 EDT–2019-09-06 02:00:00 EDT |
| Horry County | 45051 | Dorian, 2019 | 2019-09-05 11:00:00 EDT–2019-09-06 05:00:00 EDT |
| Marion County | 45067 | Dorian, 2019 | 2019-09-05 11:00:00 EDT–2019-09-06 02:00:00 EDT |
| Orangeburg County | 45075 | Dorian, 2019 | 2019-09-05 11:00:00 EDT–2019-09-05 17:00:00 EDT |
| Williamsburg County | 45089 | Dorian, 2019 | 2019-09-05 11:00:00 EDT–2019-09-05 20:00:00 EDT |
| Chesapeake City | 51550 | Dorian, 2019 | 2019-09-06 08:00:00 EDT–2019-09-06 14:00:00 EDT |
| Suffolk City | 51800 | Dorian, 2019 | 2019-09-06 08:30:00 EDT–2019-09-06 11:00:00 EDT |
| Virginia Beach City | 51810 | Dorian, 2019 | 2019-09-06 08:00:00 EDT–2019-09-06 14:00:00 EDT |
Session Information
si <- sessionInfo()
print(si,loadedOnly=T, locale = FALSE)## R version 4.1.3 (2022-03-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lutz_0.3.1 kableExtra_1.3.4 rnaturalearth_0.1.0
## [4] curl_4.3.2 janitor_2.1.0 units_0.8-0
## [7] geosphere_1.5-14 tigris_1.6.1 sf_1.0-7
## [10] pins_1.0.1 lubridate_1.8.0 forcats_0.5.1
## [13] stringr_1.4.0 dplyr_1.0.9 purrr_0.3.4
## [16] readr_2.1.2 tidyr_1.2.0 tibble_3.1.7
## [19] ggplot2_3.3.6 tidyverse_1.3.1 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.2 bit64_4.0.5 filelock_1.0.2 webshot_0.5.3
## [5] httr_1.4.3 tools_4.1.3 backports_1.4.1 bslib_0.3.1
## [9] utf8_1.2.2 rgdal_1.5-32 R6_2.5.1 KernSmooth_2.23-20
## [13] DBI_1.1.3 colorspace_2.0-3 withr_2.5.0 sp_1.5-0
## [17] tidyselect_1.1.2 bit_4.0.4 compiler_4.1.3 cli_3.3.0
## [21] rvest_1.0.2 xml2_1.3.3 sass_0.4.1 scales_1.2.0
## [25] classInt_0.4-7 proxy_0.4-27 rappdirs_0.3.3 systemfonts_1.0.4
## [29] digest_0.6.29 foreign_0.8-82 rmarkdown_2.14 svglite_2.1.0
## [33] pkgconfig_2.0.3 htmltools_0.5.2 highr_0.9 dbplyr_2.2.0
## [37] fastmap_1.1.0 rlang_1.0.2 readxl_1.4.0 rstudioapi_0.13
## [41] farver_2.1.0 jquerylib_0.1.4 generics_0.1.2 jsonlite_1.8.0
## [45] vroom_1.5.7 magrittr_2.0.3 s2_1.0.7 Rcpp_1.0.8.3
## [49] munsell_0.5.0 fansi_1.0.3 lifecycle_1.0.1 stringi_1.7.6
## [53] yaml_2.3.5 snakecase_0.11.0 grid_4.1.3 maptools_1.1-4
## [57] parallel_4.1.3 crayon_1.5.1 lattice_0.20-45 haven_2.5.0
## [61] hms_1.1.1 knitr_1.39 pillar_1.7.0 uuid_1.1-0
## [65] wk_0.6.0 reprex_2.0.1 glue_1.6.2 evaluate_0.15
## [69] modelr_0.1.8 vctrs_0.4.1 tzdb_0.3.0 cellranger_1.1.0
## [73] gtable_0.3.0 assertthat_0.2.1 xfun_0.31 broom_0.8.0
## [77] e1071_1.7-11 class_7.3-20 viridisLite_0.4.0 ellipsis_0.3.2All of the Code
Sys.setenv(TMPDIR=tempdir())
if (!require("pacman")) install.packages("pacman")
pacman::p_install_version("tidyverse", "1.3.1")
pacman::p_load(tidyverse)
pacman::p_load(lubridate)
pacman::p_install_version("pins", "1.0.1")
pacman::p_load(pins)
pacman::p_install_version("sf", "1.0-7")
pacman::p_load(sf)
pacman::p_install_version("tigris", "1.6.1")
pacman::p_load(tigris)
pacman::p_install_version("geosphere", "1.5-14")
pacman::p_load(geosphere)
pacman::p_install_version("units", "0.8-0")
pacman::p_load(units)
pacman::p_install_version("janitor", "2.1.0")
pacman::p_load(janitor)
pacman::p_install_version("curl", "4.3.2")
pacman::p_load(curl)
pacman::p_install_version("rnaturalearth", "0.1.0")
pacman::p_load(rnaturalearth)
pacman::p_install_version("kableExtra", "1.3.4")
pacman::p_load(kableExtra)
pacman::p_install_version("lutz", "0.3.1")
pacman::p_load(lutz)
theme_set(theme_bw())
counties_t <- counties(
state = c(
'01',
'05',
'09',
'10',
'12',
'13',
'15',
'22',
'23',
'24',
'25',
'28',
'33',
'34',
'36',
'37',
'42',
'44',
'45',
'48',
'50',
'51'
),
cb = FALSE,
year = 2020,
progress_bar = FALSE,
class = "sf"
) %>%
mutate(
center_point = st_centroid(geometry),
county_timezone = tz_lookup(center_point, method = "accurate")
) %>%
janitor::clean_names() %>%
select(
namelsad,
geoid,
geometry,
county_timezone
) %>%
rename(
county_name = namelsad
)
counties_4326 <- st_transform(counties_t, crs=4326) %>%
st_sf(crs=4326)
remove(counties_t)
## Altering the Monroe County Florida Shapefile
## "Any part of Monroe County, Florida south of Latitude 25˚ North, as identified by the United States Census Bureau, will not be used for the purpose of determining an adjacent county loss trigger under the Hurricane Insurance Protection-Wind Index Endorsement."
counties_4326_12087 <- counties_4326 %>%
filter(
geoid == "12087"
)
bbox_original_12087 <-st_bbox(counties_4326_12087)
bbox_crop_12087 <- bbox_original_12087
bbox_crop_12087[[2]] <- 25
counties_4326_12087_crop <- st_crop(
counties_4326_12087 ,
bbox_crop_12087
) %>%
mutate(
geometry = geometry %>% st_sfc(crs=4326)
) %>%
st_sf(crs=4326)
counties_4326_12087_psudo <- st_difference(
counties_4326_12087 %>% st_make_valid(),
counties_4326_12087_crop %>% st_make_valid()
) %>%
select(
county_name,
geoid,
geometry,
county_timezone
) %>%
mutate(
county_name = "Monroe Keys County",
geoid = "12087.2",
geoid_neighbor = "12087"
)
## Altering the Honolulu County Shapefile
## "The trailing islands in the state of Hawaii, defined by the United States Census Bureau as Census Tract 114.98 of Honolulu County, shall not constitute a triggered county or county adjacent to a triggered county for purposes of determining a county loss trigger under the Hurricane Insurance Protection-Wind Index Endorsement."
counties_4326_15003 <- counties_4326 %>%
filter(
geoid == "15003"
)
bbox_original_15003 <-st_bbox(counties_4326_15003)
bbox_crop_15003 <- bbox_original_15003
bbox_crop_15003[[4]] <- 22
counties_4326_15003_crop <- st_crop(
counties_4326_15003,
bbox_crop_15003
)
## Updating the County File
counties_4326_update <- counties_4326 %>%
mutate(
geometry = ifelse(
geoid == "12087",
counties_4326_12087_crop$geometry,
ifelse(
geoid == "15003",
counties_4326_15003_crop$geometry,
geometry
)
) %>%
st_sfc(crs=4326)
) %>%
st_sf(crs=4326)
## US Census Bureau County Adjacency File
county_adjacency_pin <- pin("https://www2.census.gov/geo/docs/reference/county_adjacency.txt")
county_adjacency <- read_delim(
county_adjacency_pin,
"\t",
escape_double = FALSE,
trim_ws = TRUE,
col_names = FALSE
) %>%
select(
X2,
X4
) %>%
rename(
geoid = X2,
geoid_neighbor = X4
) %>%
fill(geoid, .direction="down") %>%
group_by(geoid) %>%
nest() %>%
ungroup()
counties_4326_update <- left_join(
counties_4326_update %>%
select(
county_name,
geoid,
county_timezone
),
county_adjacency
) %>%
full_join(
counties_4326_12087_psudo %>%
tibble() %>%
nest(data=c("geoid_neighbor"))
) %>%
drop_na(
geometry,
data
)
counties_5070_update <- counties_4326_update %>% st_transform(crs=5070)
# IBTrACS
ibtracs_shape_zip <- tempfile()
curl_download(
'https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/shapefile/IBTrACS.since1980.list.v04r00.points.zip',
ibtracs_shape_zip
)
ibtracs_shape_unzip <- tempfile()
unzip(
zipfile = ibtracs_shape_zip,
exdir = ibtracs_shape_unzip
)
remove(ibtracs_shape_zip)
ibtracs_shape <- list.files(
ibtracs_shape_unzip,
pattern = ".shp$",
full.names = TRUE
)
remove(ibtracs_shape_unzip)
ibtracs_4326 <- read_sf(ibtracs_shape) %>%
clean_names() %>%
select(
usa_atcfid,
usa_wind,
usa_pres,
name,
season,
basin,
iso_time,
starts_with("usa_r64"),
track_type
) %>%
filter(
basin %in% c("NA", "EP")
) %>%
mutate(
name = name %>% str_to_title(),
storm_name = str_c(
name,
usa_atcfid %>% str_trunc(4, side = "left", ellipsis = ""),
sep = ", "
) %>%
str_trim(),
season = season %>% as.numeric(),
iso_time = iso_time %>% as_datetime()
) %>%
select(-name) %>%
mutate(
across(
starts_with("usa_r"),
as.numeric
)
) %>%
filter(
storm_name=="Dorian, 2019"
)
remove(ibtracs_shape)
# Constructing Hurricane Corridors
bearing_map <- function(geometry, lon_lead, lat_lead) {
bearing(
st_coordinates(geometry),
c(lon_lead, lat_lead)
) %>%
set_units(degrees)
}
## Another storm can be selected here
points_4326 <- ibtracs_4326 %>%
arrange(iso_time) %>%
group_by(usa_atcfid) %>%
mutate(
lat = map(geometry, 2) %>% unlist(),
lon = map(geometry, 1) %>% unlist(),
lat_lead = lead(lat),
lon_lead = lead(lon),
lat_lag = lag(lat),
lon_lag = lag(lon),
usa_wind_lead = lead(usa_wind),
usa_wind_lag = lag(usa_wind),
point_prior_hu = if_else(
usa_wind_lead >= 64 & usa_wind < 64,
1,
0
),
point_after_hu = if_else(
usa_wind_lag >= 64 & usa_wind < 64,
1,
0
),
hu = if_else(usa_wind >= 64, 1, 0) %>% as_factor,
# Nautical Miles to Meters
wind_extent_hu_m = pmax(
usa_r64_ne,
usa_r64_se,
usa_r64_nw,
usa_r64_sw,
na.rm = T
)*1852,
bearing_accurate = pmap(
.l = list(
geometry=geometry,
lon_lead=lon_lead,
lat_lead=lat_lead
),
.f = bearing_map
)
) %>%
ungroup() %>%
rowwise() %>%
mutate(
distance_to_next = distGeo(
c(lon, lat),
c(lon_lead, lat_lead)
),
wind_extent_hu_m = wind_extent_hu_m %>% replace_na(0)
) %>%
ungroup()
buffers_5070 <- points_4326 %>%
st_transform(crs=5070) %>%
ungroup() %>%
mutate(
geometry = st_buffer(
st_geometry(.),
dist = wind_extent_hu_m,
nQuadSegs = 180
)
) %>%
ungroup()
## Interpolating Hurricane Status Transition Points
### Hurricane Dissolution (end points)
points_4326_e1 <- points_4326 %>%
arrange(iso_time) %>%
group_by(usa_atcfid) %>%
mutate(
lag_wind_extent_hu_m = lag(wind_extent_hu_m),
lag_distance_to_next = lag(distance_to_next),
bearing_accurate = lag(bearing_accurate),
lag_iso_time = lag(iso_time)
) %>%
ungroup() %>%
filter(point_after_hu==1) %>%
mutate(
mid_distance_to_lead = lag_distance_to_next*((usa_wind_lag - 64)/(usa_wind_lag - usa_wind)),
wind_extent_hu_m_guess_last =
pmax(
lag_wind_extent_hu_m/2,
lag_wind_extent_hu_m*(1 -((usa_wind_lag - 64)/(usa_wind_lag - usa_wind)))
),
iso_time = iso_time - minutes(
round(
difftime(
iso_time,
lag_iso_time,
units = "mins"
)*(mid_distance_to_lead/lag_distance_to_next)
)
),
wind_extent_hu_m = wind_extent_hu_m_guess_last,
interpolated_hu = T,
hu = 1 %>% as_factor,
distance_to_next = mid_distance_to_lead,
usa_wind = 64
) %>%
arrange(iso_time) %>%
select(
usa_atcfid,
storm_name,
season,
iso_time,
wind_extent_hu_m,
bearing_accurate,
lon,
lon_lag,
lat,
lat_lag,
usa_wind,
interpolated_hu,
hu,
distance_to_next,
basin
) %>%
ungroup()
points_4326_e2 <- points_4326_e1
for (i in 1:dim(points_4326_e2)[1]) {
points_4326_e2$geometry[i] <- destPoint(
c(
points_4326_e2$lon_lag[i],
points_4326_e2$lat_lag[i]
),
points_4326_e2$bearing_accurate[i],
points_4326_e2$distance_to_next[i]
) %>% st_point() %>% st_sfc()
}
interpolated_end_hu_5070 <- points_4326_e2 %>%
select(
-lat,
-lat_lag,
-lon,
-lon_lag
) %>%
st_transform(crs=5070) %>%
ungroup() %>%
mutate(
geometry = st_buffer(
st_geometry(.),
dist = wind_extent_hu_m,
nQuadSegs = 180
)
) %>%
ungroup()
### Hurricane Inception (start points)
points_4326_s1 <- points_4326 %>%
arrange(iso_time) %>%
group_by(usa_atcfid) %>%
mutate(
lead_wind_extent_hu_m = lead(wind_extent_hu_m),
lead_iso_time = lead(iso_time)
) %>%
filter(point_prior_hu==1) %>%
ungroup() %>%
mutate(
mid_distance_to_lead = distance_to_next*(1-(usa_wind_lead - 64)/(usa_wind_lead - usa_wind)),
wind_extent_hu_m_guess_first =
pmax(
lead_wind_extent_hu_m/2,
lead_wind_extent_hu_m*(1 -((usa_wind_lead - 64)/(usa_wind_lead - usa_wind)))
),
iso_time = iso_time + minutes(
round(
difftime(
lead_iso_time,
iso_time,
units = "mins"
) *(mid_distance_to_lead/distance_to_next)
)
),
wind_extent_hu_m = wind_extent_hu_m_guess_first,
interpolated_hu = T,
hu = 1 %>% as_factor,
distance_to_next = mid_distance_to_lead,
usa_wind = 64,
) %>%
arrange(iso_time) %>%
select(
usa_atcfid,
storm_name,
season,
iso_time,
wind_extent_hu_m,
bearing_accurate,
lon,
lat,
distance_to_next,
usa_wind,
interpolated_hu,
hu,
basin
)
points_4326_s2 <- points_4326_s1
for (i in 1:dim(points_4326_s2)[1]) {
points_4326_s2$geometry[i] <- destPoint(
c(
points_4326_s2$lon[i],
points_4326_s2$lat[i]
),
points_4326_s2$bearing_accurate[i],
points_4326_s2$distance_to_next[i]
) %>% st_point() %>% st_sfc()
}
interpolated_start_hu_5070 <- points_4326_s2 %>%
select(
-lat,
-lon
) %>%
st_transform(crs=5070) %>%
ungroup() %>%
mutate(
geometry = st_buffer(
st_geometry(.),
dist = wind_extent_hu_m,
nQuadSegs = 180
)
) %>%
ungroup()
remove(points_4326)
### Combine Interpolated Data
buffers_5070_new <- buffers_5070 %>%
as_tibble() %>%
mutate(
interpolated_hu = F
) %>%
full_join(
interpolated_start_hu_5070 %>%
as_tibble()
) %>%
full_join(
interpolated_end_hu_5070 %>%
as_tibble()
) %>%
st_sf(crs=5070)
remove(buffers_5070)
### Creating Polygon Pairs
status_change_function <- function(x) {
x <- rle(x)$lengths
rep(seq_along(x), times=x)
}
poly_5070 <- buffers_5070_new %>%
ungroup() %>%
arrange(usa_atcfid, iso_time) %>%
mutate(
hu_status_change = status_change_function(hu %>% as.numeric())
) %>%
group_by(usa_atcfid, hu, hu_status_change) %>%
mutate(
geometry_lag = lag(geometry) %>% st_sfc()
) %>%
slice(2:n()) %>%
ungroup()
geo_union_function <- function(x, y, ...) st_union(x, y)
poly_5070$geometry <- map2(
poly_5070$geometry,
poly_5070$geometry_lag,
geo_union_function
)
poly_5070 <- poly_5070 %>%
select(-geometry_lag) %>%
mutate(
geometry = st_sfc(geometry, crs = 5070)
)
### Convex Hull
hull_5070 <- poly_5070 %>%
group_by(usa_atcfid, iso_time) %>%
mutate(
geometry = st_convex_hull(geometry)
) %>%
ungroup()
remove(poly_5070)
# Finding Triggered Counties
triggered_direct <- st_intersection(
hull_5070,
counties_5070_update %>%
st_make_valid()
) %>%
tibble() %>%
select(-geometry) %>%
arrange(usa_atcfid, iso_time) %>%
group_by(usa_atcfid, geoid) %>%
ungroup() %>%
select(
usa_atcfid,
geoid,
county_name,
data,
storm_name,
season,
iso_time,
interpolated_hu
)
triggered_direct_4326_t <- triggered_direct %>%
select(
usa_atcfid,
geoid,
iso_time,
county_name
) %>%
group_by(usa_atcfid, geoid) %>%
slice_min(iso_time) %>%
ungroup() %>%
left_join(
counties_4326_update
) %>%
st_sf(crs=4326)
triggered <- triggered_direct %>%
unnest("data")%>%
mutate(
trigger = if_else(
geoid == geoid_neighbor,
"direct",
"indirect"
),
county_name = ifelse(
trigger == "indirect",
NA,
county_name
),
geoid = if_else(
geoid == geoid_neighbor,
geoid,
geoid_neighbor
)
) %>%
select(
-geoid_neighbor,
-county_name
) %>%
left_join(
counties_4326_update %>%
tibble() %>%
select(
county_name,
geoid,
county_timezone
)
) %>%
mutate(
local_time = iso_time %>% with_tz(county_timezone)
) %>%
select(
-iso_time
) %>%
group_by(
usa_atcfid,
geoid,
storm_name,
county_name,
county_timezone
) %>%
summarize(
start_local_time = min(local_time),
end_local_time = max(local_time),
) %>%
ungroup() %>%
mutate(
trigger_interval_local_time = interval(start_local_time, end_local_time)
)
triggered_4326 <- triggered %>%
left_join(
counties_4326 %>%
tibble() %>%
select(
geoid,
geometry
)
) %>%
filter(
geoid != 12087.2
) %>%
st_sf(crs=4326)
# Map
north_america_4326 <- ne_countries(
continent = 'north america',
returnclass = "sf"
) %>%
st_transform(crs=4326)
ggplot() +
geom_sf(
data = north_america_4326,
fill = "#f7f7f7"
) +
geom_sf(
data = counties_4326,
fill = "#f7f7f7"
) +
geom_sf(
data = triggered_4326,
fill = "#67a9cf"
) +
geom_sf(
data = hull_5070 %>%
st_transform(crs=4326) %>%
st_union() %>%
st_crop(
xmin = triggered_4326 %>%
st_bbox() %>%
pluck(1) - .3,
ymin = triggered_4326 %>%
st_bbox() %>%
pluck(2) - .3,
xmax = triggered_4326 %>%
st_bbox() %>%
pluck(3) + .3,
ymax = triggered_4326 %>%
st_bbox() %>%
pluck(4) + .3,
),
color = "#ef8a62",
alpha = .5,
fill = "#ef8a62",
size = 1,
show.legend = FALSE
) +
geom_sf(
data = triggered_direct_4326_t,
fill = NA
) +
geom_sf(
data = hull_5070 %>%
st_transform(crs=4326) %>%
st_union() %>%
st_crop(
xmin = triggered_4326 %>%
st_bbox() %>%
pluck(1) - .3,
ymin = triggered_4326 %>%
st_bbox() %>%
pluck(2) - .3,
xmax = triggered_4326 %>%
st_bbox() %>%
pluck(3) + .3,
ymax = triggered_4326 %>%
st_bbox() %>%
pluck(4) + .3,
),
color = "#ef8a62",
fill = NA,
size = 1,
show.legend = FALSE
) +
coord_sf(
xlim =
c(
triggered_4326 %>%
st_bbox() %>%
pluck(1),
triggered_4326 %>%
st_bbox() %>%
pluck(3)
),
ylim =
c(
triggered_4326 %>%
st_bbox() %>%
pluck(2),
triggered_4326 %>%
st_bbox() %>%
pluck(4)
)
) +
theme(
panel.background = element_rect(fill = "aliceblue")
) +
ggtitle(
paste0(
"Hurricane ",
hull_5070$storm_name %>% first()
),
subtitle = paste0("HIP-WI County Triggers")
)
## Triggered Counties
triggered_4326 %>%
tibble() %>%
mutate(
county_name = county_name %>% str_to_title()
) %>%
select(
county_name,
geoid,
storm_name,
trigger_interval_local_time
) %>%
kbl(
col.names = c("County", "FIP", "Storm Name", "Trigger Interval Local Time")
) %>%
kable_styling(
bootstrap_options = c(
"condensed",
"responsive"
)
)
si <- sessionInfo()
print(si,loadedOnly=T, locale = FALSE)