#load libraries
library(tidyverse)
library(tmap)
library(sf)
library(tigris)

#bolin creek daily data from 2015-2025
bolin_creek_data <-
  read_csv("https://drive.google.com/uc?export=download&id=1eEzriJ5XzR-aS74nJmyiQFVrmV4T3cel")

#NC asos station data from 2000-2020
asos_jan_data <-
  read_csv("https://drive.google.com/uc?export=download&id=1_zE1kcC7YY083R6rv20sqeN1l2dCdmhi")

#Get daily max
max_bolin_creek <- bolin_creek_data |> group_by(time) |> summarise(max_water = max(value))

Executive Summary

This briefing addresses two hazard-preparedness questions for NCEM: (1) whether repeated flood-risk days at Bolin Creek are common enough to justify an enhanced monitoring protocol, and (2) where in North Carolina the state’s next warming shelter should be located. In short: the raw probability of a bad 25-day flood stretch is low, but the historical record shows these events cluster around named storms far more than chance alone would predict — so we recommend an event-triggered monitoring protocol rather than a constant one. For cold-weather risk, we recommend Guilford County (Piedmont region) for the new shelter, based on meaningful freeze risk combined with a large and growing unsheltered population.

Section 1: Bolin Creek Flood Risk

Defining flood risk and how often it happens

A “flood risk day” is any day where the water height at Bolin Creek exceeds the median daily height for the 2015–2025 period of record by more than 1.5 feet. We calculate that threshold and the resulting empirical probability below.

flood_median <- median(max_bolin_creek$max_water)
flood_threshold <- flood_median + 1.5

max_bolin_creek <- max_bolin_creek |>
  mutate(flood_risk = max_water > flood_threshold)

n_flood_days <- sum(max_bolin_creek$flood_risk)
n_total_days <- nrow(max_bolin_creek)
p_flood <- n_flood_days / n_total_days

tibble(
  `Median daily height (ft)` = round(flood_median, 2),
  `Flood-risk threshold (ft)` = round(flood_threshold, 2),
  `Flood-risk days observed` = n_flood_days,
  `Total days on record` = n_total_days,
  `Empirical P(flood-risk day)` = round(p_flood, 4)
) |> knitr::kable()
Median daily height (ft) Flood-risk threshold (ft) Flood-risk days observed Total days on record Empirical P(flood-risk day)
1.83 3.33 118 3628 0.0325

Over eleven years of record, Bolin Creek crossed the flood-risk threshold on 118 of 3628 days — an empirical probability of 3.3% on any given day. The chart below shows every day of record, with flood-risk days flagged in red.

ggplot(max_bolin_creek, aes(x = time, y = max_water)) +
  geom_line(color = "grey45", linewidth = 0.3) +
  geom_point(data = filter(max_bolin_creek, flood_risk), color = "firebrick", size = 1.3) +
  geom_hline(yintercept = flood_threshold, color = "firebrick", linetype = "dashed") +
  annotate("text", x = min(max_bolin_creek$time), y = flood_threshold + 0.5,
           label = paste0("Flood-risk threshold: ", round(flood_threshold, 2), " ft"),
           hjust = 0, color = "firebrick", size = 3.3) +
  labs(title = "Daily Maximum Water Height at Bolin Creek (2015\u20132025)",
       subtitle = paste0(n_flood_days, " flood-risk days out of ", n_total_days, " days on record"),
       x = NULL, y = "Daily max water height (ft)") +
  theme_minimal(base_size = 12)

Notice that flood-risk days are not spread evenly across the eleven years — they arrive in tight clusters (visible spikes around 2016, 2018, and 2019), almost certainly tied to specific tropical storm systems moving through the watershed, rather than being scattered randomly day to day. This pattern matters a great deal for the next question.

How many flood-risk days should we expect in a 25-day window?

Using the empirical probability above as a per-day probability, we can model the number of flood-risk days in a 25-day period as a binomial random variable with \(n = 25\) and \(p = 0.0325\).

