hackathon_presentation.Rmd

---
title: "`qaqcr` Demo"
author: "Celeste, Steve, John, Kristin, Kaelin, Jamie, Amanda, and Patrick"
output: html_notebook
---

# Overview

A suite of tools addressing common yet critical data curating tasks aimed at the particular needs of the LTER research and Information Managers community, and that also that applies to the QA/QC needs of the broader research community. These tools are designed to automated quality control of sensor data, but may also be useful for non-sensor data sets.

```{r echo=FALSE}
library(dplyr)
j_test <- read.csv("data/sample_corrupted_data.csv",
                   stringsAsFactors = FALSE)

j_test <- j_test %>%
  select(Station, sampleDateTime, Temperature_C) %>%
  filter(Station == "CB") 

j_test <- j_test[1:1000,] 
```


## Function demos

### `numeric_range_checker(dat, parms, ...)`

```{r echo=FALSE}

numeric_range_checker <- function(dat,
                                  parms,
                                  ...) {
  # Description: checks the range to determine within acceptible realistic ranges, 
  # specified by user (e.g., relative humidity must be between 0 and 100, range of 
  # callibrated temperatures). 
  
  # Inputs: dat a data.frame with 
  # Returns: 0 = pass; 1 = fail; -1 = cannot be used
  # Tests: range_test1, range_test2
  # Examples: 
  # parms <- data.frame()
  # numeric_range_checker(dat = range_test1, 
  #                       parms = parms)
  # Team members: Celeste
  
  if(!("range_checks" %in% names(parms))) 
    stop("range_checks must be a column in parms table")
  if( length(parms$range_checks) != 2 ) 
    stop("range_checks must be a length 2 with a lower and upper bound")
  
  if( parms$range_checks[1] >= parms$range_checks[2] ) 
    stop("lower bound of range_checks must be less than upper bound")
  
  if (!is.numeric(dat$values)) {
    attr(dat, "range_fails") <- rep(-1, length(dat$values))
  } else {
    lower_lim <- parms$range_checks[1]
    upper_lim <- parms$range_checks[2]
    
    attr(dat, "range_fails") <- which( ( dat$values < lower_lim) | (dat$values > upper_lim)  )
  }
  return(dat)
}

```

#### Test catch values that exceed abs threshold

```{r range-exceed-test}

# Test data
set.seed(42)
t <- 1:300
values <- 25 * sin(t) + rnorm(length(t), 0, 10)
dates <- seq.Date(as.Date("2017-01-01"), by = "day", length.out = length(t))
range_test1 <- data.frame(dates = dates, 
                          values = values)
print("test data")
print(head(range_test1))

# Parms magnitude between -40 and 40
parms <- data.frame(range_checks = c(-40, 40))
print(parms)

# Test it
test1 <- numeric_range_checker(dat = range_test1, 
                               parms = parms)

# Result (TRUE???)
print(all.equal(attributes(test1)$range_fails, 
          c(18, 36, 246, 269, 300)) )
```

#### Cast attribute flags

```{r}

test1$range_fails <- 0
test1$range_fails[attributes(test1)$range_fails] <- 1

print(paste0( round(100*sum(test1$range_fails)/nrow(test1), 2), "% of rows contained errors"))

```

### `numeric_spike_checker(x, spike_min = NA, verbose = FALSE)`


```{r echo=FALSE}
numeric_spike_checker <- function(x, spike_min = NA, verbose = FALSE) {
  # Description: Tests for spikes in data
  # Inputs: vector to be tested, optional threshold for how large a spike 
  # must be, verbose = TRUE/FALSE
  # Returns: result contains 0 for no spike, 1 for spike detected and -1 
  # for test could not be performed.  If verbose is selected,  the 
  # threshold value and a summary of the frequency of the codes is printed. 
  # # Team members: John, Patrick, Kaelin
  
  # if spike_min isn't specified, try to define a spike minimum based on 
  # the variance of the data stream. Default is 1.96 * standard deviation
  # of the sequential differences
  # NOTE if the default step is used, about 5% of the data will be flagged
  # even if there are no real problems with the data
  spike_min <- ifelse(is.na(spike_min), 
                      (1.96 * sd(diff(x, lag=1), na.rm=TRUE)), 
                      spike_min)
  # use 3 element moving window to detect spikes
  lag1 <- c(NA, head(x, -1)) # lag-1
  lag_minus1 <- c(tail(x, -1), NA) # lag+1
  s <- cbind(lag1, lag_minus1, x)
  #calculate the median
  median_vec <- apply(s[, c("lag1", "lag_minus1", "x")], 1, median)
  is_spike_val<- ((abs(x - median_vec) > spike_min) 
                  & (!is.na(lag1)) & (!is.na(lag_minus1)))
  result <- ifelse(is_spike_val, 1, 0)
  result <- ifelse(is.na(result), -1, result)
  if(verbose) {
    print(paste0("Looking for jumps or drops in the data > ",
                 spike_min))
    print("Frequency of result codes:")
    print(summary(as.factor(result)))
  }
  return(result)
}


```

#### Identify outliers

```{r}

out <- numeric_spike_checker(j_test$Temperature_C,
                             spike_min = 10, 
                             verbose = TRUE)

attr(j_test, "check_spikes") <- out

plot(as.Date(j_test$sampleDateTime), j_test$Temperature_C, type = "l",
     xlab = "Date",
     ylab = "Temperature C")
points(as.Date(j_test$sampleDateTime[which(out == 1)]),
       j_test$Temperature_C[which(out == 1)],
       col = "red", pch = 19)

```


### `time_gap_checker(dat, datetimeVar, gapThreshold)`

```{r}
time_gap_checker <- function(dat, datetimeVar, gapThreshold) {
  
  # confirm that the date time variable is of type posixct or posixlt
  if (!(is.POSIXct(tower[['timestamp']]) | is.POSIXlt(tower[['timestamp']]))) {
    stop("variable to check must be of type POSIXct or type POSIXlt")
  }
  
  datetimeVec <- dat[[datetimeVar]]
  timeGap <- datetimeVec - lag(datetimeVec)
  
  # set gap threshold to the standard deviation of the time step if not provided
  if (missing(gapThreshold)) {
    gapThreshold <- sd(timeGap, na.rm = T)
  }
  
  sequenceGap <- timeGap > gapThreshold
  
  # set flags according to:
  # -1 not run
  # 0 passed
  # 1 fail
  sequenceGap <- ifelse(is.na(sequenceGap), -1,
                        ifelse(sequenceGap == TRUE, 1,
                               ifelse(sequenceGap == FALSE, 0, NA)
                        )
  )
  
  return(sequenceGap)
  
}

```

### Gap checker test

```{r echo=FALSE}

library(lubridate)
library(dplyr)
library(readr)

tower <- read_csv('data/636_tower_data_dbg_26824efddd275586c075df1c01326ba0.csv') %>% 
  slice(1:500) %>% # pare down to a reasonable size
  slice(-c(10:50)) # introduce an intentional data gap



```

### demo with threshold
```{r time_gap_checker1 }

out <- time_gap_checker(dat = tower,
                 datetimeVar = 'timestamp',
                 gapThreshold = 30)

print(head(out))
```

### demo sans threshold

```{r time_gap_checker2 }

out_noT <- time_gap_checker(dat = tower,
                 datetimeVar = 'timestamp')

print(head(out_noT))

```


# Next steps

- See issues