############################################################################
# Title           	:   BrazilClim: script to gauge-calibrate the surfaces 
# Task          	:   Create gauge-calibrated montlhy climate surfaces
# Reference 	:   Ramoni-Perazzi, P. et al. BrazilClim: The overcoming of limitations 
#                                       of preexisting bioclimate data
# Written by 	:   P. Ramoni-Perazzi 
# Institution    	:   Federal University of Lavras
# Input         	:   Climatic information measured on-field (Supplementry Material A; DOI:
#                                       10.4121/14932947), surfaces from TRMM 3B43 v7 and NOAA, as 
#                                        well as the GMTED2010 digital elevation model  (please,
#                                        download them from the respective sources).
# Output      	:   Produces the monthly surfaces for each one of the covariates 
#                                       required to create the bioclimatic dataset: maximum and 
#                                       minimum temperatures, as well as precipitation. To select which
#                                       interpolation technique should be applied in each case, a previous 
#                                       selection is recommendable, based on some goodness-of-fit 
#                                      measurement such as RMSE (calculated through leave-one-out or 
#                                     any other cross-validation technique). These surfaces are created 
#                                     in GeoTiff format and 30 arc-sec resolution. 

rm(list=ls())
library(raster)
library(gstat)
library(doParallel)
library(fields)
library(spm) # Must be installed from source in newer versions of R
library(quantregForest)

# Activates the cluster 
cl <- parallel::makeCluster(detectCores())
doParallel::registerDoParallel(cl)

# Sets root path
pathx <- "C:/To/Documents/Dir/"

# Sets the methods to be tested and mode
metodos <- c("OK",  "UK",  "GBM", "TPS", "LM", "RF")
modo <- c("Dirct", "Error")

WGS84 <- "+proj=longlat +datum=WGS84" # WGS84 projection

# Locates the satellite images
trmm <- list.files(paste0(pathx, "GeoDataBases/TRMM_3B43/ppt"), full.names = T)[1:12]
noaaTmax <- list.files(paste0(pathx, "GeoDataBases/NOAA/tmax"), full.names = T)[1:12]
noaaTmin <- list.files(paste0(pathx, "GeoDataBases/NOAA/tmin"), full.names = T)[1:12]    
n <- c(1, 5:12, 2:4) # Fix the order of the layers
e <- extent(-73.99181, -34.6, -33.75014, 5.266527) # The extent of the study area

# Sets the output dir
op <- paste0(pathx, "To/Output/Directory/")

# Measuraments from gauges (BrazilClim_SupplMat_A)
dta <- read.csv(paste0("To/Dir/cointaining/SupplMat_A.csv"))

# In TPS:
# Z is a separate argument in Krig.predict, so we need a new function
# Internally (in interpolate) a matrix is formed of x, y, and elev (Z)
pfun <- function(model, x, ...) {
        predict(model, x[,1:2], Z=x[,3], ...)
}

#=========== PPT ===========
ppt <- dta[which(dta$Variable == "Ppt"), - c(1,2)]
gauges <- ppt[,c(1,2)]
coordinates(gauges) <- ~Lon + Lat
proj4string(gauges) <- WGS84

