snow_data_assignment.Rmd

---
title: "Snow Data Assignment: Web Scraping, Functions, and Iteration"
author: "Nathan Mueller, Devin Hunt"
date: "2-7-2022"
output: 
  html_document:
    theme: cerulean
    toc: TRUE
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
library(rvest)
library(tidyverse)
library(lubridate)
library(readxl)

```


# Simple web scraping

R can read html using either rvest, xml, or xml2 packages. Here we are going to navigate to the Center for Snow and Avalance Studies  [Website](https://snowstudies.org/archived-data/) and read a table in. This table contains links to data we want to programatically download for three sites. We don't know much about these sites, but they contain incredibly rich snow, temperature, and precip data. 


## Reading an html 

### Extract CSV links from webpage

```{r, message=F, warning=F}
site_url <- 'https://snowstudies.org/archived-data/'

#Read the web url
webpage <- read_html(site_url)

#See if we can extract tables and get the data that way
tables <- webpage %>%
  html_nodes('table') %>%
  magrittr::extract2(3) %>%
  html_table(fill = TRUE)
#That didn't work, so let's try a different approach

#Extract only weblinks and then the URLs!
links <- webpage %>%
  html_nodes('a') %>%
  .[grepl('24hr',.)] %>%
  html_attr('href')

```

## Data Download

### Download data in a for loop

```{r, message=F, warning=F}

#Grab only the name of the file by splitting out on forward slashes
splits <- str_split_fixed(links,'/',8)

#Keep only the 8th column
dataset <- splits[,8] 

#generate a file list for where the data goes
file_names <- paste0('data/',dataset)

for(i in 1:3){
  download.file(links[i],destfile=file_names[i])
}

downloaded <- file.exists(file_names)

evaluate <- !all(downloaded)

```


### Download data in a map

```{r, message=F, warning=F}

#Map version of the same for loop (downloading 3 files)
if(evaluate == T){
  map2(links[1:3],file_names[1:3],download.file)
}else{print('data already downloaded')}

```

## Data read-in 

### Read in just the snow data as a loop

```{r, message=F, warning=F}
#Pattern matching to only keep certain files
snow_files <- file_names %>%
  .[!grepl('SG_24',.)] %>%
  .[!grepl('PTSP',.)]

empty_data <- list()

snow_data <- for(i in 1:length(snow_files)){
  empty_data[[i]] <- read_csv(snow_files[i]) %>%
  select(Year,DOY,Sno_Height_M)
}

snow_data_full <- do.call('rbind',empty_data)

summary(snow_data_full)
```


### Read in the data as a map function

```{r, message=F, warning=F}

our_snow_reader <- function(file){
  name = str_split_fixed(file,'/',2)[,2] %>%
    gsub('_24hr.csv','',.)
  df <- read_csv(file) %>%
    select(Year,DOY,Sno_Height_M) %>%
    mutate(site = name)
}

snow_data_full <- map_dfr(snow_files,our_snow_reader)

summary(snow_data_full)
```


### Plot snow data

```{r, message=F, warning=F}
snow_yearly <- snow_data_full %>%
  group_by(Year,site) %>%
  summarize(mean_height = mean(Sno_Height_M,na.rm=T))

ggplot(snow_yearly,aes(x=Year,y=mean_height,color=site)) + 
  geom_point() +
  ggthemes::theme_few() + 
  ggthemes::scale_color_few()
```


# Assignment:

1. Extract the meteorological data URLs. Here we want you to use the `rvest` package to get the URLs for the `SASP forcing` and `SBSP_forcing` meteorological datasets.

```{r, message=F, warning=F}

links_forcing <- webpage %>%
  html_nodes('a') %>%
  .[grepl('forcing',.)] %>%
  html_attr('href')

```
2. Download the meteorological data. Use the `download_file` and `str_split_fixed` commands to download the data and save it in your data folder. You can use a for loop or a map function. 

```{r, message=F, warning=F}

#Grab only the name of the file by splitting out on forward slashes
splits2 <- str_split_fixed(links_forcing,'/',8)

#Keep only the 8th column
dataset2 <- splits2[,8] 

#generate a file list for where the data goes
file_names2 <- paste0('data/', dataset2)

for(i in 1:2){
  download.file(links_forcing[i], destfile=file_names2[i])
}

```
3. Write a custom function to read in the data and append a site column to the data. 

```{r, message=F, warning=F}

# this code grabs the variable names from the metadata pdf file
library(pdftools)
headers <- pdf_text('https://snowstudies.org/wp-content/uploads/2022/02/Serially-Complete-Metadata-text08.pdf') %>%
  readr::read_lines(.) %>%
  trimws(.) %>%
  str_split_fixed(.,'\\.',2) %>%
  .[,2] %>%
  .[1:26] %>%
  str_trim(side = "left")