n_period <- 25
p_hat <- p_flood

binom_probs <- dbinom(0:n_period, size = n_period, prob = p_hat)
most_likely_days <- which.max(binom_probs) - 1
p_more_than_3 <- 1 - pbinom(3, size = n_period, prob = p_hat)
expected_days <- n_period * p_hat

tibble(
  `Expected (mean) flood-risk days in 25 days` = round(expected_days, 2),
  `Most likely number of days (mode)` = most_likely_days,
  `P(more than 3 flood-risk days in 25 days)` = round(p_more_than_3, 4)
) |> knitr::kable()
Expected (mean) flood-risk days in 25 days Most likely number of days (mode) P(more than 3 flood-risk days in 25 days)
0.81 0 0.0082

Under this model, the single most likely outcome for any 25-day stretch is 0 flood-risk days, and the probability of seeing more than 3 flood-risk days in a 25-day window is only 0.8%.

binom_df <- tibble(k = 0:n_period, prob = binom_probs) |>
  mutate(group = ifelse(k > 3, "More than 3 days", "3 or fewer days")) |>
  filter(prob > 0.0002)

ggplot(binom_df, aes(x = factor(k), y = prob, fill = group)) +
  geom_col() +
  scale_fill_manual(values = c("More than 3 days" = "firebrick", "3 or fewer days" = "steelblue")) +
  labs(title = "Binomial Model: Flood-Risk Days in a 25-Day Period",
       subtitle = paste0("P(X > 3) = ", round(p_more_than_3, 4), " under an independence assumption"),
       x = "Number of flood-risk days (k)", y = "Probability", fill = NULL) +
  theme_minimal(base_size = 12) + theme(legend.position = "top")

Does the binomial model tell the whole story?

Taken at face value, a 0.8% chance of a rough 25-day stretch would not obviously justify a standing enhanced-monitoring protocol. But that binomial model assumes each day’s flood risk is independent of the day before — and the time series above already suggests that’s not realistic for a creek whose worst days are storm-driven. We can check this directly by sliding a real 25-day window across the actual historical record, day by day, and counting how many flood-risk days fall inside each window:

ordered_days <- max_bolin_creek |> arrange(time)
flood_vec <- as.integer(ordered_days$flood_risk)
n_days <- length(flood_vec)
window_counts <- sapply(1:(n_days - 24), function(i) sum(flood_vec[i:(i + 24)]))

empirical_p_more_than_3 <- mean(window_counts > 3)
max_observed_in_window <- max(window_counts)

tibble(
  `25-day windows checked` = length(window_counts),
  `Max flood-risk days seen in any real window` = max_observed_in_window,
  `Share of real 25-day windows with more than 3 days` = round(empirical_p_more_than_3, 4)
) |> knitr::kable()
25-day windows checked Max flood-risk days seen in any real window Share of real 25-day windows with more than 3 days
3604 8 0.0794

This is the key finding: in the actual historical record, 7.9% of all realized 25-day periods had more than 3 flood-risk days — roughly 10 times higher than the independence-based binomial estimate. One stretch in 2016 saw 8 flood-risk days within a single 25-day window. In other words, the low binomial probability understates real risk because it doesn’t account for storm-driven clustering.

Recommendation: A continuous, always-on enhanced monitoring posture is hard to justify statistically — the baseline day-to-day risk is genuinely low (3.3% per day), and most 25-day periods will see zero flood-risk days. However, given how strongly flood-risk days cluster around specific storm systems, we recommend NCEM adopt an event-triggered enhanced monitoring protocol: elevate monitoring specifically when a tropical system, remnant storm, or heavy rain event is forecast to track through the Bolin Creek watershed, rather than maintaining elevated monitoring year-round. This targets resources at the periods when clustered, multi-day flood risk is actually likely to materialize.

Section 2: Statewide Cold-Temperature Risk & Warming Shelter Siting

