lab1-pt2.R

# Part 2
# Loading Required Libraries
library(dataRetrieval)

## Warning: package 'dataRetrieval' was built under R version 4.3.3

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

#Adding site data
USGS_ID <- "11042000" 
DA <- 1443
gauge_lat <- 33.21809119
mmw_file <- "san_luis_rey_climate.csv"

#Adding function used for later
#Function for PET by Hamon method. Adapted from JP Gannnon, VT Hydroinformatics
hamonPET <- function(Month, Temp1, lat) {
  latrad <- (lat/360) * 2 * pi #convert to radians
  DOY <- yday(ymd(paste("2000", as.integer(Month), "15", sep= ' ')))
  tempvar = (2 * pi / 365) * DOY
  #declination of the sun above the celestial equator in radians on JulDay of the year
  delta_h = 0.4093 * sin(tempvar - 1.405)
  #day length in hours
  daylen = (2 * acos(-tan(delta_h) * tan(latrad)) / 0.2618)
  PET = 29.8 * daylen * 0.611 * exp(17.3 * Temp1 / (Temp1 + 237.3)) / (Temp1 + 273.2)
  return(PET) # PET in mm/day
}

#Downloading Monthly Stream flow
dfraw <- readNWISdata(
  site = c(11042000),
  parameterCd = "00060",
  service = "stat",
  statReportType = "monthly",
  startDate = "1988-10",
  endDate = "2024-04"
)

## Please be aware the NWIS data service feeding this function is in BETA.
## 
##             Data formatting could be changed at any time, and is not guaranteed

#Reformatting USGS Data
qm <- data.frame(date=as.Date(with(dfraw,paste(year_nu,month_nu,"15",sep="-")),"%Y-%m-%d"),
                 flow_cfs = dfraw$mean_va)

#Plotting Data to make sure it looks reasonable
plot(qm, type="l", xlab = NA, ylab = "Flow, cfs")

#Calculating monthly average values
qm_avg <- aggregate(list(flow = qm$flow_cfs),
                    by = list(month = format(qm$date, '%m')),
                    FUN = mean)

#Converting to cm/month
dpm <- days_in_month(as.integer(qm_avg$month))

#Calculating Average runoff in cm/month
qm_avg$runoff_cm <- qm_avg$flow * 86400 * dpm * (30.48^3) / (DA * 10^10)

#Importing data
pt <- read_csv(mmw_file)

## Rows: 13 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Month
## dbl (2): Mean Precip. (cm), Mean Temp. (°C)
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

pt <- rename(pt, prcp_cm = `Mean Precip. (cm)`, temp_c = `Mean Temp. (°C)`)

#Combing climate with flow data
monavg <- cbind(qm_avg, pt[1:12,])

#Coverts output from mm/d to cm/month
monavg$PET_cm <- hamonPET(monavg$month,monavg$temp_c,gauge_lat) * dpm / 10.0

#
runoff_sum <- sum(monavg$runoff_cm)
prcrp_sum <- (monavg$prcp_cm)

#Creating a plot
xx <- barplot(t(cbind(monavg$runoff_cm, monavg$PET_cm)),
              names = month.abb[as.numeric(qm_avg$month)],
              ylab = "Average runoff or precip, cm",
              cex.names = 0.7,
              col=c("steelblue","orange"),
              ylim = c(0, max(monavg$runoff_cm, monavg$prcp_cm))
)
lines(xx,monavg$prcp_cm, xaxt="n", col="darkblue",lty=1, lwd=2)
legend(x="top",legend = c("runoff", "PET", "precip"), fill = c("steelblue", "orange", NA),
       border = c("black", "black", NA), lty = c(NA, NA, 1),col = c(NA, NA, "darkblue"),
       lwd = c(NA, NA, 2), bty = "n", horiz = FALSE)
grid()