for(i in 3:ncol(ppt)){
   print(paste("PPT", i-2, "/", Sys.time()))
   trmmTemp <- raster(trmm[n[i-2]])
   trmmTemp <- crop(trmmTemp, e)
   proj4string(trmmTemp) <- WGS84
   valTemp <- extract(trmmTemp, gauges)
   br <- rasterToPoints(trmmTemp)
   br <- as.data.frame(br)
   colnames(br) <- c("LON", "LAT", "TRMM")
   coordinates(br) <- ~LON+LAT
   proj4string(br) <- WGS84
   names(trmmTemp) <- "TRMM"
   #-------------------------- CK --------------------------
   # Direct interpolations using TRMM as covariate --------------------------
   train <- ppt[,c(1,2,i)]
   colnames(train) <- c("LON", "LAT", "VAR")
   train$TRMM <- valTemp
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   mCK <- gstat(formula=VAR~TRMM, locations=train, nmax=floor(2.5*nrow(ppt)/100))  
   temp <- predict(mCK, newdata=br, debug.level=0)$var1.pred
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(trmmTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "CK_prec_Dirct_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "CK_prec_Dirct_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mCK, temp, temp2, r)
   # Interpolations using error (obs - TRMM) --------------------------
   train@data$VAR <- train@data$VAR-train@data$TRMM # Calculates error: Obs - TRMM
   mCK <- gstat(formula=VAR~TRMM, locations=train, nmax=floor(2.5*nrow(ppt)/100))  
   temp <- predict(mCK, newdata=br, debug.level=0)$var1.pred
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(trmmTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "CK_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "CK_prec_ERROR_", (i-2), ".tif")
   }
   r <-  r + trmmTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mCK, temp, temp2, r)

   #-------------------------- GBM --------------------------
   # Direct interpolations using TRMM as covariate --------------------------
   train <- ppt[,c(1,2,i)]
   colnames(train) <- c("LON", "LAT", "VAR")
   train$TRMM <- valTemp
   trainGBM <- as.data.frame(train)
   testGBM <- as.data.frame(br)
   temp <- gbmpred(trainGBM[, -3], trainGBM[, 3], testGBM[, c(1:2)],
                                         testGBM, family = "gaussian")[3]
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp$Predictions
   r <- setValues(trmmTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "GBM_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "GBM_prec_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(trainGBM, temp, temp2, r)
   # Interpolations using error (obs - TRMM) --------------------------
   train[,3] <- train[,3]-train[,4] # Calculates error: Obs - TRMM
   trainGBM <- as.data.frame(train)
   temp <- gbmpred(trainGBM[, -3], trainGBM[, 3], testGBM[, c(1:2)],
                                         testGBM, family = "gaussian")[3]
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp$Predictions
   r <- setValues(trmmTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "GBM_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "GBM_prec_ERROR_", (i-2), ".tif")
   }
   r <-  r + trmmTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(trainGBM, testGBM, temp, temp2, r)

   #-------------------------- LM --------------------------
   # Direct interpolations using TRMM as covariate --------------------------
   train <- ppt[,c(1,2,i)]
   colnames(train) <- c("LON", "LAT", "VAR")
   train$TRMM <- valTemp
   mLM <- lm(VAR~TRMM, train)
   temp <- predict(mLM, br)
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(trmmTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "LM_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "LM_prec_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mLM, temp, temp2, r)
   # Interpolations using error (obs - TRMM) --------------------------
   train[,3] <- train[,3]-train[,4] # Calculates error: Obs - TRMM
   mLM <- lm(VAR~TRMM, train)
   temp <- predict(mLM, br)
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(trmmTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "LM_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "LM_prec_ERROR_", (i-2), ".tif")
   }
   r <-  r + trmmTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mLM, temp, temp2, r)

   #-------------------------- RF --------------------------
   # Direct interpolations using TRMM as covariate --------------------------
   train <- ppt[,c(1,2,i)]
   colnames(train) <- c("LON", "LAT", "VAR")
   train$TRMM <- valTemp
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   mRF <- quantregForest(x=train@data[c("TRMM")], y=train$VAR)
   br1 <- br[c(1:3404811),]
   tempA <- predict(mRF, br1@data, y=br1$VAR)[,2]
   br2 <- br[c(3404812:6809623),]
   tempB <- predict(mRF, br2@data, y=br2$VAR)[,2]
   br3 <- br[c(6809624:nrow(br)),]
   tempC <- predict(mRF, br3@data, y=br3$VAR)[,2]
   temp <- c(tempA, tempB, tempC)
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(trmmTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "RF_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "RF_prec_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mRF, tempA, tempB, tempC, br1, br2, br3, temp, temp2, r)
   # Interpolations using error (obs - TRMM) --------------------------
   train@data$VAR <- train@data$VAR - train@data$TRMM
   mRF <- quantregForest(x=train@data[c("TRMM")], y=train$VAR)
   br1 <- br[c(1:3404811),]
   tempA <- predict(mRF, br1@data, y=br1$VAR)[,2]
   br2 <- br[c(3404812:6809623),]
   tempB <- predict(mRF, br2@data, y=br2$VAR)[,2]
   br3 <- br[c(6809624:nrow(br)),]
   tempC <- predict(mRF, br3@data, y=br3$VAR)[,2]
   temp <- c(tempA, tempB, tempC)
   temp2 <- getValues(trmmTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(trmmTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "RF_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "RF_prec_ERROR_", (i-2), ".tif")
   }
   r <-  r + trmmTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mRF, tempA, tempB, tempC, br1, br2, br3, temp, temp2, r)

   #-------------------------- TPS --------------------------
   # Direct interpolations using TRMM as covariate --------------------------
   train <- ppt[,c(1,2,i)]
   colnames(train) <- c("LON", "LAT", "VAR")
   train$TRMM <- valTemp
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   temp <- trmmTemp
   trainTPS <- list()
   trainTPS[[1]] <- train@coords
   trainTPS[[2]] <- train@data$VAR
   trainTPS[[3]] <- train@data$TRMM
   names(trainTPS) <- c("Coords", "VAR", "TRMM")
   mTPS <- Tps(trainTPS$Coords, trainTPS$VAR, Z = trainTPS$TRMM)
   r <- raster::interpolate(temp, mTPS, xyOnly=FALSE, fun=pfun)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "TPS_prec_Dirct_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "TPS_prec_Dirct_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(temp, trainTPS, mTPS, r)
   # Interpolations using error (obs - TRMM) --------------------------
   train@data$VAR <- train@data$VAR - train@data$TRMM
   temp <- trmmTemp
   trainTPS <- list()
   trainTPS[[1]] <- train@coords
   trainTPS[[2]] <- train@data$VAR
   trainTPS[[3]] <- train@data$TRMM
   names(trainTPS) <- c("Coords", "VAR", "TRMM")
   mTPS <- Tps(trainTPS$Coords, trainTPS$VAR, Z = trainTPS$TRMM)
   r <- raster::interpolate(temp, mTPS, xyOnly=FALSE, fun=pfun)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "TPS_prec_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "TPS_prec_ERROR_", (i-2), ".tif")
   }
   r <-  r + trmmTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(temp, trainTPS, mTPS, r)
}

