Introduction

This document demonstrates how to use the soilDB package to download climate data from the SCAN/SNOTEL network. Stay tuned for updates and more detailed examples.

Note

This interface is very much a work in progress. There are some SCAN / SNOTEL site with multiple (above-ground) sensors per sensor prefix, which can lead to confusing results. For now the first above-ground sensor is selected for each above-ground sensor prefix: c('TAVG', 'PRCP', 'PREC', 'SNWD', 'WTEQ', 'WDIRV', 'WSPDV', 'LRADT').

There are some sites with multiple below-ground sensors installed at the same depth. A message is printed when this happens. You can use the sensor.id column to access specific sensors. We are working on a solution.

Installation

With a recent version of R (>= 2.15), it is possible to get all of the packages that this tutorial depends on via:

# run these commands in the R console
install.packages('latticeExtra', dep=TRUE)
install.packages('plyr', dep=TRUE)
install.packages('rvest', dep=TRUE)
install.packages('httr', dep=TRUE)
install.packages('reshape2', dep=TRUE)
install.packages('soilDB', dep=TRUE)

You will also need the latest version of soilDB and aqp:

install.packages('devtools', dep=TRUE)
devtools::install_github("ncss-tech/aqp", dependencies=FALSE, upgrade_dependencies=FALSE)
devtools::install_github("ncss-tech/soilDB", dependencies=FALSE, upgrade_dependencies=FALSE)

Quick Example: getting sensor metadata

Each SCAN / SNOTEL station may have a different set of sensors.

# load required packages
library(aqp)
library(soilDB)
library(reshape2)
library(plyr)
library(latticeExtra)

# get basic sensor metadata for several SCAN/SNOTEL sites
m <- SCAN_sensor_metadata(site.code=c(2072,356,2148,2187))

# print a list of sensor types by station
dlply(m, 'site.code', function(i) unique(i$Ecode))
## $`356`
##  [1] "WTEQ" "PREC" "TOBS" "TMAX" "TMIN" "TAVG" "SNWD" "SMS"  "STO"  "SAL"  "RDC"  "BATT"
## 
## $`2072`
##  [1] "PREC"  "TOBS"  "TMAX"  "TMIN"  "TAVG"  "PRCP"  "SMS"   "STO"   "SAL"   "RDC"   "BATT"  "WDIRV"
## [13] "WSPDV" "RHUM"  "RHENC" "LRADT"
## 
## $`2148`
##  [1] "PREC"  "TOBS"  "TMAX"  "TMIN"  "TAVG"  "PRCP"  "SMS"   "STO"   "SAL"   "RDC"   "BATT"  "WDIRV"
## [13] "WSPDV" "RHUM"  "RHENC"
## 
## $`2187`
##  [1] "PREC"  "TOBS"  "TMAX"  "TMIN"  "TAVG"  "PRCP"  "SMS"   "STO"   "SAL"   "RDC"   "BATT"  "WDIRV"
## [13] "WSPDV" "RHUM"  "RHENC" "LRADT"
## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##   site.code
## 1       356
## 2      2072
## 3      2148
## 4      2187

Getting Data from a Single Station

Lets get some recent (daily) climate data from the Rogers Farm #1 SCAN station; site number 2001. The fetchSCAN() function does all of the work. Notice that the resulting object is a list, each element is a suite of sensor data. For example, the SMS element contains soil moisture from several depths

# fetchSCAN can accept multiple years
x <- fetchSCAN(site.code=2001, year=c(2013,2014,2015,2016))

# check the results
str(x, 1)
## List of 11
##  $ SMS  :'data.frame':   7325 obs. of  5 variables:
##  $ STO  :'data.frame':   7325 obs. of  5 variables:
##  $ SAL  :'data.frame':   0 obs. of  0 variables
##  $ TAVG :'data.frame':   1465 obs. of  5 variables:
##  $ PRCP :'data.frame':   1465 obs. of  5 variables:
##  $ PREC :'data.frame':   1465 obs. of  5 variables:
##  $ SNWD :'data.frame':   0 obs. of  0 variables
##  $ WTEQ :'data.frame':   0 obs. of  0 variables
##  $ WDIRV:'data.frame':   0 obs. of  0 variables
##  $ WSPDV:'data.frame':   0 obs. of  0 variables
##  $ LRADT:'data.frame':   0 obs. of  0 variables
# tabulate number of soil moisture measurements per depth (cm)
table(x$SMS$depth)
5 10 20 51 102
1465 1465 1465 1465 1465
# get unique set of soil moisture sensor depths
sensor.depths <- unique(x$SMS$depth)

