Consider manually measured GW sites

[lindsayplatt]

One of the review comments (see #68) mentioned that if we are only using continuous data, we "may be missing most of the data". We wanted to explore this idea and we are certainly missing some of the data, but it really depends on the state. If you compare this data to our current map, the best gap the manually measured GW sites fill is in Utah & Louisiana. If we do add manual measurements to our data, we will need to consider how that integrates with continuous data (per another reviewers comments, see #68 this would skew statistics to compare discrete with continuous)

# Inventory and understand the manually measured GW sites
# Only query for sites that had a measurement in the last year

library(dataRetrieval)
library(tidyverse)
library(sf)

hucs <- zeroPad(1:21, 2) # all hucs

sites <- data <- tibble()
for(huc_i in hucs) {
  gwl_sites_i <- whatNWISsites(huc = huc_i, service = "gw", 
                               # Had data within the last year ("52 weeks")
                               period="P52W",
                               parameterCd = c("72019", "62610", "72150"))
  
  gwl_data_i <- whatNWISdata(siteNumber = unique(gwl_sites_i$site_no), 
                             service = "gw",
                             parameterCd = c("72019", "62610", "72150"))
  
  sites <- bind_rows(sites, gwl_sites_i)
  data <- bind_rows(data, gwl_data_i)
}

# Filter and plot the results

# Keep only sites with at least 120 measurements, spanning more than 10 years,
# which means an average of 1 measurement per month (12 per year for 10 years = 120 obs)
sites_from_gwl <- data %>% 
  filter(count_nu > 120 & end_date - begin_date > 10) %>% 
  pull(site_no) %>% 
  unique()

# Now map these sites
site_info <- readNWISsite(sites_from_gwl)
sites_sf <- site_info %>% st_as_sf(coords = c("dec_long_va", "dec_lat_va"))
us_sf <- st_as_sf(maps::map('state', plot=F, fill=T)) %>% st_transform(4326)

plot(st_geometry(us_sf))
plot(st_geometry(sites_sf), add=T)