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

SNOTEL Location

Analysis in R

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

setwd("/Users/bgatts/Desktop/WWU/ESCI597")
CRLA.SNOTEL=read.csv("CRLA_SNOTEL_slim_corrected.csv", header=T, stringsAsFactors=F)

Head and Tail of the data:

head(CRLA.SNOTEL)
tail(CRLA.SNOTEL)

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)))
}

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 each variable 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)

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)

Splines

Using splines with raw data:

spl.8 <- smooth.spline(CRLA.SNOTEL[6],spar=0.8)$y
spl.6 <- smooth.spline(CRLA.SNOTEL[6],spar=0.6)$y
spl.4 <- smooth.spline(CRLA.SNOTEL[6],spar=0.4)$y
spl.7 <- smooth.spline(CRLA.SNOTEL[6],spar=0.7)$y

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(CRLA.SNOTEL[[6]],type="l",xlab="Index",ylab="Average Temperature",col="grey")
lines(spl.8, col = "blue", lwd = 2)
axis(3);axis(4)

Or, perhaps we could convert things to zoo objects.

CRLA.Tavgspl.8 <- zoo(spl.8,order.by = foo)
CRLA.Tavgspl.6 <- zoo(spl.6,order.by = foo)
CRLA.Tavgspl.4 <- zoo(spl.4,order.by = foo)
CRLA.Tavgspl.7 <- zoo(spl.7,order.by = foo)

par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(CRLA.Tavgzoo,type="l",xlab="Year",ylab="Average Temperature",col="grey")
lines(CRLA.Tavgspl.8, col = "blue", lwd = 2)
lines(CRLA.Tavgspl.6, col = "purple", lwd = 2)
lines(CRLA.Tavgspl.4, col = "red", lwd = 2)
lines(CRLA.Tavgspl.7, col = "green", lwd = 2)
axis(4);axis(4)

PDO and ENSO Periodicities

load(file = "TimeSeries/ENSO_PDO.Rdata")
par(mfcol=c(2,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
require(dplR)
## Loading required package: dplR

creating a 24-month spline with the {dplR} funciton ffcsaps()

enso.spl24 <- ffcsaps(enso.monthly,nyrs=24)
pdo.spl24 <- ffcsaps(pdo.monthly,nyrs=24)

Next convert the spline to a time series object to plot the overlay nicely.

enso.spl24ts=ts(enso.spl24, start=1900, frequency=12)
pdo.spl24ts=ts(pdo.spl24, start=1900, frequency=12)

par(mfcol=c(2,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(enso.monthly,ylab="Nino 3.4 SST Index w/ spl24")
abline(h=0)
lines(enso.spl24ts, col="blue", lwd=2)
plot(pdo.monthly,ylab="PDO Index w/ spl24")
abline(h=0)
lines(pdo.spl24ts, type="l", col="blue", lwd=2)

Let’s zoom in to take a better look at the results using the plot par xlim=c(min,max)

par(mfcol=c(2,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(enso.monthly,ylab="Nino 3.4 SST Index w/ spl24", xlim=c(2000,2014))
abline(h=0)
lines(enso.spl24ts, col="blue", lwd=2)
plot(pdo.monthly,ylab="PDO Index w/ spl24", xlim=c(2000,2014))
abline(h=0)
lines(pdo.spl24ts, type="l", col="blue", lwd=2)

par(mfcol=c(2,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
spectrum(enso.monthly,log="no", main="ENSO")
spectrum(pdo.monthly,log="no", main="PDO")

We can see that the main periodicities are frequency <1 (period of greater than a year), so let’s zoom in.

par(mfcol=c(2,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
enso.spec= spectrum(enso.monthly,log="no", main="ENSO", xlim=c(0,1))
abline(v=0.25)
pdo.spec= spectrum(pdo.monthly,log="no", main="PDO", xlim=c(0,1))

Okay, what about dividing the series by the spline? I don’t really understand this plot.

enso.ratio=enso.monthly/enso.spl24
plot(enso.ratio)

How about subtracting it?

enso.diff=enso.monthly-enso.spl24
plot(enso.diff)

par(mfcol=c(2,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(enso.monthly,ylab="Nino 3.4 SST Index w/ spl24")
plot(enso.diff)

It looks like subtracting the 24-month spline leaves us with only the frequencies less than 24-months.

par(mfcol=c(2,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
spectrum(enso.monthly,log="no", main="ENSO")
spectrum(enso.diff,log="no", main="ENSO - 24-month spline")
abline(v=1/4)

Here we separate the three components: the original ENSO data, the 24-month spline, and the difference of the two.

par(mfcol=c(3,1),tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(enso.monthly)
plot(enso.spl24, type="l")
plot(enso.diff)

Here we zoom in fron 2000-2014 and we can really see the effect the subtracting the spline on the data.

plot(enso.monthly, xlim=c(2000,2014))
abline(h=0, lty=2)
lines(enso.spl24ts, col="blue", lwd=2)
lines(enso.diff, col="dark green", lwd=2)

What about a high frequency spline?

enso.spl4 <- ffcsaps(enso.monthly,nyrs=4)
enso.spl4ts=ts(enso.spl4, start=1900, frequency=12)
enso.diff4=enso.monthly-enso.spl4

plot(enso.monthly, xlim=c(2000,2014))
abline(h=0, lty=2)
lines(enso.spl4ts, col="blue", lwd=2)
lines(enso.diff4, col="dark green", lwd=2)

Leaves only very high frequency variability, might be a good way to filter out noise.

anyways..

Frequency Analysis for CRLA.SNOTEL data

CRLA.Tavgts <- ts(CRLA.SNOTEL[6],frequency = 365)
CRLA.Tavglowess=lowess(CRLA.Tavgts, f = 1/30)


spectrum(CRLA.Tavgzoo,log="no", main="CRLA Average Temperature")

We can see that by far the primary periodicity is about 1/365, or an annual cycle:

spectrum(CRLA.Tavglowess$y,log="no", main="CRLA Average Temperature Lowess Spline", xlim=c(0,.005))
abline(v=1/365)