# generate a better axis for dates
date.axis <- seq.Date(as.Date(min(x$SMS$Date)), as.Date(max(x$SMS$Date)), by='3 months')

# plot soil moisture time series, panels are depth
xyplot(value ~ Date | factor(depth), data=x$SMS, as.table=TRUE, type=c('l','g'), strip=strip.custom(bg=grey(0.80)), layout=c(1,length(sensor.depths)), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), ylab='Volumetric Water Content', main='Soil Moisture at Several Depths')

Getting Data from Multiple Stations

That was interesting, but most of the time we want to make comparisons between stations. Note that the format of data returned makes it possible "stack" data from several stations that do not share the same collection of sensors.

# get more data
x <- fetchSCAN(site.code=c(356, 2072), year=c(2015, 2016))

# same format as before, sensor data are "stacked" within each list element
str(x, 1)
## List of 11
##  $ SMS  :'data.frame':   5864 obs. of  5 variables:
##  $ STO  :'data.frame':   5864 obs. of  5 variables:
##  $ SAL  :'data.frame':   0 obs. of  0 variables
##  $ TAVG :'data.frame':   1466 obs. of  5 variables:
##  $ PRCP :'data.frame':   733 obs. of  5 variables:
##  $ PREC :'data.frame':   1466 obs. of  5 variables:
##  $ SNWD :'data.frame':   733 obs. of  5 variables:
##  $ WTEQ :'data.frame':   733 obs. of  5 variables:
##  $ WDIRV:'data.frame':   0 obs. of  0 variables
##  $ WSPDV:'data.frame':   0 obs. of  0 variables
##  $ LRADT:'data.frame':   0 obs. of  0 variables
# deeper look
str(x$SMS)
## 'data.frame':    5864 obs. of  5 variables:
##  $ Site     : int  356 356 356 356 356 356 356 356 356 356 ...
##  $ Date     : Date, format: "2015-01-01" "2015-01-02" "2015-01-03" ...
##  $ value    : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ depth    : num  5 5 5 5 5 5 5 5 5 5 ...
##  $ sensor.id: Factor w/ 5 levels "SMS.I_2","SMS.I_8",..: 1 1 1 1 1 1 1 1 1 1 ...
# get unique set of soil moisture sensor depths
sensor.depths <- unique(x$SMS$depth)

# generate a better axis for dates
date.axis <- seq.Date(as.Date(min(x$SMS$Date)), as.Date(max(x$SMS$Date)), by='3 months')

# plot soil moisture time series, panels are depth
xyplot(value ~ Date | factor(depth), groups=factor(Site), data=x$SMS, as.table=TRUE, type=c('l','g'), auto.key=list(columns=2, lines=TRUE, points=FALSE), strip=strip.custom(bg=grey(0.80)), layout=c(1,length(sensor.depths)), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), ylab='Volumetric Water Content', main='Soil Moisture at Several Depths')

Example Figures

# combine sensors
g <- make.groups('soil moisture'=x$SMS, 'soil temperature'=x$STO)

# get unique set of soil moisture sensor depths
sensor.depths <- unique(g$depth)

# generate a better axis for dates
date.axis <- seq.Date(as.Date(min(g$Date)), as.Date(max(g$Date)), by='3 months')

# plot soil moisture time series, panels are depth
p <- xyplot(value ~ Date | factor(Site) + which, groups=factor(depth), data=g, as.table=TRUE, type=c('l','g'), auto.key=list(columns=length(sensor.depths), lines=TRUE, points=FALSE), strip=strip.custom(bg=grey(0.80)), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y"), y=list(relation='free', rot=0)), ylab='', main='Soil Moisture and Temperature at Several Depths')

useOuterStrips(p)

# plot a single suite of sensors
# compare sites / depths
# no SAL sensors for these sites
# xyplot(value ~ Date | factor(depth), groups=factor(Site), data=x$SAL, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), layout=c(1,length(sensor.depths)), scales=list(alternating=3, y=list(relation='free', rot=0)))

xyplot(value ~ Date, groups=factor(Site), data=x$TAVG, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Average Air Temperature by Site')

xyplot(value ~ Date, groups=factor(Site), data=x$PRCP, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Daily Precipitation by Site')

xyplot(value ~ Date, groups=factor(Site), data=x$PREC, as.table=TRUE, type=c('l','g'), auto.key=list(title='site', columns=2, lines=TRUE, points=FALSE), scales=list(alternating=3, x=list(at=date.axis, format="%b\n%Y")), main='Accumulated Precipitation by Site')