Here is code to get just the TP numbers. I am tempted to alter your base R code into a dplyr pipe, but didn't.
#install.packages("dataRetrieval")
library("dataRetrieval")
#> Warning: package 'dataRetrieval' was built under R version 4.0.4
site_nums <- c("04087170", "04101500", "04108660", "04119400", "04121970")
params <- c("00665", "00060")
get_one_site <- function(site_num) {
TP <- readNWISqw(siteNumber = site_num, parameterCd = "00665")
TP <- TP[c("site_no", "sample_dt", "parm_cd", "result_va")]
names(TP)[names(TP) == "site_no"] <- "Site_Number"
names(TP)[names(TP) == "sample_dt"] <- "Sample_Date"
names(TP)[names(TP) == "parm_cd"] <- "Parameter_Code"
names(TP)[names(TP) == "result_va"] <- "mg_per_L"
return(TP)
}
combos_df <- tidyr::expand_grid(site_nums, params)
TPs <- purrr::map_dfr(combos_df$site_nums, get_one_site)
head(TPs)
#> Site_Number Sample_Date Parameter_Code mg_per_L
#> 1 04087170 1994-03-29 00665 0.08
#> 2 04087170 1994-04-08 00665 0.06
#> 3 04087170 1994-04-19 00665 0.07
#> 4 04087170 1994-05-24 00665 0.04
#> 5 04087170 1994-06-14 00665 0.04
#> 6 04087170 1994-07-08 00665 0.10
Created on 2021-03-25 by the reprex package (v0.3.0)
If you want to start pulling in parameters, you could use map2_dfr and pass the params column. I tried that but got an error (it appears that that "00060" discharge parameter is not reported for some sites, which threw and error in the dataRetrieval function that I don't have time to track down.