Task1_Update.Rmd

---
title: "Task 1 Report"
author: "Katie Willi and Caitlin Mothes"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    theme: paper
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE,
                      warning = FALSE,
                      error = FALSE,
                      message = FALSE)
                      #cache = TRUE)

source("setup.R")
```

This report aims to update NPS on our progress related to **Task 1 in our Statement of Work:**

*Compile available data and information describing the hydrology, climate, and water supplies of priority NPS-identified parks (e.g., DEVA, JOTR, and MOJA).*

-   *Examples: Hydrologic data may include stream and spring flow data, groundwater level data, and water quality data associated with specific sources of NPS water supplies or suitable analogs. This may also include information describing contributing basins (e.g., catchment areas, groundwater recharge areas, surface water/groundwater interaction, and output from water balance models). Climate data may include historical precipitation and temperature data and locations and monitored parameters of weather stations and SNOTEL sites. Water supply information may include types of supplies (i.e., surface water or groundwater), point-of-diversion locations, diversion amounts, NPS-observed supply challenges, etc.*

### Technical Approach

We have developed a suite of self-contained R functions (i.e., no local data downloads necessary) that pull relevant data sets for any specified park unit(s) in the contiguous US (CONUS) (see Table 1). This workflow allows for easy reproducibility across all parks and quickly pulls data directly from R, saving local storage space and reducing the need for complicated file management if the option to save the data locally is not selected. Moreover, all of our code and project tasks/timelines can be found in the ROSSyndicate GitHub [nps_watershed_vulnerability repository](https://github.com/rossyndicate/nps_water_vulnerability).

```{r}
readxl::read_xlsx("data/conus_functions.xlsx") %>%
  mutate(Source = paste("<a href=\"", Link, "\">", Data_Source,"</a>")) %>%
  dplyr::select(-c("Data_Source", "Link")) %>%
  DT::datatable(filter = "none",
                caption = 'Table 1: List of current functions developed to mine relevant CONUS-wide hydrologic and climatic data. Asterisks indicate data that is not publicly available yet.',
                rownames = FALSE,
                fillContainer = T,
                escape = FALSE,
                options = list(dom = 't',
                               pageLength = nrow(.),
                               scrollY = '300px'))
```

<br>

Because the end-goal of this project is to assess the water supply and climatic challenges of many park units, our primary focus has been to mine data sets that are available across CONUS so the bulk of our analyses can be reproducible and comparable across our assessments. Nationally available field sampling/instrumentation data sets that we are utilizing include the National Water Quality Portal, the National Water Information System, and NPS Aquarius. Modeled/interpolated data sets include gridded water balance models (Tercek et al., 2021), freshwater use projections (Warziniack et al., 2022), and groundwater recharge estimates (Wollock, 2003). Additional national-scale data sets include the National Land Cover Dataset's Land Use Change Index, infrastructure layers from the National Park Service, and the National Inventory of Dams.

Though our focus is centered around utilizing nationally-available data sets, we would lose key information by not exploring more localized sources of data, including state-run water monitoring programs and state water rights data. Therefore, we have also pulled relevant state-specific data sets associated with our test parks: Death Valley, Mojave Desert, and Joshua Tree National Parks (see Table 2).

```{r}
readxl::read_xlsx("data/park_specific_functions.xlsx") %>%
  mutate(Source = paste("<a href=\"", Link, "\">", Data_Source,"</a>")) %>%
  dplyr::select(-c("Data_Source", "Link")) %>% 
  DT::datatable(filter = "none",
                caption = "Table 2: List of current functions developed to mine park-/state-specific hydrologic data.",
                rownames = FALSE,
                fillContainer = TRUE,
                
                escape = FALSE,
                options = list(dom = 't',
                               scrollY = '300px'))
```

<br>

Data can be collected using two separate spatial extents (Figure 1):

1.  Data within a specified distance of the park boundaries

2.  Data contained within the parks' surface water contributing area.

Park unit boundaries are downloaded from the NPS-IRMA DataStore. For the purpose of this report we have been selecting data that is within 33.3 kilometers of each park unit. For surface water contributing areas, we use the R package `nhdplusTools` to select NHDPlus Version 2 features that are then used to delineate watersheds for all surface water features within each park unit.

![**Figure 1**: Maps showing the two different approaches used for pulling data. Here, we are downloading water right point-of-diversion data associated with Death Valley National Park's park boundary (left) and its surface water contributing area (right).](data/1_report/DEVA_example.png){fig-align="center"}

<br>

#### The following maps and tables represent just some of the data that we have currently pulled for each focus park.

```{r, error = FALSE, warning = FALSE, message = FALSE, include=FALSE}
# import park units
park <- c("BRCA","ZION")

