fieldAnomaly: Calculate the daily (julian or month/day date) or monthly...

In marchtaylor/sinkr: Collection of functions with emphasis in multivariate data analysis

Calculate the daily (julian or month/day date) or monthly anomaly of a data field

Description

The fieldAnomaly function calculates an anomaly field by subtracting the daily or monthly means from each spatial location (columns) of a field.

Usage

fieldAnomaly(y, x, level = "month")

Arguments

y	A matrix of the spatio-temporal field with rows being the temporal dimension and columns being the spatial dimension
x	A vector of the class "Date", "POSIXct", or "POSIXlt" containing the dates that correspond to the rows of matrix y.
level	Character string. Anomaly to compute ("julian", "day" or "month"). "julian" will consider unique as.POSIXlt(x)$yday values, while "day" will consider unique format(x,"%m%d") strings.

Examples

set.seed(1)
Time <- seq.Date(as.Date("1990-01-01"), as.Date("1996-12-31"), by="day")
Space <- seq(pi,2*pi,,10)
Signal1 <- sin(as.POSIXlt(Time)$yday/365*(2*pi)) # annual signal
Signal2 <- sin(as.POSIXlt(Time)$yday/365*(24*pi)) # monthly signal
cumSignal <- 1*Signal1+0.3*Signal2
plot(Time, cumSignal, t="l")
abline(v=seq(as.Date("1900-01-01"), as.Date("2100-01-01"), by="year"), lty=2, col=8)
Z <- outer(cumSignal, sin(Space))
Noise <- Z * 0 * array(rnorm(length(Z)), dim=dim(Z))
Z <- Z + Noise

# anomaly by month (removes "Signal1")
Z.anom <- fieldAnomaly(y=Z, x=Time, level="month")
zran <- c(-1,1) * max(abs(range(Z, Z.anom)))
pal <- colorPalette(c("blue4", "cyan", "grey90", "yellow", "red4"), c(20,1,1,20))
op <- par(no.readonly=TRUE)
layout(matrix(c(1,2,3,3), nrow=2, ncol=2), widths=c(4,1), heights=c(2,2))
par(mar=c(5,5,3,1))
image(Time, Space, Z, col=pal(101), zlim=zran, main="Original")
image(Time, Space, Z.anom, col=pal(101), zlim=zran, main="Anomaly")
par(mar=c(5,0,3,5))
imageScale(Z, col=pal(101), zlim=zran, axis.pos=4)
mtext("Value", side=4, line=3)
par(op)


# anomaly by day of year (i.e. date, removes "Signal1" and "Signal2")
Z.anom <- fieldAnomaly(y=Z, x=Time, level="day")
zran <- c(-1,1) * max(abs(range(Z, Z.anom)))
pal <- colorPalette(c("blue4", "cyan", "grey90", "yellow", "red4"), c(20,1,1,20))
op <- par(no.readonly=TRUE)
layout(matrix(c(1,2,3,3), nrow=2, ncol=2), widths=c(4,1), heights=c(2,2))
par(mar=c(5,5,3,1))
image(Time, Space, Z, col=pal(101), zlim=zran, main="Original")
image(Time, Space, Z.anom, col=pal(101), zlim=zran, main="Anomaly")
par(mar=c(5,0,3,5))
imageScale(Z, col=pal(101), zlim=zran, axis.pos=4)
mtext("Value", side=4, line=3)
par(op)


# anomaly by julian day (day of year, removes "Signal1" and "Signal2")
Z.anom <- fieldAnomaly(y=Z, x=Time, level="julian")
zran <- c(-1,1) * max(abs(range(Z, Z.anom)))
pal <- colorPalette(c("blue4", "cyan", "grey90", "yellow", "red4"), c(20,1,1,20))
op <- par(no.readonly=TRUE)
layout(matrix(c(1,2,3,3), nrow=2, ncol=2), widths=c(4,1), heights=c(2,2))
par(mar=c(5,5,3,1))
image(Time, Space, Z, col=pal(101), zlim=zran, main="Original")
image(Time, Space, Z.anom, col=pal(101), zlim=zran, main="Anomaly")
par(mar=c(5,0,3,5))
imageScale(Z, col=pal(101), zlim=zran, axis.pos=4)
mtext("Value", side=4, line=3)
par(op)


# test that anomalies are removed correctly
Zalt <- as.matrix(apply(scale(Z),1,mean))
plot(Time, Zalt, t="l")

tmp <- aggregate(Zalt ~ as.POSIXlt(Time)$yday, FUN=mean)
plot(strptime(paste("2000",tmp[,1], sep="-"), format="%Y-%j"), tmp[,2],
     t="o", xlab="", ylab="n values", main="level = 'julian'")

tmp <- aggregate(Zalt ~ format(Time, "%m-%d"), FUN=mean)
plot(strptime(paste("2000",tmp[,1], sep="-"), format="%Y-%m-%d"), tmp[,2], 
     t="o", xlab="", ylab="n values", main="level = 'day'")


Z.anom <- fieldAnomaly(y=Zalt, x=Time, level="month")
sqrt(mean(Z.anom^2))
plot(Time, Z.anom, t="l")

Z.anom <- fieldAnomaly(y=Zalt, x=Time, level="day")
sqrt(mean(Z.anom^2))
plot(Time, Z.anom, t="l")

Z.anom <- fieldAnomaly(y=Zalt, x=Time, level="julian")
sqrt(mean(Z.anom^2))
plot(Time, Z.anom, t="l")

marchtaylor/sinkr documentation built on July 4, 2022, 5:48 p.m.