This package includes functions to read in and analyze accelerometer data generated by Oregon Research Electronic AL100 Acceleration Loggers to estimate tree phenology.
#Install from github devtools::install_github("agougher/accel")
#Load packages require(accel) require(forecast) require(pracma) require(gtools)
To get a sense of what raw accelerometer data looks like, there are several files available here. If you download these files to your working directory, you can analyze them with the code below. The season-length example data below was processed in the same way, but for about 250 files.
#List files in order and calculates dominant period x <- mixedsort(list.files("./", pattern="CSV", full.names=TRUE)) ex1 <- lapply(x,readAndDom) ex1 <- data.frame(do.call(rbind, ex1)) #Plotting plot(1:5, ex1$peakY,pch=20, col="red") points(1:5, ex1$peakZ, pch=20)
The example data was generated by analyzing 250 acceleration files collected between March 1, 2016 and November 5, 2016, using the code below. The raw data was not included here because of its size, so the code is for illustrative purposes only. readAndDom
reads in accelerometer files and calculates the dominant period of the three axes, as well as a weight. Reading in and calculating the dominant periods will take a while, so it is best to run in parallel when possible.
#NOT RUN #x is a list of ordered csv file locations al3 <- lapply(x,readAndDom) al3 <- do.call(rbind, al3) #Assign day of year al3["day"] <- 61:310
The example data includes the dominant period values for each of three dimensional axes (x, y, z), and weights reflective of the difference between the dominant and second-most dominant periods for each axis, for a single accelerometer.
#Load example data and make long form dataframe for functions below data(al3) head(al3) #Extract y and z axes, and bind them for the functions below domY <- al3[,c("day","peakY","wpeakY")] domZ <- al3[,c("day","peakZ","wpeakZ")] colnames(domY) <- c("day","value","weights") colnames(domZ) <- c("day","value","weights") dat <- rbind(domY,domZ) #Plot data, point size is proportional to the weight plot(dat$day, dat$value, cex=dat$weights*2, xlab="Day of year",ylab="Dominant period")
This function, rmOut
, identifies and removes points that have residuals beyond 1.5x the interquartile range of a loess curve. In the plot below, points in purple are kept after removing outliers.
#The 'rmOut' function requires a dataframe with three columns named "day","value", and "weights" dat2 <- rmOut(dat) plot(dat$day, dat$value, xlab="Day of year",ylab="Dominant period") points(dat2$day, dat2$value, pch=20, col="purple")
The model is a dual logistic curve, with an additional parameter that controls the slope of the line between the two logistic curves. In the plot below, black points are cleaned dominant period values, and red points are fitted values from the model. The two vertical lines are the spring and autumn inflection points of the model.
#Remove NA's from data frame dat2 <- na.omit(dat2) #Set starting parameters m <- c(20, 19, 145, 6, 275, 3, 1) #Fit non-linear least squares model mod = nls(value ~ a + (b- g*day)*((1/(1+exp((c-day)/d))) - (1/(1+exp((e-day)/f)))), start = list(a = m[1],b = m[2],c = m[3], d = m[4], e = m[5], f = m[6], g = m[7] ), algorithm="port",control = list(maxiter = 500), data=dat2, weights=dat2$weights^2) #Plot data, and fitted values #Units are in tenths of seconds, so divide by 10 to convert to seconds plot(dat2$day, dat2$value/10, pch=20, xlab="Day of year",ylab="Dominant period (s)") points(dat2$day, fitted(mod)/10, pch=20, col="red", lwd=4) abline(v=c(coef(mod)[3],coef(mod)[5]))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.