Hello,

I'm trying to fit a sine curve over successive temperature readings (i.e. minimum and maximum temperature) over several days and for many locations. The code below shows a hypothetical example of 5000 locations with  7 days of temperature data. Not very efficient when you have many more locations and days.

The linear interpolation takes 0.7 seconds, and the sine interpolations take 2 to 4 seconds depending on the approach.

Any ideas on how to speed this up? Thanks in advance.

Ariel

### R Code ######

# 1- Prepare data fake data
  days<- 7
  n <- 5000*days
  tmin <- matrix(rnorm(n, mean=0) , ncol=days, nrow=5000)
  tmax <- matrix(rnorm(n, mean=10), ncol=days, nrow=5000)
  m <- matrix(NA, ncol=days*2, nrow=5000)
  m[,seq(1,ncol(m),2)]  <- tmin
  m[,seq(2,ncol(m)+1,2)]<- tmax
  # check first row
  plot(1:ncol(m), m[1,], type="l")

# 2 -linear interpolation: 0.66 seconds
  xout <- seq(0,ncol(m),0.25/24*2)[-1] # time step = 0.25 hours or 15 minutes
  system.time( m1 <- t(apply(m,1, function(y) approx(x=1:ncol(m), y=y, xout=xout, method="linear")$y)) )
  # Check first row
  plot(1:ncol(m), m[1,], type="l")
  points(xout, m1[1,], col="red", cex=1)


# 3- sine interpolation
  sine.approx1 <- function(index, tmin, tmax) {
    b <- (2*pi)/24  # period = 24 hours
    c <- pi/2 # horizontal shift
    xout <- seq(0,24,0.25)[-1]
    yhat <- apply(cbind(tmin[index,],tmax[index,]), 1, function(z) diff(z)/2 * sin(b*xout-c) + mean(z))
    #yhat <- yhat[-nrow(yhat),]
    yhat <- c(yhat)
    #plot(yhat, type="l")
  }
  sine.approx2 <- function(index, tmin, tmax) {
    b <- (2*pi)/24  # period = 24 hours
    c <- pi/2 # horizontal shift
    xout1 <- seq(0 ,12,0.25)
    xout2 <- seq(12,24,0.25)[-1]
    xout2 <- xout2[-length(xout2)]
    yhat1 <- apply(cbind(tmin[index,]                       ,tmax[index,]    ), 1, function(z) diff(z)/2 * sin(b*xout1-c) + mean(z))
    yhat2 <- apply(cbind(tmax[index,][-length(tmax[index,])],tmin[index,][-1]), 1, function(z) diff(z)/2 * sin(b*xout2+c) + mean(z))
    yhat2 <- cbind(yhat2,NA)
    yhat3 <- rbind(yhat1,yhat2)
    #yhat3 <- yhat3[-nrow(yhat3),]
    yhat3 <- c(yhat3)
    yhat <- yhat3
    #plot(c(yhat1))
    #plot(c(yhat2))
    #plot(yhat, type="l")
  }

  # Single sine: 2.23 seconds
  system.time( m2 <- t(sapply(1:nrow(m), function(i) sine.approx1(i, tmin=tmin, tmax=tmax))) )

  # Double sine: 4.03 seconds
  system.time( m3 <- t(sapply(1:nrow(m), function(i) sine.approx2(i, tmin=tmin, tmax=tmax))) )

  # take a look at approach 1
  plot(seq(-1,ncol(m)-1,1)[-1], m[1,], type="l")
  points(xout, m2[1,], col="red", cex=1)

  # take a look at approach 2
  plot(seq(-1,ncol(m)-1,1)[-1], m[1,], type="l")
  points(xout, m3[1,], col="blue", cex=1)


---
Ariel Ortiz-Bobea
Fellow
Resources for the Future
1616 P Street, N.W.
Washington, DC 20036
202-328-5173


	[[alternative HTML version deleted]]