Estimating freeze-day probability at each station

For each of the state’s ASOS stations, we treat January daily maximum temperatures (2000–2019) as approximately normally distributed and use each station’s mean and standard deviation to estimate the probability that a January day’s maximum temperature falls below 32°F.

station_stats <- asos_jan_data |>
  group_by(station) |>
  summarise(
    n_obs = n(),
    mean_temp = mean(max_temp),
    sd_temp = sd(max_temp),
    lat = first(lat),
    lon = first(lon),
    .groups = "drop"
  ) |>
  mutate(p_below_32 = pnorm(32, mean = mean_temp, sd = sd_temp))

station_stats |>
  arrange(desc(p_below_32)) |>
  select(station, n_obs, mean_temp, sd_temp, p_below_32) |>
  mutate(across(c(mean_temp, sd_temp), \(x) round(x, 1)),
         p_below_32 = scales::percent(p_below_32, accuracy = 0.1)) |>
  head(10) |>
  knitr::kable(caption = "Ten highest-risk stations by estimated freeze-day probability")
Ten highest-risk stations by estimated freeze-day probability
station n_obs mean_temp sd_temp p_below_32
GEV 608 42.7 11.4 17.4%
AVL 620 47.7 10.9 7.5%
ILM 36 53.2 14.0 6.5%
GSO 620 49.0 11.0 6.2%
INT 620 49.0 11.0 6.1%
ASJ 612 51.4 12.4 5.8%
MWK 618 49.0 10.6 5.4%
UKF 597 49.5 10.9 5.4%
BUY 619 50.3 11.3 5.3%
RWI 616 51.3 11.9 5.3%

A few of these stations (n < 150 observations) have noticeably thinner records than the typical ~600-observation station and should be treated cautiously — a short record can produce a noisy mean/SD and an unreliable probability estimate, even if the number itself looks high.

Mapping risk across the state

To visualize risk statewide, we assign each North Carolina county the estimated freeze-day probability of its nearest ASOS station and map the result as a classed choropleth, with station locations overlaid as points.

# Try to get official Census county boundaries via tigris; if unavailable
# (e.g., no internet access at knit time), fall back to base R map boundaries
# so the report still renders end-to-end.
nc_counties <- tryCatch({
  options(tigris_use_cache = TRUE)
  tigris::counties(state = "NC", cb = TRUE, year = 2020, progress_bar = FALSE) |>
    sf::st_transform(4326) |>
    dplyr::rename(county_name = NAME) |>
    dplyr::select(county_name)
}, error = function(e) {
  message("tigris/Census download unavailable this session; using maps-package boundaries instead.")
  sf::st_as_sf(maps::map("county", "north carolina", plot = FALSE, fill = TRUE)) |>
    dplyr::mutate(county_name = stringr::str_to_title(stringr::str_remove(ID, "north carolina,"))) |>
    dplyr::select(county_name) |>
    sf::st_set_crs(4326) |>
    sf::st_make_valid()
})

stations_sf <- sf::st_as_sf(station_stats, coords = c("lon", "lat"), crs = 4326, remove = FALSE)

sf::sf_use_s2(FALSE)
county_centroids <- sf::st_centroid(nc_counties)
nearest_idx <- sf::st_nearest_feature(county_centroids, stations_sf)
nc_counties$nearest_station <- stations_sf$station[nearest_idx]
nc_counties$p_below_32 <- stations_sf$p_below_32[nearest_idx]
breaks <- classInt::classIntervals(nc_counties$p_below_32, n = 5, style = "jenks")$brks

use_tmap <- requireNamespace("tmap", quietly = TRUE) &&
  requireNamespace("tmaptools", quietly = TRUE)