################# TEMPERATURES #################
lats <- dem <- raster("To/DEM?Directory/dem.asc")
for(i in 1:nrow(lats)){
   lats[i,] <- yFromRow(dem, i)
}
lats <- crop(lats, e)
dem <- crop(dem, e)
proj4string(lats) <- proj4string(dem) <- WGS84

#=========== TMAX ===========
tmax <- dta[which(dta$Variable == "Tmax"), - c(1,2,17)]
gauges <- tmax[,c(1,2)]
coordinates(gauges) <- ~Lon + Lat
proj4string(gauges) <- WGS84
tmax$LAT1 <- extract(lats, gauges)
tmax$ELEV <- extract(dem, gauges)

for(i in 3:ncol(tmax)){
   print(paste("TMAX", i-2, "/", Sys.time()))
   noaaTemp <- raster(noaaTmax[n[i-2]])
   noaaTemp <- crop(noaaTemp, e)
   proj4string(noaaTemp) <- WGS84
   valTemp <- extract(noaaTemp, gauges)
   br <- rasterToPoints(noaaTemp)
   br <- as.data.frame(br)
   colnames(br) <- c("LON", "LAT", "NOAA")
   coordinates(br) <- ~LON+LAT
   proj4string(br) <- WGS84   
   br$LAT1 <- extract(lats, br)
   br$ELEV <- extract(dem, br)
   names(br) <- c("NOAA", "LAT1", "ELEV")
   #-------------------------- CK --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmax[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   mCK <- gstat(formula=VAR~NOAA+LAT1+ELEV, locations=train, nmax=floor(2.5*nrow(tmax)/100))  
   temp <- predict(mCK, newdata=br, debug.level=0)$var1.pred
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "CK_tmax_Dirct_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "CK_tmax_Dirct_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mCK, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train@data$VAR <- train@data$VAR-train@data$NOAA # Calculates error: Obs - NOAA
   mCK <- gstat(formula=VAR~NOAA+LAT1+ELEV, locations=train, nmax=floor(2.5*nrow(tmax)/100))  
   temp <- predict(mCK, newdata=br, debug.level=0)$var1.pred
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "CK_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "CK_tmax_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mCK, temp, temp2, r)

   #-------------------------- GBM --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmax[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   trainGBM <- as.data.frame(train)
   testGBM <- as.data.frame(br)
   temp <- gbmpred(trainGBM[, -3], trainGBM[, 3], testGBM[, c(1:2)],
                                         testGBM, family = "gaussian")[3]
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp$Predictions
   r <- setValues(noaaTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "GBM_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "GBM_tmax_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(trainGBM, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train[,3] <- train[,3]-train[,6] # Calculates error: Obs - NOAA
   trainGBM <- as.data.frame(train)
   temp <- gbmpred(trainGBM[, -3], trainGBM[, 3], testGBM[, c(1:2)],
                                         testGBM, family = "gaussian")[3]
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp$Predictions
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "GBM_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "GBM_tmax_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(trainGBM, testGBM, temp, temp2, r)

   #-------------------------- LM --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmax[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   mLM <- lm(VAR~NOAA+LAT1+ELEV, train)
   temp <- predict(mLM, br)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "LM_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "LM_tmax_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mLM, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train[,3] <- train[,3]-train[,6] # Calculates error: Obs - NOAA
   mLM <- lm(VAR~NOAA+LAT1+ELEV, train)
   temp <- predict(mLM, br)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "LM_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "LM_tmax_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mLM, temp, temp2, r)

   #-------------------------- RF --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmax[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   mRF <- quantregForest(x=train@data[c("NOAA", "LAT1", "ELEV")], y=train$VAR)
   br1 <- br[c(1:3404811),]  # Split in chunks to allow running it in less powerful computers
   tempA <- predict(mRF, br1@data, y=br1$VAR)[,2]
   br2 <- br[c(3404812:6809623),]
   tempB <- predict(mRF, br2@data, y=br2$VAR)[,2]
   br3 <- br[c(6809624:nrow(br)),]
   tempC <- predict(mRF, br3@data, y=br3$VAR)[,2]
   temp <- c(tempA, tempB, tempC)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "RF_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "RF_tmax_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mRF, tempA, tempB, tempC, br1, br2, br3, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train@data$VAR <- train@data$VAR - train@data$NOAA
   mRF <- quantregForest(x=train@data[c("NOAA", "LAT1", "ELEV")], y=train$VAR)
   br1 <- br[c(1:3404811),]
   tempA <- predict(mRF, br1@data, y=br1$VAR)[,2]
   br2 <- br[c(3404812:6809623),]
   tempB <- predict(mRF, br2@data, y=br2$VAR)[,2]
   br3 <- br[c(6809624:nrow(br)),]
   tempC <- predict(mRF, br3@data, y=br3$VAR)[,2]
   temp <- c(tempA, tempB, tempC)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "RF_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "RF_tmax_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mRF, tempA, tempB, tempC, br1, br2, br3, temp, temp2, r)

   #-------------------------- TPS --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmax[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   temp <- noaaTemp
   trainTPS <- list()
   trainTPS[[1]] <- train@coords
   trainTPS[[2]] <- train@data$VAR
   trainTPS[[3]] <- train@data$NOAA
   trainTPS[[4]] <- train@data$LAT
   trainTPS[[5]] <- train@data$ELEV
   names(trainTPS) <- c("Coords", "VAR", "NOAA", "LAT", "ELEV")
   mTPS <- Tps(trainTPS$Coords, trainTPS$VAR, Z = trainTPS$NOAA+train$LAT+ train$ELEV) 
   r <- raster::interpolate(temp, mTPS, xyOnly=FALSE, fun=pfun)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "TPS_tmax_Dirct_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "TPS_tmax_Dirct_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(temp, trainTPS, mTPS, r)
   # Interpolations using error (obs - NOAA) --------------------------
   trainTPS[[2]] <- train@data$VAR - train@data$NOAA
   mTPS <- Tps(trainTPS$Coords, trainTPS$VAR, Z = trainTPS$NOAA+train$LAT+ train$ELEV) 
   r <- raster::interpolate(temp, mTPS, xyOnly=FALSE, fun=pfun)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "TPS_tmax_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "TPS_tmax_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(temp, trainTPS, mTPS, r)
}