forcing_snow_reader <- function(dfile){
  name = str_split_fixed(dfile,'/',2)[,2] %>%
    gsub('_Forcing_Data.txt','',.) %>% 
    gsub('SBB_','',.)
  df <- read_table(dfile, col_names = headers, na = "NA")
  
# When using mutate 
  # df <- mutate(site = name)
  # Continued to recieve the error: 
  # Error in UseMethod("mutate") : no applicable method for 'mutate' applied to an object of class "character"
  
# Solution:
  df %>% add_column(site = name)
}

```



4. Use the `map` function to read in both meteorological files. Display a summary of your tibble.
```{r warning = FALSE}

# Attempted renaming strategies:
  # names(df) <- headers
  # 
  # select(Year,DOY,Sno_Height_M)
  # rename_with(headers) %>% 

forcing_data_full <- map_dfr(file_names2, forcing_snow_reader)

summary(forcing_data_full)


```

5. Make a line plot of mean temp by year by site (using the `air temp [K]` variable). Is there anything suspicious in the plot? Adjust your filtering if needed.

```{r, message=F, warning=F}

mean_temp_year <- forcing_data_full %>% 
  group_by(year, site) %>% 
  summarize(mean_temp = mean(`air temp [K]`, na.rm = TRUE))

ggplot(mean_temp_year, aes(year, mean_temp, color = site)) + 
  geom_line() + labs(title = "Mean Temperature [K] By Year", x = "Year", y = "Mean Temp [Kelvin]")

```
The year 2003 produced abnormally low mean temperatures for both sites.

6. Write a function that makes line plots of monthly average temperature at each site for a given year. Use a for loop to make these plots for 2005 to 2010. Are monthly average temperatures at the Senator Beck Study Plot ever warmer than the Snow Angel Study Plot?
Hint: https://ggplot2.tidyverse.org/reference/print.ggplot.html

```{r, warning = FALSE, message = FALSE}

monthly_temp <- function(data, year_value){
  dataf <- data %>% filter(year == year_value)
  
  month_values <- dataf %>% group_by(month, site) %>% 
    summarize(mean_temp = mean(`air temp [K]`, na.rm = TRUE))
  
  plot_title <- paste("Mean Temperature [K] By Month in", year_value)
  
  p1 <- ggplot(month_values, aes(month, mean_temp, color = site)) + 
    geom_line() + labs(title = plot_title, x = "Month", y = "Mean Temp [Kelvin]")
  print(p1)
  
}


year_input <- c(2005, 2006, 2007, 2008, 2009, 2010)

for (i in 1:length(year_input)){
  monthly_temp_plots <- monthly_temp(forcing_data_full, year_input[i])
  
}

```
From the 6 graphs created, we cannot find a monthly temperature observation where SBSP exceeds the SASP site.


Bonus: Make a plot of average daily precipitation by day of year (averaged across all available years). Color each site.

```{r, message = FALSE}
# Create a date from given days
forcing_dates <- forcing_data_full %>% mutate(Date = make_date(year, month, day)) %>% 
  filter(!is.na(Date)) %>% 
  mutate(day_val = yday(Date))

daily_values <- forcing_dates %>% group_by(day_val, site) %>% 
  summarize(mean_precip = mean(`precip [kg m-2 s-1]`, na.rm = TRUE))


p2 <- ggplot(daily_values, aes(day_val, mean_precip, color = site)) + 
    geom_line() + labs(title = "Average Daily Precipitation by Day of Year", x = "Day", y = "Precip [kg m-2 s-1]")
  print(p2)

```


Bonus #2: Use a function and for loop to create yearly plots of precipitation by day of year. Color each site.
```{r, message = FALSE}

daily_precip <- function(precip_data, year_value){
  precip_dataf <- precip_data %>% filter(year == year_value)
  
  day_values <- precip_dataf %>% group_by(day_val, site) %>% 
    summarize(mean_year_precip = mean(`precip [kg m-2 s-1]`, na.rm = TRUE))
  
  plot_title <- paste("Mean Precipitation [kg m-2 s-1] By Date in", year_value)
  
  p3 <- ggplot(day_values, aes(day_val, mean_year_precip, color = site)) + 
    geom_line() + labs(title = plot_title, x = "Date", y = "Mean Temp [Kelvin]")
  print(p3)
  
}


year_input <- c(2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011)

for (i in 1:length(year_input)){
  daily_precip_plots <- daily_precip(forcing_dates, year_input[i])
  
}
```

<script>
<div class="tocify-extend-page" data-unique="tocify-extend-page" style="height: 0;"></div>
</script>