if (use_tmap) {
  tmap::tmap_mode("plot")
  tmap::tm_shape(nc_counties) +
    tmap::tm_polygons("p_below_32", breaks = breaks, palette = "Blues",
                       title = "P(max temp < 32\u00b0F)", border.col = "white") +
    tmap::tm_shape(stations_sf) +
    tmap::tm_dots(size = 0.1, col = "black") +
    tmap::tm_layout(title = "Probability of a January Day Below Freezing, by County",
                     legend.outside = TRUE) +
    tmap::tm_compass(position = c("left", "bottom")) +
    tmap::tm_scalebar(position = c("left", "bottom"))
} else {
  nc_counties$prob_class <- cut(nc_counties$p_below_32, breaks = breaks, include.lowest = TRUE,
                                 labels = paste0(round(breaks[-6] * 100, 1), "\u2013",
                                                  round(breaks[-1] * 100, 1), "%"))
  ggplot(nc_counties) +
    geom_sf(aes(fill = prob_class), color = "white", linewidth = 0.2) +
    geom_sf(data = stations_sf, shape = 21, fill = "black", color = "white", size = 1.6) +
    scale_fill_brewer(palette = "Blues", name = "P(max temp\n< 32\u00b0F)") +
    labs(title = "Probability of a January Day Below Freezing, by County",
         subtitle = "County value assigned from nearest ASOS station, NC, 2000\u20132019",
         caption = "Black dots = ASOS station locations") +
    theme_minimal(base_size = 12) +
    theme(axis.text = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank())
}

The pattern is exactly what elevation would predict: the western mountains (Ashe, Buncombe, and neighboring counties) show the highest freeze-day probabilities in the state, in some cases exceeding 15–17%. The Piedmont forms a second, moderate-risk band (roughly 5–6%), and the Coastal Plain is consistently the lowest-risk region (commonly 2–3%). Because the western mountains already have relatively strong cold-weather infrastructure in place, we look to the Piedmont for the new shelter site.

Warming shelter recommendation: Guilford County

station_stats |>
  filter(station == "GSO") |>
  mutate(p_below_32 = scales::percent(p_below_32, accuracy = 0.1)) |>
  select(station, n_obs, mean_temp, sd_temp, p_below_32) |>
  knitr::kable(caption = "Piedmont Triad International Airport (GSO), Guilford County")
Piedmont Triad International Airport (GSO), Guilford County
station n_obs mean_temp sd_temp p_below_32
GSO 620 48.96855 11.04478 6.2%

We recommend Guilford County (served by the Piedmont Triad International Airport station, GSO) for NCEM’s new warming shelter, for three reasons:

  1. Meaningful, well-supported cold risk. Guilford’s estimated freeze-day probability is about 6.2%, well above the statewide Piedmont/Coastal Plain baseline, and it’s backed by a full 20-year, ~620-day station record (not one of the thin, less reliable records flagged above).
  2. Not in the west. Guilford is squarely in the Piedmont, so this recommendation doesn’t duplicate cold-weather infrastructure the western mountain counties already have.
  3. A large and growing vulnerable population. Guilford County’s annual Point-in-Time count of people experiencing homelessness rose from 452 in 2023 to 641 in 2024, and the county’s poverty rate (13.4%) runs above both the state (12.6%) and national (12.2%) rates. A meaningfully cold Piedmont county with a fast-growing unsheltered population is exactly the profile a new warming shelter should target.

As a secondary option, Winston-Salem/Forsyth County (station INT, ~6.1% freeze probability) shows a very similar risk profile and sits in the same region, and could be considered if siting logistics in Guilford prove difficult.

Summary of Recommendations

Question Recommendation
Enhanced flood monitoring at Bolin Creek? Don’t run it continuously — baseline day-to-day risk is low. Do trigger it around forecast tropical/heavy-rain events, since flood-risk days cluster strongly around storms (historical 25-day windows exceeded the binomial estimate by roughly 10x).
Where should the new warming shelter go? Guilford County — meaningful, well-supported freeze risk outside the already-prepared western mountains, paired with a large and rapidly growing unsheltered population.