park_boundary <- getParkBoundary(park = c("BRCA","ZION"))

park_buffer <- park_boundary %>% 
    sf::st_buffer(., dist = 0.3) # dist in dd; approximately 33.3 kilometers

# delineate surface wawter contributing areas per park
watershed <- getWatersheds(aoi = park_boundary, save = TRUE, path = 'data/1_report/watersheds_utah.RDS')

watershed <- readRDS('data/1_report/watershed_utah.RDS') %>%
  # dissolve for data pulling:
  group_by(REGION) %>%
  summarize()

watershed_split <- readRDS('data/1_report/watershed.RDS')

# put all of these areas of interest together for a single data pull

aoi_split <- bind_rows(watershed, park_buffer)

aoi <- aoi_split %>%
  summarize()
```

<br>

<br>

# Water Rights Data

### Point-of-diversion Locations

We pulled water rights point-of-diversion locations for the states of California and Nevada using the `getPODCalifornia()` and `getPODNevada()` functions.

```{r}
# Download water right PODs from Nevada database:
# nv_water_rights <- getPODNevada(aoi = aoi, save = T, path = 'data/1_report/nv_water_rights.RDS')
nv_water_rights <- readRDS('data/1_report/nv_water_rights.RDS') %>%
  st_join(select(aoi_split, UNIT_CODE), left = TRUE) %>%
  distinct(., .keep_all=TRUE)

# Download water right PODs from California database:
# ca_water_rights <- getPODNevada(aoi = aoi, save = T, path = 'data/1_report/ca_water_rights.RDS')
ca_water_rights <- readRDS('data/1_report/ca_water_rights.RDS') %>%
  st_join(select(aoi_split, UNIT_CODE), left = TRUE) %>%
  distinct(., .keep_all=TRUE)
```

```{r}
mapview(nv_water_rights, 
        col.regions = c("#56B4E9","#F0E442"), 
        zcol = 'PARK_RIGHT', 
        alpha.regions = 1, 
        alpha = 0.75, 
        cex = 4, 
        layer.name = "PODs") +
mapview(ca_water_rights,
        col.regions = c("#56B4E9","#F0E442"), 
        zcol = 'PARK_RIGHT', 
        alpha.regions = 1, 
        alpha = 0.75, 
        cex = 4,
        legend = FALSE) +
