Crater Lake SNOTEL Site

The Annie Springs SNOTEL site at Crater Lake National Park has been in operation since September 2000 and provides daily measurements of Precipitation, Snow Depth, and Temperature variables. More information can be found on the USDA’s Natural Resource Conservation Service webpage here.

Snowy Crater Lake

[] (http://upload.wikimedia.org/wikipedia/commons/b/b5/Crater_lake_oregon.jpg) How can I make this image a link? Tried a few things and it didn’t work.

SNOTEL Location

Analysis in R

For this analysis we will be using several packages that are not native to R.

#install.packages("moments")
#install.packages("lubridate")
#install.packages("zoo")

First, we will read in the file from the *.csv

#setwd("U:/ESCI597/")
setwd("~/Desktop/WWU/ESCI597/TimeSeries")
CRLA.SNOTEL=read.csv("CRLA_SNOTEL_slim_corrected.csv", header=T, stringsAsFactors=F)

Take a look at the beginning and end of this record:

head(CRLA.SNOTEL)
##      Date SWE PPTacc Tmax Tmin Tavg PPTinc
## 1 10/1/00   0    0.0   57   40   50    0.1
## 2 10/2/00   0    0.1   59   31   43    0.0
## 3 10/3/00   0    0.1   62   28   43    0.0
## 4 10/4/00   0    0.1   61   34   46    0.0
## 5 10/5/00   0    0.1   62   29   45    0.0
## 6 10/6/00   0    0.1   65   35   48    0.1
tail(CRLA.SNOTEL)
##         Date  SWE PPTacc Tmax Tmin Tavg PPTinc
## 5291 3/27/15 11.7   46.8   61   32   45    0.0
## 5292 3/28/15 11.4   46.8   50   26   36    0.1
## 5293 3/29/15 11.3   46.9   58   29   42    0.0
## 5294 3/30/15 11.2   46.9   57   26   42    0.0
## 5295 3/31/15 11.0   46.9   35   26   31    0.5
## 5296  4/1/15 11.5   47.4   37   19   28    0.0

The variables of interest are:

names(CRLA.SNOTEL)
## [1] "Date"   "SWE"    "PPTacc" "Tmax"   "Tmin"   "Tavg"   "PPTinc"

which refer to the Date, SWE (Snow Water Equivalent - inches), PPTacc (accumulated precipitation for the snow year, which begins October 1 - inches), Tmax (Daily Maximum Temperature - F), Tmin (Daily minimum temperature - F), Tavg (daily average temperture - F), and PPTinc (daily precipitation increments - inches)

Summary Statistics

Due to the different types of measurement here, we exclude the accumulated precipitation because it abruptly ends on October 1 when it re-sets.

library(moments)

for (ctr in c(2,4:7)){

kt <-kurtosis(CRLA.SNOTEL[[ctr]])
sk <- skewness(CRLA.SNOTEL[[ctr]])
my.xlim <- range(CRLA.SNOTEL[[ctr]])
h<-hist(CRLA.SNOTEL[[ctr]], breaks=10, col="lightblue", xlab=names(CRLA.SNOTEL)[ctr],
        main="",xlim=my.xlim) 
xfit<-seq(min(CRLA.SNOTEL[[ctr]]),max(CRLA.SNOTEL[[ctr]]),length=100) 
yfit<-dnorm(xfit,mean=mean(CRLA.SNOTEL[[ctr]]),sd=sd(CRLA.SNOTEL[[ctr]])) 
yfit <- yfit*diff(h$mids[1:2])*length(CRLA.SNOTEL[[ctr]]) 
lines(xfit, yfit, col="darkblue", lwd=2)
boxplot(CRLA.SNOTEL[[ctr]], horizontal=TRUE, add = TRUE, outline=TRUE,  axes=FALSE, ylim=my.xlim, col = "lightgreen",  boxwex=100)
#text(x = 5, y=1000, labels = paste("Kurtosis=",round(kt,2)),pos = 4)
#text(x = 5, y=950, labels = paste("Skewness=",round(sk,2)),pos = 4)
legend("topleft", paste("Kurtosis=",round(kt,2),"
Skewness=",round(sk,2)))
}

We can learly see that the incremental precipitation is non-normal. Typically precipitatin data are log-normal, and these measurements may be transofrmed for further research.

Time Series

For our next few steps we will need the data in a Time Series format.

CRLA.SWEts <- ts(CRLA.SNOTEL[2],frequency = 1)
CRLA.PPTaccts <- ts(CRLA.SNOTEL[3],frequency = 1)
CRLA.Tmaxts <- ts(CRLA.SNOTEL[4],frequency = 1)
CRLA.Tmints <- ts(CRLA.SNOTEL[5],frequency = 1)
CRLA.Tavgts <- ts(CRLA.SNOTEL[6],frequency = 1)
CRLA.PPTincts <- ts(CRLA.SNOTEL[7],frequency = 1)

CRLA.time=time(CRLA.SWEts)

Decomposition

My data are failing to decompose, as is evidenced by the cycle(CRLA.Tavgts) function returning all “1”s. I will have to ask Andy about this. Apparently this was due to an import error, using the wrong frequency.

CRLA.Tavgts365 <- ts(CRLA.SNOTEL[6],frequency = 365, start=274/365) #Start date is 10/1/2000 or Day Of Year 274
plot(decompose(CRLA.Tavgts365))

Here time is in years from 10/1/2000.

Using the lubridate and zoo packages to match dates with readings

# lubridate - handles dates very nicely, we use the function parse_date_time() 
# to import the data (arg1) and the date format (arg2)
# install.packages("lubridate")

# zoo is Zotar's Ordered Objeczoo so we the create the CRLA.Tmax object in the zoo class form
# This makes it easy to plot the data with the appropriate labesl and handles leap years
# install.packages("zoo")

library(lubridate)
library(zoo)

We will create the lubridate object “foo” by reading in the snotel dates (which are of type string). This is basically turning the test dates in the .csv into a usable format for analysis. This also allows us to plot the dates along the x-axis easily. Technically this is a POSIXct date-time object. for more see http://www.inside-r.org/packages/cran/lubridate/docs/parse_date_time

foo <- parse_date_time(CRLA.SNOTEL$Date, "m/d/y")

Using the order.by() argument we assign these lubridate POSIXct objects to measurements for eachvariable of interest:

CRLA.Datezoo <- zoo(CRLA.SNOTEL[1],order.by = foo)
CRLA.SWEzoo <- zoo(CRLA.SNOTEL[2],order.by = foo)
CRLA.PPTacczoo <- zoo(CRLA.SNOTEL[3],order.by = foo)
CRLA.Tmaxzoo <- zoo(CRLA.SNOTEL[4],order.by = foo)
CRLA.Tminzoo <- zoo(CRLA.SNOTEL[5],order.by = foo)
CRLA.Tavgzoo <- zoo(CRLA.SNOTEL[6],order.by = foo)
CRLA.PPTinczoo <- zoo(CRLA.SNOTEL[7],order.by = foo)

This then allows us to plot the measurements nicely:

plot(CRLA.PPTacczoo)

plot(CRLA.SWEzoo)

plot(CRLA.Tmaxzoo)

plot(CRLA.Tminzoo)

plot(CRLA.Tavgzoo)

plot(CRLA.PPTinczoo)

Smoothing

Of course these plots are fairly noisy at this scale. To view them better, let’s add smoothed lines using the filter() function:

SWEma365 <- filter(x=CRLA.SWEzoo, filter=rep(x=1/365,times=365), sides=2)
SWEma30 <- filter(x=CRLA.SWEzoo, filter=rep(x=1/90,times=90), sides=2)
plot(CRLA.SWEzoo,col="grey")
lines(SWEma365,col="red",lwd=2)
lines(SWEma30,col="blue",lwd=2)

Tmaxma365 <- filter(x=CRLA.Tmaxzoo, filter=rep(x=1/365,times=365), sides=2)
Tmaxma30 <- filter(x=CRLA.Tmaxzoo, filter=rep(x=1/90,times=90), sides=2)
plot(CRLA.Tmaxzoo,col="grey")
lines(Tmaxma365,col="red",lwd=2)
lines(Tmaxma30,col="blue",lwd=2)

Tminma365 <- filter(x=CRLA.Tminzoo, filter=rep(x=1/365,times=365), sides=2)
Tminma30 <- filter(x=CRLA.Tminzoo, filter=rep(x=1/90,times=90), sides=2)
plot(CRLA.Tminzoo,col="grey")
lines(Tminma365,col="red",lwd=2)
lines(Tminma30,col="blue",lwd=2)

Tavgma365 <- filter(x=CRLA.Tavgzoo, filter=rep(x=1/365,times=365), sides=2)
Tavgma30 <- filter(x=CRLA.Tavgzoo, filter=rep(x=1/90,times=90), sides=2)
plot(CRLA.Tavgzoo,col="grey")
lines(Tavgma365,col="red",lwd=2)
lines(Tavgma30,col="blue",lwd=2)

PPTaccma365 <- filter(x=CRLA.PPTacczoo, filter=rep(x=1/365,times=365), sides=2)
PPTaccma30 <- filter(x=CRLA.PPTacczoo, filter=rep(x=1/90,times=90), sides=2)
plot(CRLA.PPTacczoo,col="grey")
lines(PPTaccma365,col="red",lwd=2)
lines(PPTaccma30,col="blue",lwd=2)

PPTincma365 <- filter(x=CRLA.PPTinczoo, filter=rep(x=1/365,times=365), sides=2)
PPTincma30 <- filter(x=CRLA.PPTinczoo, filter=rep(x=1/90,times=90), sides=2)
plot(CRLA.PPTinczoo,col="grey")
lines(PPTincma365,col="red",lwd=2)
lines(PPTincma30,col="blue",lwd=2)

Trend Lines

Trend lines in R are created using the linear model function lm()

CRLA.time=time(CRLA.Tavgts)

CRLA.TavgLM=lm(CRLA.Tavgts~CRLA.time)
plot(CRLA.Tavgts,col="grey")
# Heck, let's add the trend line from our linear model! 
abline(CRLA.TavgLM, col="black",lwd=2, lty="dashed") # So satisfying! - Thanks Andy

legend("topleft", paste("Slope (°F per year) =",round(CRLA.TavgLM$coefficients[[2]]*365, 3)))