#=========== TMIN ===========
tmin <- dta[which(dta$Variable == "Tmin"), - c(1,2,17)]
gauges <- tmin[,c(1,2)]
coordinates(gauges) <- ~Lon + Lat
proj4string(gauges) <- WGS84
tmin$LAT1 <- extract(lats, gauges)
tmin$ELEV <- extract(dem, gauges)

for(i in 3:ncol(tmin)){
   print(paste("TMIN", i-2, "/", Sys.time()))
   noaaTemp <- raster(noaaTmin[n[i-2]])
   noaaTemp <- crop(noaaTemp, e)
   proj4string(noaaTemp) <- WGS84
   valTemp <- extract(noaaTemp, gauges)
   br <- rasterToPoints(noaaTemp)
   br <- as.data.frame(br)
   colnames(br) <- c("LON", "LAT", "NOAA")
   coordinates(br) <- ~LON+LAT
   proj4string(br) <- WGS84   
   br$LAT1 <- extract(lats, br)
   br$ELEV <- extract(dem, br)
   names(br) <- c("NOAA", "LAT1", "ELEV")
   #-------------------------- CK --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmin[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   mCK <- gstat(formula=VAR~NOAA+LAT1+ELEV, locations=train, nmax=floor(2.5*nrow(tmin)/100))  
   temp <- predict(mCK, newdata=br, debug.level=0)$var1.pred
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "CK_tmin_Dirct_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "CK_tmin_Dirct_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mCK, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train@data$VAR <- train@data$VAR-train@data$NOAA # Calculates error: Obs - NOAA
   mCK <- gstat(formula=VAR~NOAA+LAT1+ELEV, locations=train, nmax=floor(2.5*nrow(tmin)/100))  
   temp <- predict(mCK, newdata=br, debug.level=0)$var1.pred
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "CK_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "CK_tmin_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mCK, temp, temp2, r)

   #-------------------------- GBM --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmin[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   trainGBM <- as.data.frame(train)
   testGBM <- as.data.frame(br)
   temp <- gbmpred(trainGBM[, -3], trainGBM[, 3], testGBM[, c(1:2)],
                                         testGBM, family = "gaussian")[3]
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp$Predictions
   r <- setValues(noaaTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "GBM_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "GBM_tmin_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(trainGBM, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train[,3] <- train[,3]-train[,6] # Calculates error: Obs - NOAA
   trainGBM <- as.data.frame(train)
   temp <- gbmpred(trainGBM[, -3], trainGBM[, 3], testGBM[, c(1:2)],
                                         testGBM, family = "gaussian")[3]
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp$Predictions
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "GBM_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "GBM_tmin_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(trainGBM, testGBM, temp, temp2, r)

   #-------------------------- LM --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmin[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   mLM <- lm(VAR~NOAA+LAT1+ELEV, train)
   temp <- predict(mLM, br)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "LM_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "LM_tmin_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mLM, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train[,3] <- train[,3]-train[,6] # Calculates error: Obs - NOAA
   mLM <- lm(VAR~NOAA+LAT1+ELEV, train)
   temp <- predict(mLM, br)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   if(i < 12){
       nombre <- paste0(op, "LM_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "LM_tmin_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mLM, temp, temp2, r)

   #-------------------------- RF --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmin[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   mRF <- quantregForest(x=train@data[c("NOAA", "LAT1", "ELEV")], y=train$VAR)
   br1 <- br[c(1:3404811),]  # Split in chunks to allow running it in less powerful computers
   tempA <- predict(mRF, br1@data, y=br1$VAR)[,2]
   br2 <- br[c(3404812:6809623),]
   tempB <- predict(mRF, br2@data, y=br2$VAR)[,2]
   br3 <- br[c(6809624:nrow(br)),]
   tempC <- predict(mRF, br3@data, y=br3$VAR)[,2]
   temp <- c(tempA, tempB, tempC)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "RF_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "RF_tmin_ERROR_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mRF, tempA, tempB, tempC, br1, br2, br3, temp, temp2, r)
   # Interpolations using error (obs - NOAA) --------------------------
   train@data$VAR <- train@data$VAR - train@data$NOAA
   mRF <- quantregForest(x=train@data[c("NOAA", "LAT1", "ELEV")], y=train$VAR)
   br1 <- br[c(1:3404811),]
   tempA <- predict(mRF, br1@data, y=br1$VAR)[,2]
   br2 <- br[c(3404812:6809623),]
   tempB <- predict(mRF, br2@data, y=br2$VAR)[,2]
   br3 <- br[c(6809624:nrow(br)),]
   tempC <- predict(mRF, br3@data, y=br3$VAR)[,2]
   temp <- c(tempA, tempB, tempC)
   temp2 <- getValues(noaaTemp)
   temp2[which(!is.na(temp2))] <- temp
   r <- setValues(noaaTemp, temp2)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "RF_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "RF_tmin_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(mRF, tempA, tempB, tempC, br1, br2, br3, temp, temp2, r)

   #-------------------------- TPS --------------------------
   # Direct interpolations using NOAA as covariate --------------------------
   train <- tmin[,c(1,2,i, 15, 16)]
   colnames(train) <- c("LON", "LAT", "VAR", "LAT1", "ELEV")
   train$NOAA <- valTemp
   train <- train[complete.cases(train),]
   coordinates(train) <- ~LON+LAT
   proj4string(train) <- WGS84
   temp <- noaaTemp
   trainTPS <- list()
   trainTPS[[1]] <- train@coords
   trainTPS[[2]] <- train@data$VAR
   trainTPS[[3]] <- train@data$NOAA
   trainTPS[[4]] <- train@data$LAT
   trainTPS[[5]] <- train@data$ELEV
   names(trainTPS) <- c("Coords", "VAR", "NOAA", "LAT", "ELEV")
   mTPS <- Tps(trainTPS$Coords, trainTPS$VAR, Z = trainTPS$NOAA+train$LAT+ train$ELEV) 
   r <- raster::interpolate(temp, mTPS, xyOnly=FALSE, fun=pfun)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "TPS_tmin_Dirct_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "TPS_tmin_Dirct_", (i-2), ".tif")
   }
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(temp, trainTPS, mTPS, r)
   # Interpolations using error (obs - NOAA) --------------------------
   trainTPS[[2]] <- train@data$VAR - train@data$NOAA
   mTPS <- Tps(trainTPS$Coords, trainTPS$VAR, Z = trainTPS$NOAA+train$LAT+ train$ELEV) 
   r <- raster::interpolate(temp, mTPS, xyOnly=FALSE, fun=pfun)
   proj4string(r) <- WGS84
   if(i < 12){
       nombre <- paste0(op, "TPS_tmin_ERROR_0", (i-2), ".tif")
   } else {
       nombre <- paste0(op, "TPS_tmin_ERROR_", (i-2), ".tif")
   }
   r <-  r + noaaTemp
   writeRaster(r, nombre, format="GTiff", overwrite=TRUE)
   rm(temp, trainTPS, mTPS, r)
}
parallel::stopCluster(cl)