lab1-pt2.R

title: “Lab 1” author: “John King” date: “April 5, 2024”

# Load necessary 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

## Warning: package 'ggplot2' was built under R version 4.3.2

## Warning: package 'tibble' was built under R version 4.3.2

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

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

## Warning: package 'purrr' was built under R version 4.3.2

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

## Warning: package 'stringr' was built under R version 4.3.2

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

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

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ 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

# Input site data
USGS_ID <- "11182500"
DA <- 15.26
gauge_lat <- 37.772983
mmw_file <- "crowcanyonroad.csv"

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 USGS data from website
dfraw <- readNWISdata(
  site = c(11182500),
  parameterCd = "00060",
  service = "stat",
  statReportType = "monthly",
  startDate = "1991-01",
  endDate = "2020-12"
)

## 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

#Reformattting 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)

#Plot data
plot(qm, type="l", xlab = NA, ylab = "Flow, cfs")

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

#Adjustment for months with a different total number of days
dpm <- days_in_month(as.integer(qm_avg$month))

qm_avg$runoff_cm <- qm_avg$flow * 86400 * dpm * (30.48^3) / (DA*(10^10))

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)`)

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

#Calculating potential evapotranspiration
monavg$PET_cm <- hamonPET(monavg$month,monavg$temp_c,gauge_lat) * dpm / 10.0

#Summing runoff and precipitation from monavg
runoff_sum <- sum(monavg$runoff_cm)
prcrp_sum <- sum(monavg$prcp_cm)

#Create plot with average monthly values
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()