How do remotely sensed snow cover metrics impact discharge in the same water

year in central Colorado?

Data checking

Data read in

First we need to get our snow metric (ndsi) data and then also download discharge data from the USGS

library(tidyverse) 
library(lubridate)
library(dataRetrieval) #for downloading USGS data


#ndsi 
ndsi <- read_csv('data/hayman_ndsi.csv') %>%
  rename(burned=2,unburned=3) %>%
  filter(!is.na(burned),
         !is.na(unburned)) %>%
  gather(.,key='site',
         value='ndsi',
         -DateTime) # For this analysis we want the data in long format
  
  
#USGS gauge above cheeseman lake '00060'
q_hayman <- readNWISdata(sites=c('06700000'), #Site code 
                  parameterCd='00060', #discharge code in cfs
                  service='dv', # service = daily values (versus annual)
                  startDate='1984-10-01', #Start date for getting the data
                  endDate = '2019-9-10') %>% # End date (today) 
  rename(q_cfs = X_00060_00003,
         quality_cd = X_00060_00003_cd) %>% #rename long column name
  filter(!is.na(q_cfs)) %>% #Drop NAs which can occur when there is ice or sensor breaks
  as_tibble() #To make it act like a tibble

Data exploring

NDSI summary

##     DateTime              site                ndsi        
##  Min.   :1984-04-10   Length:3208        Min.   :-0.5727  
##  1st Qu.:1999-10-13   Class :character   1st Qu.:-0.4835  
##  Median :2006-05-01   Mode  :character   Median :-0.4307  
##  Mean   :2005-06-30                      Mean   :-0.2364  
##  3rd Qu.:2013-03-17                      3rd Qu.:-0.1352  
##  Max.   :2019-08-02                      Max.   : 0.9459

Q summary

##   agency_cd           site_no             dateTime                  
##  Length:3132        Length:3132        Min.   :2002-08-01 00:00:00  
##  Class :character   Class :character   1st Qu.:2007-04-09 18:00:00  
##  Mode  :character   Mode  :character   Median :2011-05-30 12:00:00  
##                                        Mean   :2011-04-23 04:10:06  
##                                        3rd Qu.:2015-07-20 06:00:00  
##                                        Max.   :2019-09-09 00:00:00  
##      q_cfs       quality_cd           tz_cd          
##  Min.   :  53   Length:3132        Length:3132       
##  1st Qu.: 124   Class :character   Class :character  
##  Median : 179   Mode  :character   Mode  :character  
##  Mean   : 243                                        
##  3rd Qu.: 291                                        
##  Max.   :2210

Combining the data

Adding a water year column

When analyzing water flux data, we typically break the year up into “water years” which run from October to the end of September. For this exploratory analysis, we want to group the datasets by water year and then join them to each other so we can compare winter average, max, median, etc… of snow cover versus the next water year’s water flux. So we have to add a column called water year

Q water year

q_water_year <- q_hayman %>%
  mutate(month=month(dateTime),
         year_offset = ifelse(month > 9,1,0),
         wtr_yr = year(dateTime) + year_offset)

table(q_water_year$wtr_yr)
## 
## 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 
##   61  183  183  164  183  183  183  183  183  183  183  183  183  183  183 
## 2017 2018 2019 
##  183  183  162

NDSI water year

## 
## 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 
##   14   34   44   66   44   48   26   34   50   56   60   60   64   60   60 
## 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 
##   78  106  120  134  116  128  132  132  130  114  112  128  130   90  108 
## 2014 2015 2016 2017 2018 2019 
##  138  116  130  114  130  102

Filtering and summarizing

Now that we have our matched datasets we want to do a couple filtering operations. First, we want to make sure that we are only analyzing complete water years from the Q dataset. Second, we want to make sure we are only summarizing the snow data over months where snow cover is possible, which I would guess is between october and may. Once we have these filtering operations done, we want to summarize each dataset by water year so we can eventually join them and see if snow cover predicts Q.

Snow water year summary statistics

## # A tibble: 36 x 4
##    wtr_yr mean_ndsi max_ndsi median_ndsi
##     <dbl>     <dbl>    <dbl>       <dbl>
##  1   1984   0.483      0.605      0.483 
##  2   1985  -0.104      0.539     -0.123 
##  3   1986  -0.130      0.436     -0.203 
##  4   1987   0.0614     0.655      0.247 
##  5   1988  -0.309      0.616     -0.425 
##  6   1989  -0.398     -0.208     -0.404 
##  7   1990  -0.269      0.524     -0.428 
##  8   1991   0.0130     0.632     -0.0550
##  9   1992  -0.0982     0.592     -0.272 
## 10   1993  -0.00697    0.626     -0.0943
## # ... with 26 more rows

Whole Q water year summaries

## # A tibble: 15 x 4
##    wtr_yr mean_q max_q median_q
##     <dbl>  <dbl> <dbl>    <dbl>
##  1   2003   121.   350      102
##  2   2004   145.   296      140
##  3   2006   222.   395      239
##  4   2007   325.   494      350
##  5   2008   294.   585      262
##  6   2009   253.   667      179
##  7   2010   251.   622      204
##  8   2011   272.   783      179
##  9   2012   149.   278      130
## 10   2013   184.   332      163
## 11   2014   280.   542      248
## 12   2015   562.  2210      291
## 13   2016   204.   349      204
## 14   2017   205.   582      155
## 15   2018   131.   281      117

Plots of Snow Cover vs Q

Mean Snow vs Median Q

Max Snow vs. Median Q