mapview(park_buffer,
        col.regions = "black",
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(watershed_split, 
        col.regions="#56B4E9", 
        alpha.regions = 0.33, 
        lwd = 0, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(park_boundary, 
        col.regions = "#74a089", 
        alpha.regions = 0, 
        lwd = 2, 
        popup = FALSE, 
        legend = F, 
        homebutton = FALSE,
        label = FALSE)
```

<br>

```{r}
ca_water_rights %>%
  sf::st_drop_geometry() %>%
  arrange((UNIT_CODE)) %>%
  # kableExtra::kable(format='html', caption = "California Water Rights") %>%
  # kableExtra::kable_styling(position = 'center') %>%
  # kableExtra::scroll_box(width = '900px', height = '500px')
  DT::datatable(filter = "none",
                caption = 'California Water Rights',
                rownames = FALSE,
                fillContainer = T,
                options = list(
                               scrollY = '300px')
  )
```

<br>

```{r}
nv_water_rights %>%
  sf::st_drop_geometry() %>%
  arrange((UNIT_CODE)) %>%
  DT::datatable(filter = "none",
                caption = 'Nevada Water Rights',
                rownames = FALSE,
                fillContainer = T,
                options = list(
                               scrollY = '300px')
  )
  # kableExtra::kable(format='html', caption = "Nevada Water Rights") %>%
  # kableExtra::kable_styling(position='center') %>%
  # kableExtra::scroll_box(width = '900px', height = '500px')

rm(ca_water_rights,nv_water_rights)
```

<br>

<br>

### Water rights dockets

We can retrieve metadata for park water rights dockets in NPS-IRMA DataStore with `inventoryDataStore()`, and download them with `getDataStore()`.

```{r}
dockets <- inventoryDataStore(park = park) %>%
  dplyr::filter(referenceType == "Docket") %>%
  dplyr::select(-citation)

dockets %>%
  #sf::st_drop_geometry() %>%
  DT::datatable(filter = "none",
                caption = 'Water Rights Dockets',
                rownames = FALSE,
                fillContainer = T,
                
                options = list(autoWidth = TRUE,
                               columnDefs = list(list(width = '250px', targets = c(8,9))),
                               scrollY = '300px')
  )
  # kableExtra::kable(format='html', caption = "Water Rights Dockets") %>%
  # kableExtra::kable_styling(position='center') %>%
  # kableExtra::scroll_box(width = '900px', height = '500px')
rm(dockets)
```

<br>

<br>

# Water and Weather Data

### [USGS flow data](https://waterdata.usgs.gov/nwis)

We are able to inventory NWIS locations, as well as the type of data being collected there, within an area of interest using `inventoryNWIS()`. We can download the raw sensor data using `getNWIS()`.

```{r}
# pull metadata on USGS locations with flow data
# nwis_all <- inventoryNWIS(aoi = aoi, save = T, path = 'data/1_report/nwis_all.RDS')
nwis_all <- readRDS('data/1_report/nwis_all.RDS') %>%
  mutate(site_type = ifelse(is.na(site_type), "Other", site_type)) %>%
  filter(flow_data == "Flow") %>%
  st_join(select(aoi_split, UNIT_CODE), left = TRUE) %>%
  distinct(., .keep_all=TRUE) %>%
  select(UNIT_CODE, site_no, site_name, site_type, parameter_type, date_range, n_obs)
```

```{r}
mapview(nwis_all, 
        col.regions = c("#E69F00",  # light orange
                        "#009E73",  # green
                        "#F0E442",  # yellow
                         "#56B4E9",  # light blue
                        "#0072B2",  # dark blue
                        "#D55E00"), # orange
        zcol = 'site_type',
        alpha.regions = 0.5,
        alpha = 0.4,
        cex = 4, 
        layer.name = "NWIS Flow Sites")+
mapview(park_buffer,
        col.regions = "black",
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(watershed_split, 
        col.regions="#56B4E9", 
        alpha.regions = 0.5, 
        lwd = 0, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(park_boundary, 
        col.regions = "#74a089", 
        alpha.regions = 0, 
        lwd = 2, 
        popup = FALSE, 
        legend = F, 
        homebutton = FALSE,
        label = FALSE)
```

<br>

```{r}
nwis_all %>%
  sf::st_drop_geometry() %>%
  arrange((UNIT_CODE)) %>%
  DT::datatable(filter = "none",
                caption = 'NWIS Flow Data',
                rownames = FALSE,
                fillContainer = T,
                options = list(
                               scrollY = '300px')
  )
  # kableExtra::kable(format='html', caption = "NWIS Flow Data") %>%
  # kableExtra::kable_styling(position='center') %>%
  # kableExtra::scroll_box(width='900px',height='500px')
rm(nwis_all)
```

<br>

<br>

### [NPS Aquarius data](https://irma.nps.gov/AQWebPortal/Data/Map)

We can download data within NPS-Aquarius for any park of interest using `getAquarius()`.

```{r}
# aquarius <- getAquarius(park = park, save = T, path = 'data/1_report/all_aquarius.RDS') 
aquarius <- readRDS('data/1_report/all_aquarius.RDS') %>%
  distinct(.keep_all=TRUE) %>%
  mutate(timestamp = lubridate::ymd_hms(timestamp),
         year = lubridate::year(timestamp)) %>%
  group_by(siteID, Location.Type, parameter, Long, Lat, EPSG) %>%
  summarize(n_obs = n(),
            range = ifelse(min(year) == max(year), 
                           paste0(min(year)),
                           paste0(min(year), "-", max(year)))) %>%
  mutate(flow_data = ifelse(grepl("velocity|WaterPressure|discharge|flow|height|level|stage|depthtowater", 
                                  parameter, ignore.case=T) &! grepl("sediment", parameter, ignore.case = T), "Flow", "Other")) %>%
  dplyr::filter(flow_data == "Flow") %>%
  select(siteID, Location.Type, parameter, n_obs, range, Lat, Long) %>%
  sf::st_as_sf(coords=c("Lat","Long"), crs=4326) %>%
  sf::st_transform(st_crs(aoi)) %>%
  st_join(select(aoi_split, UNIT_CODE), left = TRUE) %>%
  distinct(., .keep_all=TRUE)
```

```{r}
mapview(aquarius,
        zcol = 'Location.Type',
        col.regions = c("#E69F00", # light orange
                        "#56B4E9", # light blue
                        "#009E73", # green
                        "#F0E442"), # yellow
        layer.name = "Aquarius Sites") + 
mapview(park_buffer,
        col.regions = "black",
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(watershed_split, 
        col.regions="#56B4E9", 
        alpha.regions = 0.33, 
        lwd = 0, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(park_boundary, 
        col.regions = "#74a089", 
        alpha.regions = 0, 
        lwd = 2, 
        popup = FALSE, 
        legend = F, 
        homebutton = FALSE,
        label = FALSE)
```

<br>

```{r}
aquarius %>%
  sf::st_drop_geometry() %>%
  arrange((UNIT_CODE)) %>%
  DT::datatable(filter = "none",
                caption = 'Aquarius Continuous Data',
                rownames = FALSE,
                fillContainer = T,
                options = list(
                               scrollY = '300px')
  )
  # kableExtra::kable(format='html', caption = "Aquarius Continuous Data") %>%
  # kableExtra::kable_styling(position='center') %>%
  # kableExtra::scroll_box(width='900px',height='500px')
rm(aquarius)
```

<br>

<br>

### NOAA weather stations

We can download NOAA weather station data within a specified area of interest and selected time period (here, from 1979-2022) using `getNOAA()`.

```{r}
 # noaa_all <- getNOAA(aoi = aoi, startDate = "1979-10-01", endDate = "2022-09-30",
 #                     save = T, path = "data/1_report/all_noaa.RDS")
noaa_all <- readRDS('data/1_report/noaa_all.RDS') %>%
  mutate(year = lubridate::year(lubridate::as_date(date))) %>%
  group_by(id) %>%
  summarize(range_1980 = ifelse(min(year) == max(year), 
                         paste0(min(year)), 
                         paste0(min(year), "-", max(year)))) %>%
  st_join(select(aoi_split, UNIT_CODE), left = TRUE) %>%
  filter(!is.na(UNIT_CODE))
```

```{r}
mapview(noaa_all,
        col.regions = "#CC79A7",
        cex = 4,
        layer.name = "NOAA Weather Stations") + 
mapview(park_buffer,
        col.regions = "black",
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(watershed_split, 
        col.regions="#56B4E9", 
        alpha.regions = 0.33, 
        lwd = 0, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(park_boundary, 
        col.regions = "#74a089", 
        alpha.regions = 0, 
        lwd = 2, 
        popup = FALSE, 
        legend = F, 
        homebutton = FALSE,
        label = FALSE)
```

<br>

```{r}
noaa_all  %>%
  sf::st_drop_geometry() %>%
  arrange((UNIT_CODE)) %>%
   DT::datatable(filter = "none",
                caption = 'NOAA Weather Stations',
                rownames = FALSE,
                fillContainer = T,
                options = list(
                               scrollY = '300px')
   )
  # kableExtra::kable(format='html', caption = "NOAA Weather Stations") %>%
  # kableExtra::kable_styling(position='center') %>%
  # kableExtra::scroll_box(width='900px',height='500px')
rm(noaa_all)
```

<br>

<br>

# Gridded/Interpolated Data

### [Principal aquifers](https://water.usgs.gov/GIS/metadata/usgswrd/XML/aquifers_us.xml) in park units

We have explored the option of pulling data associated with *Principal Aquifers of the 48 Conterminous United States, Hawaii, Puerto Rico, and the U.S. Virgin Islands* ([USGS, 2003](https://water.usgs.gov/GIS/metadata/usgswrd/XML/aquifers_us.xml#stdorder)); namely, selecting only aquifers that intersect park boundaries. Because many of these aquifers extend across large swaths of the US, we embedded an option to clip aquifer extents by a chosen distance outside of the park boundary (here, clipping associated aquifer extents to only 33.3 km around the park boundary).

```{r, message=FALSE, warning=FALSE, error=FALSE}
aquifers <- getAquifers(aoi = park_boundary, dist = 0.3)

mapview(aquifers, 
        zcol = "OBJECTID_1",
        col.regions = wesanderson::wes_palette("Zissou1", type = "continuous"),
        layer.name= "Aquifers within Park Units",
        popup = TRUE,
        legend = F) + 
mapview(park_boundary, 
        col.regions = "black", 
        alpha.regions = 0, 
        lwd = 2, 
        popup = FALSE, 
        legend = F, 
        homebutton = FALSE,
        label = FALSE)
rm(aquifers)
```

<br>

<br>

### [Mean Annual Ground-Water Recharge Index](https://water.usgs.gov/GIS/metadata/usgswrd/XML/rech48grd.xml)

We can import the mean annual ground-water recharge index grid (Wollock, 2003) and clip to the area(s) of interest using `getGWRecharge()`.

```{r}
recharge <- getGWRecharge(aoi = aoi, dist = 0)

mapview(recharge[[1]], 
        col.regions = c("#F21A00","#E1AF00","#EBCC2A","#78B7C5","#3B9AB2"),
        layer.name= "Monthly Groundwater Recharge Index") + 
mapview(park_buffer,
        col.regions = "black",
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(watershed_split, 
        col.regions="black", 
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(park_boundary, 
        col.regions = "#74a089", 
        alpha.regions = 0, 
        lwd = 2, 
        popup = FALSE, 
        legend = F, 
        homebutton = FALSE,
        label = FALSE)
rm(recharge)
```

<br>

<br>

### Layers from the [NPS Gridded Water Balance Model](http://screenedcleanedsummaries.s3-website-us-west-2.amazonaws.com/)

There are \~1,700 rasters associated with this database. Here we are showing just one, though all (or any combination thereof) of the summary rasters can be downloaded with the `getWBM()` function. However, the workflows linked to customizing outputs for the [*What water balance product do you need?*](https://screenedcleanedsummaries.s3.us-west-2.amazonaws.com/which_water_balance.html) data set requires NPS VPN access.

```{r}
# getWBM(aoi = aoi, type = "PET")

pet_historical <- raster::raster('data/1_report/historical_pet_1980_1999.tif')

mapview(pet_historical, 
        col.regions = c("#F21A00","#E1AF00","#EBCC2A","#78B7C5","#3B9AB2"),
        layer.name= "PET (1980-1999)") + 
mapview(park_buffer,
        col.regions = "black",
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(watershed_split, 
        col.regions="black", 
        alpha.regions = 0, 
        lwd = 1, 
        popup = FALSE, 
        legend = FALSE, 
        homebutton = FALSE,
        label = FALSE) +
mapview(park_boundary, 
        col.regions = "#74a089", 
        alpha.regions = 0, 
        lwd = 2, 
        popup = FALSE, 
        legend = F, 
        homebutton = FALSE,
        label = FALSE)
rm(pet_historical)
```

<br>

<br>

# Miscellaneous Data

### Trends in park visitation

We can import park visitation data using `getVisitation()`.

```{r}
# download visitation data from NPS-Irma STATS
visitor_stats <- getVisitation(type = 'park', units = park, 
                               startYear = 1900, endYear = 2022) %>%
  group_by(UnitCode,Year) %>%
  summarize(Annual=sum(RecreationVisitors))
```

```{r}
plotly::ggplotly(ggplot() +
  geom_point(data = visitor_stats, aes(x = Year, y = Annual, color = UnitCode)) +
  geom_line(data = visitor_stats, aes(x = Year, y = Annual, color = UnitCode)) +
  scale_color_manual(values = c("#56B4E9",    # light blue
                                "#009E73",    # green
                                "#0072B2")) + # dark blue
  theme_bw() +
  labs(color="Park Unit") +
  ylab("Annual Visitors") +
  scale_y_continuous(labels = scales::comma))
rm(visitor_stats)
```