# Have to load the knitr to use hooks library(knitr) ## This vignettes name vignettename <- "setup-data" # REMEMBER: IF CHANGING IN THE shared-init (next block), then copy to the others!
# Width will scale all figwidth <- 12 # Scale the wide figures (100% out.width) figheight <- 4 # Heights for stacked time series plots figheight1 <- 5 figheight2 <- 6.5 figheight3 <- 8 figheight4 <- 9.5 figheight5 <- 11 # Set the size of squared figures (same height as full: figheight/figwidth) owsval <- 0.35 ows <- paste0(owsval*100,"%") ows2 <- paste0(2*owsval*100,"%") # fhs <- figwidth * owsval # Set for square fig: fig.width=fhs, fig.height=fhs, out.width=ows} # If two squared the: fig.width=2*fhs, fig.height=fhs, out.width=ows2 # Check this: https://bookdown.org/yihui/rmarkdown-cookbook/chunk-styling.html # Set the knitr options knitr::opts_chunk$set( collapse = TRUE, comment = "##!! ", prompt = FALSE, cache = TRUE, cache.path = paste0("../tmp/vignettes/tmp-",vignettename,"/"), fig.align="center", fig.path = paste0("../tmp/vignettes/tmp-",vignettename,"/"), fig.height = figheight, fig.width = figwidth, out.width = "100%" ) options(digits=3) # For cropping output and messages cropfun <- function(x, options, func){ lines <- options$output.lines ## if (is.null(lines)) { ## return(func(x, options)) # pass to default hook ## } if(!is.null(lines)){ x <- unlist(strsplit(x, "\n")) i <- grep("##!!",x) if(length(i) > lines){ # truncate the output, but add .... #x <- c(x[-i[(lines+1):length(i)]], "```", "## ...output cropped", "```") x <- x[-i[(lines+1):length(i)]] x[i[lines]] <- pst(x[i[lines]], "\n\n## ...output cropped") } # paste these lines together x <- paste(c(x, ""), collapse = "\n") } x <- gsub("!!","",x) func(x, options) } hook_chunk <- knit_hooks$get("chunk") knit_hooks$set(chunk = function(x, options) { cropfun(x, options, hook_chunk) })
This vignette explains how to setup data consisting of observations and forecasts, such that it can be used for onlineforecast models. A generic introduction and description is in available in onlineforecasting. The code is available here. More information on onlineforecasting.org.
First load the package:
# Load the package library(onlineforecast)
In the package different data sets are included. The
Dbuilding
holds the data used for the example of
heat load forecasting in the building-heat-load-forecasting vignette.
When the package is loaded the data is also loaded, so we can access it directly. Let's start out by:
# Keep it in D to simplify notation D <- Dbuilding
The class is 'data.ĺist':
# The class of D class(D)
Actually, a 'data.list' is simply a 'list', but we made the class 'data.list' in order to have functions for the particular format of data - the format is explained in this document.
It consists of vectors of time, vectors of observations (model output) and data.frames of forecasts (model input):
# Print the names to see the variables in the data names(D)
An overview of the content can be generated by:
summary.default(D)
where it can be seen that t
is a time vector, heatload
is a vector, and Ta
and I
are data.frames.
A function giving a summary, including checks of the format of the 'data.list' is:
summary(D)
The 'NA' columns indicate the proportion of NAs. If there is a ok
in a column,
then the check of the variables format is passed. See the help with
?summary.data.list
to learn which checks are performed.
First, lets have a look at D$t
, which is the vector of time points:
# The time class(D$t) head(D$t) tail(D$t)
Hence, the vector is of the class POSIXct
. It is not a necessity, t
can also simply be a numeric, but for plotting and many operations, its
very useful to use the 'POSIXct' class (see ?POSIXt
).
Rules for the time vector:
It must be named t
.
There must be no gaps or NA values in t
, since only equidistant time series
can be used in the models (the other variables can have NAs).
Its best to keep the time zone in UTC
or GMT
(not providing any time
zone tz
can give rise to problems).
Use the basic R functions for handling the time class. Most needed operations can be done with:
?as.POSIXct ?strftime
A helper function is provided with the ct
function which can be called using ?
, or ?ct
. See example below:
# Convert from a time stamp (tz="GMT" per default) ct("2019-01-01 11:00") # Convert from unix time ct(3840928387)
Note that for all functions where a time value as a character is given, the time
zone is always "GMT" (or "UTC", but this can result in warnings, but they can be
ignored). For some operations the package lubridate
can
be very helpful.
Note the rules for observations:
In a data.list
observations must be vectors.
The vectors must have the same length as the time t
vector.
Observation as numerical vectors can be used directly as model output (if observations are to used as model inputs, they must be setup in a data.frame as explained below in Section [Forecasts]).
In the current data, a time series of hourly heat load observations is included:
str(D$heatload)
It must have the same length as the time vector:
# Same length as time length(D$t) length(D$heatload)
A simple plot can be generated by:
plot(D$t, D$heatload, type="l", xlab="Time", ylab="Headload (kW)")
The convention used in all examples is that the time points are always set to the time interval end point, e.g.:
# The observation D$heatload[2] # Represents the average load between D$t[1] # and D$t[2]
The main idea behind setting the time point at the end of the interval is: Working with values averaged over the time interval, such values are available at the end of the time interval, not before. Especially, in real-time applications this is a useful convention.
As described in onlineforecasting the setup of forecasts for model inputs always follows the same format - as presented in the following. This is also the format of the forecasts generated by functions in the package. Hence all forecasts must follow this format.
The rules are:
All values at row i
are available at the i
'th value in time t
.
All columns must be named with k
followed by an integer indicating the
horizon in steps (e.g. the column named k8
hold the 8-step forecasts).
Have a look at the forecasts of the global radiation:
# Global radiation forecasts head(D$I)
At the first time point:
# First time point D$t[1]
the available forecast ahead in time is at the first row:
# The forecast available ahead in time is in the first row D$I[1, ]
We can plot that by:
i <- 1:ncol(D$I) plot(D$t[i], D$I[1, ], type="l", xlab="Time", ylab="Global radiation forecast (I in W/m²)")
So this is the forecast available ahead in time at r D$t[1]
.
The column in I
named k8
holds the 8-step horizon forecasts, which, since
the steps are hourly, is an equi-distant time series. Picking out the
entire series can be done by D$I$k8
- hence a plot (together with the
observations) can be generated by:
# Just pick some points by i <- 200:296 plot(D$t[i], D$I$k8[i], type="l", col=2, xlab="Time", ylab="Global radiation (W/m²)") # Add the observations lines(D$t[i], D$Iobs[i]) legend("topright", c("8-step forecasts","Observations"), bg="white", lty=1, col=2:1)
Notice how the are not aligned, since the forecasts are 8 hours ahead. To align them the forecasts must be lagged 8 steps by:
plot(D$t[i], lagvec(D$I$k8[i], 8), type="l", col=2, xlab="Time", ylab="Global radiation (W/m²)") lines(D$t[i], D$Iobs[i]) legend("topright", c("8-step forecasts lagged","Observations"), bg="white", lty=1, col=2:1)
A few simple plotting functions are included in the package.
The plot function provided with the package actually does this lagging with plotting forecasts:
plot_ts(D, patterns=c("^I"), c("2010-12-15","2010-12-18"), kseq=c(1,8,24,36))
The argument patterns
is vector of a regular expressions (see ?regex
),
which is used to match the variables to include in the plot. See the help with ?plot_ts
for more details.
An interactive plot can be generated using (first install the package plotly
):
plotly_ts(D, patterns=c("heatload$","^I"), c("2010-12-15","2010-12-18"), kseq=c(1,8,24,36))
L <- plotly_ts(D, patterns=c("heatload$","^I"), c("2010-12-15","2010-12-18"), kseq=c(1,8,24,36), plotit=FALSE) plotly::subplot(L, shareX=TRUE, nrows=length(L), titleY = TRUE)
Note that the patterns
argument is a vector of regular expressions, which
determines which variables from D
to plot.
When modelling with the objective of forecasting, it's always a good start to have a look at scatter plots between the model inputs and the model output. For example the heatload vs. ambient temperature 8-step forecast:
par(mfrow=c(1,2)) plot(D$Ta$k8, D$heatload) plot(lagvec(D$Ta$k8, 8), D$heatload)
So lagging (thus aligning in time) makes less slightly less scatter.
A wrapper for the pairs
function is provided for a data.list
, which can
generate very useful explorative plots:
pairs(D, nms=c("heatload","Taobs","Ta","t"), kseq=c(1,8,24))
Note how the sequence of included horizons are specified in the kseq
argument,
and note that the forecasts are lagged to be aligned in time. See ?pairs.data.list
for more details.
Just as a quick side note: This is the principle used for fitting onlineforecast models, simply shift forecasts to align with the observations:
# Lag the 8-step forecasts to be aligned with the observations x <- lagvec(D$I$k8, 8) # Take a smaller range x <- x[i] # Take the observations y <- D$Iobs[i] # Fit a linear regression model fit <- lm(y ~ x) # Plot the result plot(x, y, xlab="8-step forecasts (W/m²)", ylab="Obsservations (W/m²)", main="Global radiation") abline(fit)
Seen over time the 8-step forecasts are:
plot(D$t[i], predict.lm(fit, newdata=data.frame(x)), type="l", ylim=c(0,max(y)), xlab="Time", ylab="Global radiation (W/m^2)", col=2) lines(D$t[i], y) legend("topright", c("8-step forecasts lagged","Observations"), lty=1, col=2:1)
Of course that model was very simple, see how to make a better model in [building-heat-load-forecasting] and more information on the [website].
Taking a subset of a data.list
is very useful and it can easily be done in
different ways using the subset
function (i.e. it's really the
subset.data.list
function called when:
# Take the 1 to 4 values of each variable in D Dsub <- subset(D, 1:4) summary(Dsub)
Another useful function for taking data in a time range is:
which(in_range("2010-12-20",D$t,"2010-12-21"))
always check the help of function for more details (i.e. ?in_range
)
Actually, it's easy to take subset from a period by:
Dsub <- subset(D, c("2010-12-20","2010-12-21")) summary(Dsub) Dsub$t
It can be really useful to bring the data.list on a format of a data.frame
or
equivalently data.table
for processing.
Bringing to data.frame
can easily be done by:
Df <- as.data.frame(Dsub) names(Df)
So the forecasts are just bind with the time and observations, and .kxx
is
added to the column names.
It can be converted to a data.table
by:
library(data.table) setDT(Df) class(Df)
After processing it is easily converted back to the data.list
again by:
# Set back to data.frame setDF(Df) # Convert to a data.list Dsub2 <- as.data.list(Df) # Compare it with the original Dsub summary(Dsub2) summary(Dsub)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.