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)
install.packages('sharpshootR', dep=TRUE)

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

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

Sensor Metadata

Metadata for Select Stations

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

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

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

knitr::kable(m, row.names = FALSE)
Name Site State Network County Elevation_ft Latitude Longitude HUC climstanm upedonid pedlabsampnum
Blue Lakes 356 California SNOTEL Alpine 8067 38.60800 -119.92437 180400120101 NA NA NA
Deep Springs 2187 California SCAN Inyo 5399 37.37222 -117.97383 180902010102 Deep Springs S2012CA027002
Eros Data Center 2072 South Dakota SCAN Minnehaha 1602 43.73835 -96.61455 101702031402 Eros Data Center S2003SD099001 03N0688
Jordan Valley Cwma 2148 Idaho SCAN Owyhee 4508 42.94863 -117.01060 170501080405 NA NA NA
# result is a list
x <- fetchSCAN(site.code=c(574), year=c(2016))

# print a list of sensor types
names(x[-1])
##  [1] "SMS"   "STO"   "SAL"   "TAVG"  "TMIN"  "TMAX"  "PRCP"  "PREC"  "SNWD"  "WTEQ"  "WDIRV" "WSPDV"
## [13] "LRADT"

Metadata for all Stations

data('SCAN_SNOTEL_metadata', package = 'soilDB')

map('state', lwd=2)
title(main='SCAN / SNOTEL Sites')
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, subset=Network == 'SCAN', bg='DarkRed', pch=21)
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, subset=Network == 'SNOTEL', bg='RoyalBlue', pch=21)

legend('bottomleft', legend=c('SCAN', 'SNOTEL'), pch=21, pt.bg=c('DarkRed', 'RoyalBlue'), bty='n')

map('state', lwd=2)
title(main='SCAN / SNOTEL Sites')
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, bg='RoyalBlue', pch=21, cex=0.85)
points(Latitude ~ Longitude, data=SCAN_SNOTEL_metadata, subset=!is.na(pedlabsampnum), bg='DarkRed', pch=22)

legend('bottomleft', legend=c('All Stations', 'Linked to Lab Data'), pch=c(21, 22), pt.bg=c('RoyalBlue', 'DarkRed'), bty='n')

Sensor and Site Data for 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(2014,2015,2016,2017))

# check the results
str(x, 1)
## List of 14
##  $ metadata:'data.frame':    1 obs. of  12 variables:
##  $ SMS     :'data.frame':    6656 obs. of  7 variables:
##  $ STO     :'data.frame':    7296 obs. of  7 variables:
##  $ SAL     :'data.frame':    0 obs. of  0 variables
##  $ TAVG    :'data.frame':    1459 obs. of  7 variables:
##  $ TMIN    :'data.frame':    1459 obs. of  7 variables:
##  $ TMAX    :'data.frame':    1459 obs. of  7 variables:
##  $ PRCP    :'data.frame':    1319 obs. of  7 variables:
##  $ PREC    :'data.frame':    1323 obs. of  7 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 
## 1458 1458 1434 1429  877
# 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')

Visualizing sensor depths within the context of a soil profile

One way to get a visual diagram of where the SCAN soil moisture and temperature sensors occur within the context of the soil profile horizons is to use the linkage in the attached metadata to join the sensor depths to the profile horizonation.

# get more data
x <- fetchSCAN(site.code=c(574), year=c(2016))
s <- fetchKSSL(pedlabsampnum = na.omit(x$metadata$pedlabsampnum), returnMorphologicData = FALSE)

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

# slice the KSSL data associated with the sensor depths
s1 <- slice(s, sensor.depths ~ .)

# plot profile of lab data horizons then over plot with location of sensor depths within them
par(mar=c(0,0,3,1))
plotSPC(s, name='hzn_desgn', label='pedon_id', id.style='top', cex.names=0.75, width=0.05, x.idx.offset=0.1, axis.line.offset=-15, space=2)
title(paste('Sensor Depths for SCAN site', unique(x$SMS$Site), sep=" "), line=0.5, cex.main=1.5)

# over plot the sensor depths on the previous plot
plotSPC(s1, width=0.005, name='hzn_top', label='pedon_id', print.id=FALSE, id.style='top', cex.names=0.75, x.idx.offset=0.1, lwd=2, fill='blue', plot.depth.axis=FALSE, add=TRUE)

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 14
##  $ metadata:'data.frame':    2 obs. of  12 variables:
##  $ SMS     :'data.frame':    5752 obs. of  7 variables:
##  $ STO     :'data.frame':    5843 obs. of  7 variables:
##  $ SAL     :'data.frame':    0 obs. of  0 variables
##  $ TAVG    :'data.frame':    1460 obs. of  7 variables:
##  $ TMIN    :'data.frame':    1460 obs. of  7 variables:
##  $ TMAX    :'data.frame':    1460 obs. of  7 variables:
##  $ PRCP    :'data.frame':    581 obs. of  7 variables:
##  $ PREC    :'data.frame':    1316 obs. of  7 variables:
##  $ SNWD    :'data.frame':    731 obs. of  7 variables:
##  $ WTEQ    :'data.frame':    731 obs. of  7 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':    5752 obs. of  7 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" ...
##  $ water_year: num  2015 2015 2015 2015 2015 ...
##  $ water_day : int  93 94 95 96 97 98 99 100 101 102 ...
##  $ 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)