Flash floods seem to be increasing in many areas. This post will show how to download local USGS flow and precipitation data and generate a 3-panel chart of flow, gage height and precipitation.

There was a tragic flash flood incident in Berks County Pennsylvania on July 11-12, 2019 that caused the death of 3 people, including a pregnant woman and her young son. This newspaper article (link) provides some details on this tragic event.

The timing of this flash flooding was interesting to me because the woman and her son were trapped in their car by a small tributary to the Schuylkill River at 4:30 PM on Thursday, July 11 and flooding in Philadelphia started about 1:30 Am on Friday, July 12th. The tragic drowning occurred on the Manatawny Creek , about 35 miles Northwest of Philadelphia.

I plan to write a series of posts to document the rainfall patterns and timing of the rainfall and subsequent flooding.

This initial post will present show the flow hydrograph, cumulative precipitation and gage height at for the USGS’s Schuylkill river gaging station #01474500 at Fairmount Dam.

Follow-up posts will review the upstream USGS data as well as national weather service rainfall data.

Here is the 3 panel chart showing the hydrograph, precipitation data and gage height data for the period July 9 – 16, 2019.

Here’s the R script for analyzing and plotting csv file downloaded from USGS site.

#################################################################################################################

## get source data file link
link <- "C:\\Users\\Kelly O'Day\\Desktop\\Schuylkill_flood\\CSV_files\\Schuylkil_Fairmont_Dam_7_days_july_2019.csv"
#(link <- choose.files())
df <- read.csv(link, as.is=T)
df$dt <- mdy_hm(df$datetime)
dt_min <- df$dt[1]
dt_max <- df$dt[nrow(df)]

peak_flow <- max(df$cfs, na.rm=T)
peak_flow_comma <- format(peak_flow, big.mark=',')
peak_dt <- df$dt[which(df$cfs == peak_flow)]
from_to <- paste(month(dt_min),"/",day(dt_min), " to ", month(dt_max),"/",day(dt_max),"/",year(dt_max),sep="")

# Create 3 panel charts for flow, precipitation and gage height
par(mfrow = c(3,1),ps=10, pty="m",xaxs="r",yaxs="r")
par(las=1,oma=c(2.5,2,3,1), mar=c(2,5,3,0))
par(mgp=c(4,1,0)) # Number of margin lines for title, axis labels and axis line

## Flow chart
plot_title <- paste("Flow - cfs")
plot(df$dt, df$cfs, type="l", xlab="", ylab = "Flow -cfs",xlim=c(dt_min,dt_max), yaxt="n",
main =plot_title)
axis(side=2, at=axTicks(2),
labels=formatC(axTicks(2), format="d", big.mark=','))

points(peak_dt, peak_flow, pch=16, col="red")
spacer <- 4*60*60
note <- paste("Peak @ ", peak_flow_comma, " cfs \n", peak_dt)
abline(h = seq(10000,50000,10000), col = "grey")
abline(v= seq(dt_min,dt_max,by="day"), col="grey")
text(peak_dt+ spacer,peak_flow * 0.975, note, adj=0, cex=1)
################################################
# Sheet title annotation
mtext("Flood Conditions @ Fairmount Dam (July 9-16, 2019)", side = 3, line = 3, adj = 0.5, cex = 1.2)

##### precipitation data analyis n& Chart
df$cum_precip <- cumsum(df$inches)
tot_precip <- df$cum_precip[nrow(df)]
precip_st_row <- length(subset(df$cum_precip, df$cum_precip == 0))
precip_st_dt <- df$dt[precip_st_row]

precip_note <- paste0("Total Precip - ",tot_precip, " inches")
precip_subset <- subset(df, df$cum_precip == tot_precip)
precip_end_dt <- mdy_hm(precip_subset[1,1])

#precip_end_dt_t <- mdy_hm(precip_end_dt)
plot(df$dt, df$cum_precip, type="l", xlab="", ylab = "Precipitation -inches",xlim=c(dt_min,dt_max),
main = "Cummulative Precipitation - Inches")
points(precip_st_dt,0, pch=12, col = "blue")
points(precip_end_dt,tot_precip, pch=12, col="blue")
abline(v= seq(dt_min,dt_max,by="day"), col="grey")
text(dt_min,tot_precip, precip_note, adj=0)
precip_st_note <- paste0(" Precipitation Starts @ ",precip_st_dt)
dur <- precip_end_dt - precip_st_dt

precip_end_note <- paste0(" Precipitation ends @ ",precip_end_dt,"\n ",dur, " hours")
text(precip_end_dt,tot_precip * 0.9, precip_end_note, adj=0, cex = 1)
text(precip_st_dt, 0, precip_st_note, adj=0, cex = 1)

#### gage height chart
gage_act_df <- subset(df, df$gage >= 10)
gage_act_dt <- gage_act_df$dt[1]
gage_act_note <- paste0("Gage Action Level @ \n",gage_act_dt)
plot(df$dt, df$gage, type="l", xlab="", ylab = "Gage Height - Ft",xlim=c(dt_min,dt_max),
main ="Gage height - ft")
abline(h = 10, col="brown")
abline(h=11, col = "red", lwd =1)
abline(v= seq(dt_min,dt_max,by="day"), col="grey")
points(gage_act_dt,gage_act_df[1,3],pch=11, col = "black")
text(gage_act_dt - 1*4*3600, 10,gage_act_note, adj=1)
min_gage <- min(df$gage, na.rm=T)
max_gage <- max(df$gage, na.rm=T)
delta_gage <- max_gage - min_gage
gage_note <- paste("Max Gage @ ", max_gage, " ft\nGage Increase ", delta_gage, " ft")
text(dt_min, max_gage*0.97, gage_note, adj=0, cex=1)

flood_act_note <- "Flood Action stage @ 10-ft"
flood_stage_note <-"Flood stage @ 11-ft"
text(dt_max-60*60,10.25, flood_act_note, adj = 1, cex=0.95)
text(dt_max - 1*1*60*60,11.25, flood_stage_note, adj = 1, cex=0.95)

##########################################
# Sheet annotation - Footer Notes
mtext("K O'Day - 7/19/19", side = 1, line = 3, adj = 0, cex = 0.9)
mtext("Data Source: USGS: https://waterdata.usgs.gov/nwis/inventory/?site_no=01474500", side=1, line = 2, adj=0, cex = 0.9)

# png(file="C:\\Users\\Kelly O'Day\\Desktop\\Schuylkill_flood\\myplot.png", bg="white")

dev.off()