Cacher is the code snippet organizer for pro developers

#thanhlv 
#May, 2015
#Inspired by the paper: http://www.parasitesandvectors.com/content/7/1/302

library(raster)
library(TSA)
library(parallel)

# Path to the folder having 12 monthly LST images
folder <- "/home/thanhlv/LST/"
#Scan tiff files inside the folder
tiffs <- list.files(folder, pattern = "*\\.tif$", full.names = TRUE)
# Create a stack raster of 12 monthly LST
lst <- stack(tiffs)

# Create an harmonic function for calc 
fun_harmonic <- function(x){
  # Convert values of a pixel over 12 months into a vector
  vt <- as.vector(x, mode = "numeric")
  # And then convert them to a timeseries structure
  vt.ts <- ts(vt, start = c(2008, 1), frequency = 12)
  # Perform harmonic regression using harmonic function of the TSA package
  har <- harmonic(vt.ts, 3)
  # Extract coef from the model of har coefs and timeseries data
  coef(lm(x ~ har))
}
#Note: This is assumed that the lst have no NA values. In case any pixels have NA, it is suggested to replace NA with a certain value (i.e -9999). Otherwise, lm function will return errors 0-NA values

#Single run of harmonic regression
s_harmonic_lst <- calc(lst_monthly, fun = fun_harmonic)
#The harmonic_lst has 7 layers which are equivalent to 7 harmonic coefficients
#
#Parallel computing to speed up the calculation

#Parallel computing
#Get available cpus for prallel. Only use no of cpu - 1, leave 1 for other tasks
cpus <- detectCores() - 1

beginCluster(n = cpus) #Initialize a cluster of cpus
# Create a custom function calc for parallel computing
f1 <- function(x) calc(x, fun_harmonic)
# Running parallel calc function
P_harmonic_lst <- clusterR(lst, f1, export = 'fun_harmonic')
# After finished, stop cluster, release occupied CPUs
endCluster()