library("fda.usc")
library("pls")
library("caret")


library(XLConnect)

# STATIONS OF CÁDIZ
wb <- loadWorkbook('C:/Datos3/Estaciones_CÁDIZ.xlsx', create = TRUE) 
Villamartin <- readWorksheet(wb, sheet = 'Villamartin')

# STATIONS OF CÓRDOBA
wb <- loadWorkbook('C:/Datos3/Estaciones_CÓRDOBA.xlsx', create = TRUE)
Adamuz <- readWorksheet(wb, sheet = 'Adamuz')
Baena <- readWorksheet(wb, sheet = 'Baena')
Belmez <- readWorksheet(wb, sheet = 'Belmez')
Cabra <- readWorksheet(wb, sheet = 'Cabra')
Cordoba <- readWorksheet(wb, sheet = 'Córdoba')
El_Carpio <- readWorksheet(wb, sheet = 'El Carpio')
Hinojosa_Duque <- readWorksheet(wb, sheet = 'Hinojosa del Duque')
Hornachuelos <- readWorksheet(wb, sheet = 'Hornachuelos')
Palma_Rio <- readWorksheet(wb, sheet = 'Palma del Río')
Santaella<- readWorksheet(wb, sheet = 'Santaella')

# STATIONS OF GRANADA
wb <- loadWorkbook('C:/Datos3/Estaciones_GRANADA.xlsx', create = TRUE)
Loja <- readWorksheet(wb, sheet = 'Loja')
Pinos_Puente <- readWorksheet(wb, sheet = 'Pinos Puente')

# STATIONS OF JAÉN
wb <- loadWorkbook('C:/Datos3/Estaciones_JAÉN.xlsx', create = TRUE)
Alcaudete <- readWorksheet(wb, sheet = 'Alcaudete')
Chiclana_Segura <- readWorksheet(wb, sheet = 'Chiclana de Segura')
Jaen <- readWorksheet(wb, sheet = 'Jaén')
Higuera_Arjona <- readWorksheet(wb, sheet = 'La Higuera de Arjona')
Mancha_Real <- readWorksheet(wb, sheet = 'Mancha Real')
Marmolejo <- readWorksheet(wb, sheet = 'Marmolejo')
Pozo_Alcon <- readWorksheet(wb, sheet = 'Pozo Alcón')
San_Jose_Propios <- readWorksheet(wb, sheet = 'San José de los Propios')
Santo_Tome <- readWorksheet(wb, sheet = 'Santo Tomé')

# STATIONS OF MÁLAGA
wb <- loadWorkbook('C:/Datos3/Estaciones_MÁLAGA.xlsx', create = TRUE)
Antequera <- readWorksheet(wb, sheet = 'Antequera')
Archidona <- readWorksheet(wb, sheet = 'Archidona')
Pizarra <- readWorksheet(wb, sheet = 'Pizarra')
Sierra_Yeguas <- readWorksheet(wb, sheet = 'Sierra Yeguas')

# STATIONS OF SEVILLA
wb <- loadWorkbook('C:/Datos3/Estaciones_SEVILLA.xlsx', create = TRUE)
Ecija <- readWorksheet(wb, sheet = 'Écija')
Osuna <- readWorksheet(wb, sheet = 'Osuna')

#############################################################
# lista_DATOS is a list of dimensions 28 x 2191 x 7, where 28 is the number of stations,
# 2191 is the number of days of the 6 years, between 2005 and 2010, and there are 7 
# agro-climatic measurements:
#############################################################

lista_DATOS <- list(Villamartín=Villamartin,
       Adamuz=Adamuz,Baena=Baena,Belmez=Belmez,Cabra=Cabra,Córdoba=Cordoba,El_Carpio=El_Carpio,
       Hinojosa_Duque=Hinojosa_Duque,Hornachuelos=Hornachuelos,Palma_Río=Palma_Rio,Santaella=Santaella,
       Loja=Loja,Pinos_Puente=Pinos_Puente,
       Alcaudete=Alcaudete,Chiclana_Segura=Chiclana_Segura,Jaén=Jaen,Higuera_Arjona=Higuera_Arjona,
       Mancha_Real=Mancha_Real,Marmolejo=Marmolejo,Pozo_Alcón=Pozo_Alcon,San_José_Propios=San_Jose_Propios,
       Santo_Tomé=Santo_Tome,Antequera=Antequera,Archidona=Archidona,Pizarra=Pizarra,
       Sierra_Yeguas=Sierra_Yeguas,Écija=Ecija,Osuna=Osuna)

#############################################################
# lista_MET is a list with 28 matrices (for each station).
# Each matrix has 72=12x6 rows (sum of measurements for each month))
# and 7 columns (one for each agro-climatic measurement)
#############################################################

matriz_MET <- matrix(0,nrow=72,ncol=7)

lista_MET <- list(Villamartín=matriz_MET,
       Adamuz=matriz_MET,Baena=matriz_MET,Belmez=matriz_MET,Cabra=matriz_MET,Córdoba=matriz_MET,El_Carpio=matriz_MET,
       Hinojosa_Duque=matriz_MET,Hornachuelos=matriz_MET,Palma_Río=matriz_MET,Santaella=matriz_MET,
       Loja=matriz_MET,Pinos_Puente=matriz_MET,
       Alcaudete=matriz_MET,Chiclana_Segura=matriz_MET,Jaén=matriz_MET,Higuera_Arjona=matriz_MET,
       Mancha_Real=matriz_MET,Marmolejo=matriz_MET,Pozo_Alcón=matriz_MET,San_José_Propios=matriz_MET,
       Santo_Tomé=matriz_MET,Antequera=matriz_MET,Archidona=matriz_MET,Pizarra=matriz_MET,
       Sierra_Yeguas=matriz_MET,Écija=matriz_MET,Osuna=matriz_MET)

for (i in 1:28)
{
  for (j in 1:7)
  {
  k <- j+3
  lista_MET[[i]][1,j] <- sum(as.numeric(lista_DATOS[[i]][1:31,k]))
  lista_MET[[i]][2,j] <- sum(as.numeric(lista_DATOS[[i]][32:59,k]))
  lista_MET[[i]][3,j] <- sum(as.numeric(lista_DATOS[[i]][60:90,k]))
  lista_MET[[i]][4,j] <- sum(as.numeric(lista_DATOS[[i]][91:120,k]))
  lista_MET[[i]][5,j] <- sum(as.numeric(lista_DATOS[[i]][121:151,k]))
  lista_MET[[i]][6,j] <- sum(as.numeric(lista_DATOS[[i]][152:181,k]))
  lista_MET[[i]][7,j] <- sum(as.numeric(lista_DATOS[[i]][182:212,k]))
  lista_MET[[i]][8,j] <- sum(as.numeric(lista_DATOS[[i]][213:243,k]))
  lista_MET[[i]][9,j] <- sum(as.numeric(lista_DATOS[[i]][244:273,k]))
  lista_MET[[i]][10,j] <- sum(as.numeric(lista_DATOS[[i]][274:304,k]))
  lista_MET[[i]][11,j] <- sum(as.numeric(lista_DATOS[[i]][305:334,k]))
  lista_MET[[i]][12,j] <- sum(as.numeric(lista_DATOS[[i]][335:365,k]))
  lista_MET[[i]][13,j] <- sum(as.numeric(lista_DATOS[[i]][366:396,k]))
  lista_MET[[i]][14,j] <- sum(as.numeric(lista_DATOS[[i]][397:424,k]))
  lista_MET[[i]][15,j] <- sum(as.numeric(lista_DATOS[[i]][425:455,k]))
  lista_MET[[i]][16,j] <- sum(as.numeric(lista_DATOS[[i]][456:485,k]))
  lista_MET[[i]][17,j] <- sum(as.numeric(lista_DATOS[[i]][486:516,k]))
  lista_MET[[i]][18,j] <- sum(as.numeric(lista_DATOS[[i]][517:546,k]))
  lista_MET[[i]][19,j] <- sum(as.numeric(lista_DATOS[[i]][547:577,k]))
  lista_MET[[i]][20,j] <- sum(as.numeric(lista_DATOS[[i]][578:608,k]))
  lista_MET[[i]][21,j] <- sum(as.numeric(lista_DATOS[[i]][609:638,k]))
  lista_MET[[i]][22,j] <- sum(as.numeric(lista_DATOS[[i]][639:669,k]))
  lista_MET[[i]][23,j] <- sum(as.numeric(lista_DATOS[[i]][670:699,k]))
  lista_MET[[i]][24,j] <- sum(as.numeric(lista_DATOS[[i]][700:730,k]))
  lista_MET[[i]][25,j] <- sum(as.numeric(lista_DATOS[[i]][731:761,k]))
  lista_MET[[i]][26,j] <- sum(as.numeric(lista_DATOS[[i]][762:789,k]))
  lista_MET[[i]][27,j] <- sum(as.numeric(lista_DATOS[[i]][790:820,k]))
  lista_MET[[i]][28,j] <- sum(as.numeric(lista_DATOS[[i]][821:850,k]))
  lista_MET[[i]][29,j] <- sum(as.numeric(lista_DATOS[[i]][851:881,k]))
  lista_MET[[i]][30,j] <- sum(as.numeric(lista_DATOS[[i]][882:911,k]))
  lista_MET[[i]][31,j] <- sum(as.numeric(lista_DATOS[[i]][912:942,k]))
  lista_MET[[i]][32,j] <- sum(as.numeric(lista_DATOS[[i]][943:973,k]))
  lista_MET[[i]][33,j] <- sum(as.numeric(lista_DATOS[[i]][974:1003,k]))
  lista_MET[[i]][34,j] <- sum(as.numeric(lista_DATOS[[i]][1004:1034,k]))
  lista_MET[[i]][35,j] <- sum(as.numeric(lista_DATOS[[i]][1035:1064,k]))
  lista_MET[[i]][36,j] <- sum(as.numeric(lista_DATOS[[i]][1065:1095,k]))
  lista_MET[[i]][37,j] <- sum(as.numeric(lista_DATOS[[i]][1096:1126,k]))
  lista_MET[[i]][38,j] <- sum(as.numeric(lista_DATOS[[i]][1127:1155,k]))
  lista_MET[[i]][39,j] <- sum(as.numeric(lista_DATOS[[i]][1156:1185,k]))
  lista_MET[[i]][40,j] <- sum(as.numeric(lista_DATOS[[i]][1186:1216,k]))
  lista_MET[[i]][41,j] <- sum(as.numeric(lista_DATOS[[i]][1217:1247,k]))
  lista_MET[[i]][42,j] <- sum(as.numeric(lista_DATOS[[i]][1248:1277,k]))
  lista_MET[[i]][43,j] <- sum(as.numeric(lista_DATOS[[i]][1278:1308,k]))
  lista_MET[[i]][44,j] <- sum(as.numeric(lista_DATOS[[i]][1309:1339,k]))
  lista_MET[[i]][45,j] <- sum(as.numeric(lista_DATOS[[i]][1340:1369,k]))
  lista_MET[[i]][46,j] <- sum(as.numeric(lista_DATOS[[i]][1370:1400,k]))
  lista_MET[[i]][47,j] <- sum(as.numeric(lista_DATOS[[i]][1401:1430,k]))
  lista_MET[[i]][48,j] <- sum(as.numeric(lista_DATOS[[i]][1431:1461,k]))
  lista_MET[[i]][49,j] <- sum(as.numeric(lista_DATOS[[i]][1462:1492,k]))
  lista_MET[[i]][50,j] <- sum(as.numeric(lista_DATOS[[i]][1493:1520,k]))
  lista_MET[[i]][51,j] <- sum(as.numeric(lista_DATOS[[i]][1521:1551,k]))
  lista_MET[[i]][52,j] <- sum(as.numeric(lista_DATOS[[i]][1552:1581,k]))
  lista_MET[[i]][53,j] <- sum(as.numeric(lista_DATOS[[i]][1582:1612,k]))
  lista_MET[[i]][54,j] <- sum(as.numeric(lista_DATOS[[i]][1613:1642,k]))
  lista_MET[[i]][55,j] <- sum(as.numeric(lista_DATOS[[i]][1643:1673,k]))
  lista_MET[[i]][56,j] <- sum(as.numeric(lista_DATOS[[i]][1674:1704,k]))
  lista_MET[[i]][57,j] <- sum(as.numeric(lista_DATOS[[i]][1705:1734,k]))
  lista_MET[[i]][58,j] <- sum(as.numeric(lista_DATOS[[i]][1735:1765,k]))
  lista_MET[[i]][59,j] <- sum(as.numeric(lista_DATOS[[i]][1766:1795,k]))
  lista_MET[[i]][60,j] <- sum(as.numeric(lista_DATOS[[i]][1796:1826,k]))
  lista_MET[[i]][61,j] <- sum(as.numeric(lista_DATOS[[i]][1827:1857,k]))
  lista_MET[[i]][62,j] <- sum(as.numeric(lista_DATOS[[i]][1858:1885,k]))
  lista_MET[[i]][63,j] <- sum(as.numeric(lista_DATOS[[i]][1886:1916,k]))
  lista_MET[[i]][64,j] <- sum(as.numeric(lista_DATOS[[i]][1917:1946,k]))
  lista_MET[[i]][65,j] <- sum(as.numeric(lista_DATOS[[i]][1947:1977,k]))
  lista_MET[[i]][66,j] <- sum(as.numeric(lista_DATOS[[i]][1978:2007,k]))
  lista_MET[[i]][67,j] <- sum(as.numeric(lista_DATOS[[i]][2008:2038,k]))
  lista_MET[[i]][68,j] <- sum(as.numeric(lista_DATOS[[i]][2039:2069,k]))
  lista_MET[[i]][69,j] <- sum(as.numeric(lista_DATOS[[i]][2070:2099,k]))
  lista_MET[[i]][70,j] <- sum(as.numeric(lista_DATOS[[i]][2100:2130,k]))
  lista_MET[[i]][71,j] <- sum(as.numeric(lista_DATOS[[i]][2131:2160,k]))
  lista_MET[[i]][72,j] <- sum(as.numeric(lista_DATOS[[i]][2161:2191,k]))
  }
}

##########################################
# READING OF FILES INFO AND NIR
##########################################

wb <- loadWorkbook('C:/Datos3/INFO.xlsx', create = TRUE)
INFO <- readWorksheet(wb, sheet = 'Hoja1')

NIR <- read.table("C:/Datos3/NIR.txt", header=T)
attach(NIR)
NIR <- as.matrix(NIR)

###################################################################################
# FUNCTION To assign the sum of METEOROLOGICAL measures:
###################################################################################

## x, station; y, harvest; y1, month1; y2, month2; z, agro-climatic measurement 1-7 
medida_MET_SUM <- function(x,y,y1,y2,z)
{
  # STATION:
  if (x==10) {i <- 1}
  if (x==20) {i <- 2}
  if (x==21) {i <- 3}
  if (x==22) {i <- 4}
  if (x==23) {i <- 5}
  if (x==24) {i <- 6}
  if (x==25) {i <- 7}
  if (x==26) {i <- 8}
  if (x==27) {i <- 9}
  if (x==28) {i <- 10}
  if (x==29) {i <- 11}
  if (x==30) {i <- 12}
  if (x==31) {i <- 13}       
  if (x==40) {i <- 14}
  if (x==41) {i <- 15}
  if (x==42) {i <- 16}
  if (x==43) {i <- 17}
  if (x==44) {i <- 18}
  if (x==45) {i <- 19}
  if (x==46) {i <- 20}
  if (x==47) {i <- 21}
  if (x==48) {i <- 22}
  if (x==50) {i <- 23}
  if (x==51) {i <- 24}
  if (x==52) {i <- 25}
  if (x==53) {i <- 26}
  if (x==60) {i <- 27}
  if (x==61) {i <- 28}
  # Average of stations:
  if (x==1053) {i <- 29}
  if (x==2027) {i <- 30}
  if (x==2129) {i <- 31}
  if (x==2329) {i <- 32}
  if (x==2427) {i <- 33}
  if (x==286061) {i <- 34}
  if (x==3031) {i <- 35}
  if (x==4148) {i <- 36}
  if (x==4243) {i <- 37}
  if (x==4345) {i <- 38}
  if (x==4748) {i <- 39}
  if (x==5051) {i <- 40}
  if (x==5053) {i <- 41}
  
  #HARVEST, MONTH:
  month1 <- y1
  if (y==2) {j1<-mes1}
  if (y==3) {j1<-mes1+12}
  if (y==4) {j1<-mes1+24}
  if (y==5) {j1<-mes1+36}
  if (y==6) {j1<-mes1+48}
  if (y==7) {j1<-mes1+60}
  dif <- y2-y1
  j2 <- j1+dif

  #AGRO-CLIMATIC MEASUREMENT:
  k <- z
  
  if (i<29)
  {
     s <- sum(as.numeric(lista_MET[[i]][j1:j2,k]))
  }
  if (i==29) {s <- mean(c(sum(as.numeric(lista_MET[[1]][j1:j2,k])),sum(as.numeric(lista_MET[[26]][j1:j2,k]))))}
  if (i==30) {s <- mean(c(sum(as.numeric(lista_MET[[2]][j1:j2,k])),sum(as.numeric(lista_MET[[9]][j1:j2,k]))))}
  if (i==31) {s <- mean(c(sum(as.numeric(lista_MET[[3]][j1:j2,k])),sum(as.numeric(lista_MET[[11]][j1:j2,k]))))}
  if (i==32) {s <- mean(c(sum(as.numeric(lista_MET[[5]][j1:j2,k])),sum(as.numeric(lista_MET[[11]][j1:j2,k]))))}
  if (i==33) {s <- mean(c(sum(as.numeric(lista_MET[[6]][j1:j2,k])),sum(as.numeric(lista_MET[[9]][j1:j2,k]))))}
  if (i==34) {s <- mean(c(sum(as.numeric(lista_MET[[10]][j1:j2,k])),sum(as.numeric(lista_MET[[27]][j1:j2,k])),sum(as.numeric(lista_MET[[28]][j1:j2,k]))))}
  if (i==35) {s <- mean(c(sum(as.numeric(lista_MET[[12]][j1:j2,k])),sum(as.numeric(lista_MET[[13]][j1:j2,k]))))}
  if (i==36) {s <- mean(c(sum(as.numeric(lista_MET[[15]][j1:j2,k])),sum(as.numeric(lista_MET[[22]][j1:j2,k]))))}
  if (i==37) {s <- mean(c(sum(as.numeric(lista_MET[[16]][j1:j2,k])),sum(as.numeric(lista_MET[[17]][j1:j2,k]))))}
  if (i==38) {s <- mean(c(sum(as.numeric(lista_MET[[17]][j1:j2,k])),sum(as.numeric(lista_MET[[19]][j1:j2,k]))))}
  if (i==39) {s <- mean(c(sum(as.numeric(lista_MET[[21]][j1:j2,k])),sum(as.numeric(lista_MET[[22]][j1:j2,k]))))}
  if (i==40) {s <- mean(c(sum(as.numeric(lista_MET[[23]][j1:j2,k])),sum(as.numeric(lista_MET[[24]][j1:j2,k]))))}
  if (i==41) {s <- mean(c(sum(as.numeric(lista_MET[[23]][j1:j2,k])),sum(as.numeric(lista_MET[[26]][j1:j2,k]))))}
  s
}

############################################################################
#To refill LIST with information of agro-climatic measurements for every month: 
############################################################################

m <- matrix(1, nrow=222, ncol=7)
lista_medida_MET <- list(m,m,m,m,m,m,m,m,m,m,m,m)

for (k in 1:12) #mes
  for (i in 1:222)
   for (j in 1:7)
   {
   lista_medida_MET[[k]][i,j] <- medida_MET_SUM(INFO$ESTACIÓN[i],INFO$AÑO[i],k,k,j) 
   }

y <- data.frame(CG)
#NIR
absorp.fdataN <- fdata(NIR,argvals = NULL ,rangeval=NULL,names=NULL,fdata2d=FALSE)
FDATOS_N <- list(absorp.fdataN = absorp.fdataN, y = y)
absorpN <- FDATOS_N$absorp.fdataN
absorpN.d1 <- fdata.deriv(absorpN, nderiv = 1)
absorpN.d2 <- fdata.deriv(absorpN, nderiv = 2)

cam <- ifelse(INFO$AÑO==2,1,ifelse(INFO$AÑO==3,2,ifelse(INFO$AÑO==4,3,ifelse(INFO$AÑO==5,4,ifelse(INFO$AÑO==6,5,6)))))

#########################################################################
# Figure 1
#########################################################################

labelsH <- c("H1","H2","H3","H4","H5","H6")
plot(absorpN[-ind,], col=cam, main = NULL, xlab="Wavelength (mm)", ylab="Absorbance",cex.main=2, cex.lab=1.5, cex.axis=1.3)
legend("topleft",col=1:6,pch=15,labelsH,cex=1.4)

#########################################################################
# Figure 2
#########################################################################

labelsH <- c("H1","H2","H3","H4","H5","H6")
par(mfcol=c(1,2));  
par(mar=c(5,4,4,2));  
par(oma=c(3,3,3,3));  
par(mar=c(4,4,2,2))  
plot(absorpN.d1[,], col=cam, main = NULL, xlab="Wavelength (mm)", ylab="d(Absorbance,1)")
legend("bottomleft",col=1:6,pch=15, labelsH)
plot(absorpN.d2[,], col=cam, main = NULL, xlab="Wavelength (mm)", ylab="d(Absorbance,2)")
legend("bottomleft",col=1:6,pch=15, labelsH)

#########################################################################
# Figure 3
#########################################################################

par(mfcol=c(2,2));  
par(mar=c(5,4,4,2));  
par(oma=c(3,3,3,3));  
par(mar=c(4,4,2,2))  
labelsH <- c("H1","H2","H3","H4","H5","H6")
plot(absorpN.d1[,40:60], legend=TRUE, col=cam, main = NULL, xlab="Wavelength (mm)", ylab="d(Absorbance,1)",cex=0.8,pt.cex=1)
legend("bottomleft",col=1:6,pch=15,labelsH)
plot(absorpN.d2[,650:660], col=cam, main = NULL, xlab="Wavelength (mm)", ylab="d(Absorbance,2)",cex=0.8,pt.cex=1)
legend("topright",col=1:6,pch=15, labelsH)
plot(absorpN.d2[,420:450], col=cam, main = NULL, xlab="Wavelength (mm)", ylab="d(Absorbance,2)",cex=0.8)
legend("bottomleft",col=1:6,pch=15, labelsH)
plot(absorpN.d1[,1090:1135], col=cam, main = NULL, xlab="Wavelength (mm)", ylab="d(Absorbance,1)",cex=0.8)
legend("bottomleft",col=1:6,pch=15,labelsH)

#########################################################################
# Figure 4
#########################################################################


media_lista_MET <- matrix(1,nrow=72,ncol=7)
for (i in 1:7)
{
   suma <- matrix(0,nrow=72,ncol=1)   
   for (j in 1:28)   
   {
     suma <- suma + lista_MET[[j]][,i]     
   }
   media_lista_MET[,i] <- suma/28
}

# To standardize 
media_lista_MET_tip <- matrix(1,nrow=72,ncol=7)
for (i in 1:72)
{
   for (j in 1:7)
   {
   media_lista_MET_tip[i,j] <- (media_lista_MET[i,j]-mean(media_lista_MET[,j]))/sd(media_lista_MET[,j])
   }
}

par(mfcol=c(3,1));  
par(mar=c(5,4,4,2));  
par(oma=c(3,3,3,3));  
par(mar=c(4,4,2,2)) 
ind <- 1
matplot(media_lista_MET_tip[,ind],axes=F,type="l",lwd=2,col=c(2),ylim=c(-2,4.5),xlab="")
legend("topleft",col=c(2),pch=18, c("Std Temp"),cex=1.5)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(-3:5))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(-3:5))
x1<- matrix(13,nrow=1,ncol=9) #C1 
x2<- matrix(25,nrow=1,ncol=9) #C2 
x3<- matrix(37,nrow=1,ncol=9) #C3
x4<- matrix(49,nrow=1,ncol=9) #C4 
x5<- matrix(61,nrow=1,ncol=9) #C5
y<-c(-3:5)
lines(x1,y,type="l",pch=22,col="black")
lines(x2,y,type="l",pch=22,col="black")
lines(x3,y,type="l",pch=22,col="black")
lines(x4,y,type="l",pch=22,col="black")
lines(x5,y,type="l",pch=22,col="black")
mtext("H1",side=1,line=1,at=c(7),cex=1.1,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.1,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.1,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.1,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.1,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.1,col=1)
ind <- 2
matplot(media_lista_MET_tip[,ind],axes=F,type="l",lwd=2,col=c(3),ylim=c(-2,4.5),xlab="")
legend("topleft",col=c(3),pch=18, c("Std Hum"),cex=1.5)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(-3:5))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(-3:5))
x1<- matrix(13,nrow=1,ncol=9) #C1 
x2<- matrix(25,nrow=1,ncol=9) #C2 
x3<- matrix(37,nrow=1,ncol=9) #C3
x4<- matrix(49,nrow=1,ncol=9) #C4 
x5<- matrix(61,nrow=1,ncol=9) #C5
y<-c(-3:5)
lines(x1,y,type="l",pch=22,col="black")
lines(x2,y,type="l",pch=22,col="black")
lines(x3,y,type="l",pch=22,col="black")
lines(x4,y,type="l",pch=22,col="black")
lines(x5,y,type="l",pch=22,col="black")
mtext("H1",side=1,line=1,at=c(7),cex=1.1,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.1,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.1,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.1,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.1,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.1,col=1)
ind <- 3
matplot(media_lista_MET_tip[,ind],axes=F,type="l",lwd=2,col=c(4),ylim=c(-2,4.5),xlab="")
legend("topleft",col=c(4),pch=18, c("Std WSpe"),cex=1.5)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(-3:5))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(-3:5))
x1<- matrix(13,nrow=1,ncol=9) #C1 
x2<- matrix(25,nrow=1,ncol=9) #C2 
x3<- matrix(37,nrow=1,ncol=9) #C3
x4<- matrix(49,nrow=1,ncol=9) #C4 
x5<- matrix(61,nrow=1,ncol=9) #C5
y<-c(-3:5)
lines(x1,y,type="l",pch=22,col="black")
lines(x2,y,type="l",pch=22,col="black")
lines(x3,y,type="l",pch=22,col="black")
lines(x4,y,type="l",pch=22,col="black")
lines(x5,y,type="l",pch=22,col="black")
mtext("H1",side=1,line=1,at=c(7),cex=1.1,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.1,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.1,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.1,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.1,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.1,col=1)


##################################################################
par(mfcol=c(3,1));  
par(mar=c(5,4,4,2));  
par(oma=c(3,3,3,3));  
par(mar=c(4,4,2,2)) 
ind <- 4
matplot(media_lista_MET_tip[,ind],axes=F,type="l",lwd=2,col=c(6),ylim=c(-2,4.5),xlab="")
legend("topleft",col=c(6),pch=18, c("Std WDir"),cex=1.5)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(-3:5))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(-3:5))
x1<- matrix(13,nrow=1,ncol=9) #C1 
x2<- matrix(25,nrow=1,ncol=9) #C2 
x3<- matrix(37,nrow=1,ncol=9) #C3
x4<- matrix(49,nrow=1,ncol=9) #C4 
x5<- matrix(61,nrow=1,ncol=9) #C5
y<-c(-3:5)
lines(x1,y,type="l",pch=22,col="black")
lines(x2,y,type="l",pch=22,col="black")
lines(x3,y,type="l",pch=22,col="black")
lines(x4,y,type="l",pch=22,col="black")
lines(x5,y,type="l",pch=22,col="black")
mtext("H1",side=1,line=1,at=c(7),cex=1.1,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.1,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.1,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.1,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.1,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.1,col=1)

ind <- 5
matplot(media_lista_MET_tip[,ind],axes=F,type="l",lwd=2,col="orange",ylim=c(-2,4.5),xlab="")
legend("topleft",col="orange",pch=18, c("Std Rad"),cex=1.5)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(-3:5))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(-3:5))
x1<- matrix(13,nrow=1,ncol=9) #C1 
x2<- matrix(25,nrow=1,ncol=9) #C2 
x3<- matrix(37,nrow=1,ncol=9) #C3
x4<- matrix(49,nrow=1,ncol=9) #C4 
x5<- matrix(61,nrow=1,ncol=9) #C5
y<-c(-3:5)
lines(x1,y,type="l",pch=22,col="black")
lines(x2,y,type="l",pch=22,col="black")
lines(x3,y,type="l",pch=22,col="black")
lines(x4,y,type="l",pch=22,col="black")
lines(x5,y,type="l",pch=22,col="black")
mtext("H1",side=1,line=1,at=c(7),cex=1.1,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.1,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.1,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.1,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.1,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.1,col=1)
ind <- 6
matplot(media_lista_MET_tip[,ind],axes=F,type="l",lwd=2,col="blue",ylim=c(-2,4.5),xlab="")
legend("topleft",col="blue",pch=18, c("Std Precip"),cex=1.5)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(-3:5))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(-3:5))
x1<- matrix(13,nrow=1,ncol=9) #C1 
x2<- matrix(25,nrow=1,ncol=9) #C2 
x3<- matrix(37,nrow=1,ncol=9) #C3
x4<- matrix(49,nrow=1,ncol=9) #C4 
x5<- matrix(61,nrow=1,ncol=9) #C5
y<-c(-3:5)
lines(x1,y,type="l",pch=22,col="black")
lines(x2,y,type="l",pch=22,col="black")
lines(x3,y,type="l",pch=22,col="black")
lines(x4,y,type="l",pch=22,col="black")
lines(x5,y,type="l",pch=22,col="black")
mtext("H1",side=1,line=1,at=c(7),cex=1.1,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.1,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.1,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.1,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.1,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.1,col=1)


##################################################################
par(mfcol=c(3,1));  
par(mar=c(5,4,4,2));  
par(oma=c(3,3,3,3));  
par(mar=c(4,4,2,2)) 

ind <- 7
matplot(media_lista_MET_tip[,ind],axes=F,type="l",lwd=2,col=c(8),ylim=c(-2,4.5),xlab="")
legend("topleft",col=c(8),pch=18, c("Std ETo"),cex=1.5)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(-3:5))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(-3:5))
x1<- matrix(13,nrow=1,ncol=9) #C1 
x2<- matrix(25,nrow=1,ncol=9) #C2 
x3<- matrix(37,nrow=1,ncol=9) #C3
x4<- matrix(49,nrow=1,ncol=9) #C4 
x5<- matrix(61,nrow=1,ncol=9) #C5
y<-c(-3:5)
lines(x1,y,type="l",pch=22,col="black")
lines(x2,y,type="l",pch=22,col="black")
lines(x3,y,type="l",pch=22,col="black")
lines(x4,y,type="l",pch=22,col="black")
lines(x5,y,type="l",pch=22,col="black")
mtext("H1",side=1,line=1,at=c(7),cex=1.1,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.1,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.1,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.1,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.1,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.1,col=1)


#########################################################################
# Figure 5
#########################################################################


# To accumulate 
media_lista_MET_ac <- matrix(0,nrow=72,ncol=7)
sum <- matrix(0,nrow=1,ncol=7)
for (i in 1:72)
{
   for (j in 1:7)
   {
   media_lista_MET_ac[i,j] <- sum[j] + media_lista_MET[i,j]
   }
   sum <- media_lista_MET_ac[i,]
}

# To accumulate for years

media_lista_MET_ac <- matrix(0,nrow=72,ncol=7)
sum <- matrix(0,nrow=1,ncol=7)
for (i in 1:12)
{
   for (j in 1:7)
   {
   media_lista_MET_ac[i,j] <- sum[j] + media_lista_MET[i,j]
   }
   sum <- media_lista_MET_ac[i,]
}
sum <- matrix(0,nrow=1,ncol=7)
for (i in 13:24)
{
   for (j in 1:7)
   {
   media_lista_MET_ac[i,j] <- sum[j] + media_lista_MET[i,j]
   }
   sum <- media_lista_MET_ac[i,]
}
sum <- matrix(0,nrow=1,ncol=7)
for (i in 25:36)
{
   for (j in 1:7)
   {
   media_lista_MET_ac[i,j] <- sum[j] + media_lista_MET[i,j]
   }
   sum <- media_lista_MET_ac[i,]
}
sum <- matrix(0,nrow=1,ncol=7)
for (i in 37:48)
{
   for (j in 1:7)
   {
   media_lista_MET_ac[i,j] <- sum[j] + media_lista_MET[i,j]
   }
   sum <- media_lista_MET_ac[i,]
}
sum <- matrix(0,nrow=1,ncol=7)
for (i in 49:60)
{
   for (j in 1:7)
   {
   media_lista_MET_ac[i,j] <- sum[j] + media_lista_MET[i,j]
   }
   sum <- media_lista_MET_ac[i,]
}
sum <- matrix(0,nrow=1,ncol=7)
for (i in 61:72)
{
   for (j in 1:7)
   {
   media_lista_MET_ac[i,j] <- sum[j] + media_lista_MET[i,j]
   }
   sum <- media_lista_MET_ac[i,]
}


matplot(media_lista_MET_ac[,6],type="b",pch=18,axes=F,ylab="Accum. Precip.",col="blue",lwd=3,cex=1.1)
axis(1,c(1,13,25,37,49,61,72))
axis(2,c(0,100,200,300,400,500,600,700,800,900))
axis(3,c(1,13,25,37,49,61,72))
axis(4,c(0,100,200,300,400,500,600,700,800,900))
x1<- matrix(13,nrow=1,ncol=901) #C1 
x2<- matrix(25,nrow=1,ncol=901) #C2 
x3<- matrix(37,nrow=1,ncol=901) #C3
x4<- matrix(49,nrow=1,ncol=901) #C4 
x5<- matrix(61,nrow=1,ncol=901) #C5
y<-c(0:900)
lines(x1,y,type="l",pch=18,col="darkgrey",lwd=2)
lines(x2,y,type="l",pch=18,col="darkgrey",lwd=2)
lines(x3,y,type="l",pch=18,col="darkgrey",lwd=2)
lines(x4,y,type="l",pch=18,col="darkgrey",lwd=2)
lines(x5,y,type="l",pch=18,col="darkgrey",lwd=2)
mtext("H1",side=1,line=1,at=c(7),cex=1.3,col=1)
mtext("H2",side=1,line=1,at=c(19),cex=1.3,col=1)
mtext("H3",side=1,line=1,at=c(31),cex=1.3,col=1)
mtext("H4",side=1,line=1,at=c(43),cex=1.3,col=1)
mtext("H5",side=1,line=1,at=c(55),cex=1.3,col=1)
mtext("H6",side=1,line=1,at=c(67),cex=1.3,col=1)

#########################################################################
# Figures 7 and 8
#########################################################################


#########################################################################################################################################################################################################
### ONLY NIR INFORMATION
#########################################################################################################################################################################################################

dtt <- read.table("C:/Data/INFO_CG_NIR_MET.txt", header=T)
attach(dtt)
dataset <- dtt
dataset[1,1:10]

###################################################################
### EVALUATE SOME ALGORITHMS
###################################################################

CG <- as.matrix(dtt[,8:10])
NIR <- as.matrix(dtt[,11:1247])

##################################### PCA
CG_NIR <- data.frame(CG,NIR)
COMP_PCA<- pcr(CG ~ NIR, ncomp=20, data=CG_NIR)
as.matrix(COMP_PCA$scores)

dataset2 <- cbind(C_PROVINCIA=as.factor(dataset$C_PROVINCIA),as.matrix(COMP_PCA$scores))
dataset2 <- as.data.frame(dataset2)
names <- c('C_PROVINCIA','PCA1','PCA2','PCA3','PCA4')
colnames(dataset2) <- names

# Run algorithms using 10-fold cross validation
control <- trainControl(method="cv", number=10)
metric <- "Accuracy"

#Let’s evaluate 5 different algorithms:
#Linear Discriminant Analysis (LDA)
#Classification and Regression Trees (CART).
#k-Nearest Neighbors (kNN).
#Support Vector Machines (SVM) with a linear kernel.
#Random Forest (RF)

ind_CV <- c(1:25)
for (j in ind_CV)
{

###Create training/validation datasets
set.seed(j)
validation_index <- createDataPartition(dataset2$C_PROVINCIA, p=0.80, list=FALSE)
validation <- dataset2[-validation_index,]
dataset2 <- dataset2[validation_index,]
dim(dataset2)
dim(validation)

# a) linear algorithms
set.seed(j)
fit.lda <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="lda", metric=metric, trControl=control)
# b) nonlinear algorithms
# CART
set.seed(j)
fit.cart <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="rpart", metric=metric, trControl=control)
# kNN
set.seed(j)
fit.knn <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="knn", metric=metric, trControl=control)
# c) advanced algorithms
# SVM
set.seed(j)
fit.svm <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="svmRadial", metric=metric, trControl=control)
# Random Forest
set.seed(j)
fit.rf <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="rf", metric=metric, trControl=control)
# summarize accuracy of models
results <- resamples(list(lda=fit.lda, cart=fit.cart, knn=fit.knn, svm=fit.svm, rf=fit.rf))
summary(results) results
results$values
ls(results$values)
mean(results$values[,2])
mean(results$values[,9])
dotplot(results,xlim=c(0,1)) #RF provides the best results

CM_Accuracy_PROVINCE <- matrix(0,nrow=5,ncol=13)
CM_Kappa_PROVINCE <- matrix(0,nrow=5,ncol=13)

# MAKE PREDICTIONS
# estimate skill of LDA on the validation dataset
predictions <- predict(fit.lda, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[1,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[1,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.cart, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[2,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[2,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.knn, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[3,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[3,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.svm, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[4,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[4,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.rf, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[5,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[5,13] <- CM$overall[[2]] # Kappa

#########################################################################################################################################################################################################
### NIR + MET INFORMATION
#########################################################################################################################################################################################################

mat_Accuracy_PROVINCE <- matrix(0,nrow=5,ncol=13)
mat_Kappa_PROVINCE <- matrix(0,nrow=5,ncol=13)

mat_Accuracy_PROVINCE[1,13] <- mean(results$values[,2])
mat_Accuracy_PROVINCE[2,13] <- mean(results$values[,4])
mat_Accuracy_PROVINCE[3,13] <- mean(results$values[,6])
mat_Accuracy_PROVINCE[4,13] <- mean(results$values[,8])
mat_Accuracy_PROVINCE[5,13] <- mean(results$values[,10])

mat_Kappa_PROVINCE[1,13] <- mean(results$values[,3])
mat_Kappa_PROVINCE[2,13] <- mean(results$values[,5])
mat_Kappa_PROVINCE[3,13] <- mean(results$values[,7])
mat_Kappa_PROVINCE[4,13] <- mean(results$values[,9])
mat_Kappa_PROVINCE[5,13] <- mean(results$values[,11])

for (i in 1:12)
{
dataset3 <- cbind(C_PROVINCIA=as.factor(dataset$C_PROVINCIA),as.matrix(COMP_PCA$scores),as.matrix(lista_medida_MET[[i]]))
dataset3 <- as.data.frame(dataset3)
names <- c('C_PROVINCIA','PCA1','PCA2','PCA3','PCA4',
            'Temp','Hum','WSped','WDic','Rad','Precip','Eto')
colnames(dataset3) <- names

# Run algorithms using 10-fold cross validation
control <- trainControl(method="cv", number=10)
metric <- "Accuracy"

#Let’s evaluate 5 different algorithms:
#Linear Discriminant Analysis (LDA)
#Classification and Regression Trees (CART).
#k-Nearest Neighbors (kNN).
#Support Vector Machines (SVM) with a linear kernel.
#Random Forest (RF)

ind_CV <- c(1:25)
for (j in ind_CV)
{

###Create training/validation datasets
set.seed(j)
validation_index <- createDataPartition(dataset3$C_PROVINCIA, p=0.80, list=FALSE)
validation <- dataset3[-validation_index,]
dataset3 <- dataset3[validation_index,]
dim(dataset3)

# a) linear algorithms
set.seed(j)
fit.lda <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="lda", metric=metric, trControl=control)
# b) nonlinear algorithms
# CART
set.seed(j)
fit.cart <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="rpart", metric=metric, trControl=control)
# kNN
set.seed(j)
fit.knn <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="knn", metric=metric, trControl=control)
# c) advanced algorithms
# SVM
set.seed(j)
fit.svm <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="svmRadial", metric=metric, trControl=control)
# Random Forest
set.seed(j)
fit.rf <- train(as.factor(C_PROVINCIA)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="rf", metric=metric, trControl=control)
# summarize accuracy of models
results2 <- resamples(list(lda=fit.lda, cart=fit.cart, knn=fit.knn, svm=fit.svm, rf=fit.rf))
mat_Accuracy_PROVINCE[1,i] <- mean(results2$values[,2])
mat_Accuracy_PROVINCE[2,i] <- mean(results2$values[,4])
mat_Accuracy_PROVINCE[3,i] <- mean(results2$values[,6])
mat_Accuracy_PROVINCE[4,i] <- mean(results2$values[,8])
mat_Accuracy_PROVINCE[5,i] <- mean(results2$values[,10])

mat_Kappa_PROVINCE[1,i] <- mean(results2$values[,3])
mat_Kappa_PROVINCE[2,i] <- mean(results2$values[,5])
mat_Kappa_PROVINCE[3,i] <- mean(results2$values[,7])
mat_Kappa_PROVINCE[4,i] <- mean(results2$values[,9])
mat_Kappa_PROVINCE[5,i] <- mean(results2$values[,11])

# MAKE PREDICTIONS
predictions <- predict(fit.lda, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[1,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[1,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.cart, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[2,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[2,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.knn, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[3,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[3,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.svm, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[4,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[4,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.rf, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_PROVINCIA))
CM_Accuracy_PROVINCE[5,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PROVINCE[5,i] <- CM$overall[[2]] # Kappa

}
} # for j

names <- c('Jan','Feb','March','April','May','June','July','Aug','Sep','Oct','Nov','Dec','No MET')
colnames(mat_Accuracy_PROVINCE) <- names
names <- c("LDA","CART","KNN","SVM","RF")
rownames(mat_Accuracy_PROVINCE) <- names

names <- c('Jan','Feb','March','April','May','June','July','Aug','Sep','Oct','Nov','Dec','No MET')
colnames(mat_Kappa_PROVINCE) <- names
names <- c("LDA","CART","KNN","SVM","RF")
rownames(mat_Kappa_PROVINCE) <- names

####################################################################################
#################### ITERATIONS AVERAGE ############################
####################################################################################

mat_Accuracy_PROVINCE <- matrix(0,nrow=5,ncol=13)
mat_Kappa_PROVINCE <- matrix(0,nrow=5,ncol=13)
CM_Accuracy_PROVINCE <- matrix(0,nrow=5,ncol=13)
CM_Kappa_PROVINCE <- matrix(0,nrow=5,ncol=13)

in1 <- c(14,15,16)
in2 <- c(14,15,16)
for (k in 1:in1){
for (i in 1:5){
for (j in 1:13){
mat_Accuracy_PROVINCE[i,j] <- mean(list_CAL_ACC_PROVINCE[[k]][i,j])
mat_Kappa_PROVINCE[i,j] <- mean(list_CAL_KAP_PROVINCE[[k]][i,j])
}}}

for (k in 1:in2){
for (i in 1:5){
for (j in 1:13){
CM_Accuracy_PROVINCE[i,j] <- mean(list_VAL_ACC_PROVINCE[[k]][i,j])
CM_Kappa_PROVINCE[i,j] <- mean(list_VAL_KAP_PROVINCE[[k]][i,j])
}}}


########################################################################################################################################
############## CALIBRATION ####################################################################

par(mfcol=c(1,2));  
par(mar=c(5,5,4,2));  
par(oma=c(3,3,3,3));
  
#### ACCURACY FIGURE

par(new=FALSE)
pointsX <- c(0,0,0,0,0)
pointsY <- c(mat_Accuracy_PROVINCE[1,13],mat_Accuracy_PROVINCE[2,13],mat_Accuracy_PROVINCE[3,13],mat_Accuracy_PROVINCE[4,13],mat_Accuracy_PROVINCE[5,13])
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)
#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

par(new=TRUE)
pointsX <- c(0.5,0.5,0.5,0.5,0.5)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(mat_Accuracy_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)
#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(mat_Accuracy_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)
legend("bottomright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)


#### KAPPA FIGURE

par(new=FALSE)

pointsX <- c(0,0,0,0,0)
pointsY <- c(mat_Kappa_PROVINCE[1,13],mat_Kappa_PROVINCE[2,13],mat_Kappa_PROVINCE[3,13],mat_Kappa_PROVINCE[4,13],mat_Kappa_PROVINCE[5,13])
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

par(new=TRUE)
pointsX <- c(0.5,0.5,0.5,0.5,0.5)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(mat_Kappa_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(mat_Kappa_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)
legend("topright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)

#mtext("PROVINCE, Calibration",side=3,line = - 2, outer = TRUE, cex=1.5)


########################################################################################################################################

############## VALIDATION ####################################################################

par(mfcol=c(1,2));  
par(mar=c(5,5,4,2));  
par(oma=c(3,3,3,3));
  
#### ACCURACY FIGURE

par(new=FALSE)
pointsX <- c(0,0,0,0,0)
pointsY <- c(CM_Accuracy_PROVINCE[1,13],CM_Accuracy_PROVINCE[2,13],CM_Accuracy_PROVINCE[3,13],CM_Accuracy_PROVINCE[4,13],CM_Accuracy_PROVINCE[5,13])
pointsY
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

par(new=TRUE)
pointsX <- c(0.9,0.9,0.9,0.9,0.9)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(CM_Accuracy_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(CM_Accuracy_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

legend("bottomright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)


#### KAPPA FIGURE

par(new=FALSE)
pointsX <- c(0,0,0,0,0)
pointsY <- c(CM_Kappa_PROVINCE[1,13],CM_Kappa_PROVINCE[2,13],CM_Kappa_PROVINCE[3,13],CM_Kappa_PROVINCE[4,13],CM_Kappa_PROVINCE[5,13])
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)


par(new=TRUE)
pointsX <- c(0.9,0.9,0.9,0.9,0.9)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(CM_Kappa_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(CM_Kappa_PROVINCE[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

legend("bottomright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)

#mtext("PROVINCE, Validation",side=3,line = - 2, outer = TRUE, cex=1.5)

####################################################################

#########################################################################
# Figures 9 and 10
#########################################################################


#########################################################################################################################################################################################################
### ONLY NIR INFORMATION
#########################################################################################################################################################################################################

dtt <- read.table("C:/Data/INFO_CG_NIR_MET.txt", header=T)
attach(dtt)
dataset <- dtt
dataset[1,1:10]

######################### SELECTION OF PDOs (without: Baena, S. Cazorla, S. Mágina, S. Segura
ind <- c(0)
for (i in 1:222)
{
 if(dataset$C_DOP[i]==2 | dataset$C_DOP[i]==10 | dataset$C_DOP[i]==11 | dataset$C_DOP[i]==12)
 {
  ind <- cbind(ind,i)
 }
}
ind 
aux <- c(1)
ind2 <- sort(ind)[-aux]
ind2 ## PDOs off

###################################################################
### EVALUATE SOME ALGORITHMS
###################################################################

CG <- as.matrix(dtt[,8:10])
NIR <- as.matrix(dtt[,11:1247])

##################################### PCA
CG_NIR <- data.frame(CG,NIR)
COMP_PCA<- pcr(CG ~ NIR, ncomp=20, data=CG_NIR)
as.matrix(COMP_PCA$scores)
COMP_PCA$scores[-ind2,]
dim(COMP_PCA$scores[-ind2,])

dataset2 <- cbind(C_DOP=as.factor(dataset$C_DOP[-ind2]),as.matrix(COMP_PCA$scores[-ind2,]))
dataset2 <- as.data.frame(dataset2)
names <- c('C_DOP','PCA1','PCA2','PCA3','PCA4')
colnames(dataset2) <- names

# Run algorithms using 10-fold cross validation
control <- trainControl(method="cv", number=10)
metric <- "Accuracy"

#Let’s evaluate 5 different algorithms:
#Linear Discriminant Analysis (LDA)
#Classification and Regression Trees (CART).
#k-Nearest Neighbors (kNN).
#Support Vector Machines (SVM) with a linear kernel.
#Random Forest (RF)

###Create training/validation datasets
set.seed(j)
validation_index <- createDataPartition(dataset2$C_DOP[-ind2], p=0.80, list=FALSE)
validation <- dataset2[-validation_index,]
dataset2 <- dataset2[validation_index,]
dim(dataset2)
dim(validation)
as.factor(dataset2$C_DOP[-ind2])
as.factor(validation$C_DOP)
# a) linear algorithms
set.seed(j)
fit.lda <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="lda", metric=metric, trControl=control)
# b) nonlinear algorithms
# CART
set.seed(j)
fit.cart <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="rpart", metric=metric, trControl=control)
# kNN
set.seed(j)
fit.knn <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="knn", metric=metric, trControl=control)
# c) advanced algorithms
# SVM
set.seed(j)
fit.svm <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="svmRadial", metric=metric, trControl=control)
# Random Forest
set.seed(j)
fit.rf <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4,
                                data=dataset2, method="rf", metric=metric, trControl=control)
# summarize accuracy of models
results <- resamples(list(lda=fit.lda, cart=fit.cart, knn=fit.knn, svm=fit.svm, rf=fit.rf))
summary(results) results
results$values
ls(results$values)
mean(results$values[,2])
mean(results$values[,9])
dotplot(results,xlim=c(0,1)) #RF provides the best results
table(dataset$C_DOP[-ind2])

CM_Accuracy_PDO <- matrix(0,nrow=5,ncol=13)
CM_Kappa_PDO <- matrix(0,nrow=5,ncol=13)

# MAKE PREDICTIONS
# estimate skill of LDA on the validation dataset
predictions <- predict(fit.lda, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[1,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[1,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.cart, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[2,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[2,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.knn, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[3,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[3,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.svm, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[4,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[4,13] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.rf, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[5,13] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[5,13] <- CM$overall[[2]] # Kappa

#########################################################################################################################################################################################################
### NIR + MET INFORMATION
#########################################################################################################################################################################################################

mat_Accuracy_PDO <- matrix(0,nrow=5,ncol=13)
mat_Kappa_PDO <- matrix(0,nrow=5,ncol=13)

mat_Accuracy_PDO[1,13] <- mean(results$values[,2])
mat_Accuracy_PDO[2,13] <- mean(results$values[,4])
mat_Accuracy_PDO[3,13] <- mean(results$values[,6])
mat_Accuracy_PDO[4,13] <- mean(results$values[,8])
mat_Accuracy_PDO[5,13] <- mean(results$values[,10])

mat_Kappa_PDO[1,13] <- mean(results$values[,3])
mat_Kappa_PDO[2,13] <- mean(results$values[,5])
mat_Kappa_PDO[3,13] <- mean(results$values[,7])
mat_Kappa_PDO[4,13] <- mean(results$values[,9])
mat_Kappa_PDO[5,13] <- mean(results$values[,11])

for (i in 1:12)
{
dataset3 <- cbind(C_DOP=as.factor(dataset$C_DOP[-ind2]),as.matrix(COMP_PCA$scores[-ind2,]),as.matrix(lista_medida_MET[[i]][-ind2,]))
dataset3 <- as.data.frame(dataset3)
names <- c('C_DOP','PCA1','PCA2','PCA3','PCA4','PCA5','PCA6','PCA7','PCA8','PCA9','PCA10',
            'PCA11','PCA12','PCA13','PCA14','PCA15','PCA16','PCA17','PCA18','PCA19','PCA20',
            'Temp','Hum','WSped','WDic','Rad','Precip','Eto')
colnames(dataset3) <- names

# Run algorithms using 10-fold cross validation
control <- trainControl(method="cv", number=10)
metric <- "Accuracy"

#Let’s evaluate 5 different algorithms:
#Linear Discriminant Analysis (LDA)
#Classification and Regression Trees (CART).
#k-Nearest Neighbors (kNN).
#Support Vector Machines (SVM) with a linear kernel.
#Random Forest (RF)

###Create training/validation datasets
set.seed(j)
validation_index <- createDataPartition(dataset$C_DOP[-ind2], p=0.80, list=FALSE)
validation <- dataset3[-validation_index,]
dataset3 <- dataset3[validation_index,]
dim(dataset3)

# a) linear algorithms
set.seed(j)
fit.lda <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="lda", metric=metric, trControl=control)
# b) nonlinear algorithms
# CART
set.seed(j)
fit.cart <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="rpart", metric=metric, trControl=control)
# kNN
set.seed(j)
fit.knn <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="knn", metric=metric, trControl=control)
# c) advanced algorithms
# SVM
set.seed(j)
fit.svm <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="svmRadial", metric=metric, trControl=control)
# Random Forest
set.seed(j)
fit.rf <- train(as.factor(C_DOP)~PCA1+PCA2+PCA3+PCA4
                                +Temp+Hum+WSped+WDic+Rad+Precip+Eto,
                                data=dataset3, method="rf", metric=metric, trControl=control)
# summarize accuracy of models
results2 <- resamples(list(lda=fit.lda, cart=fit.cart, knn=fit.knn, svm=fit.svm, rf=fit.rf))
mat_Accuracy_PDO[1,i] <- mean(results2$values[,2])
mat_Accuracy_PDO[2,i] <- mean(results2$values[,4])
mat_Accuracy_PDO[3,i] <- mean(results2$values[,6])
mat_Accuracy_PDO[4,i] <- mean(results2$values[,8])
mat_Accuracy_PDO[5,i] <- mean(results2$values[,10])

mat_Kappa_PDO[1,i] <- mean(results2$values[,3])
mat_Kappa_PDO[2,i] <- mean(results2$values[,5])
mat_Kappa_PDO[3,i] <- mean(results2$values[,7])
mat_Kappa_PDO[4,i] <- mean(results2$values[,9])
mat_Kappa_PDO[5,i] <- mean(results2$values[,11])

# MAKE PREDICTIONS
predictions <- predict(fit.lda, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[1,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[1,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.cart, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[2,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[2,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.knn, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[3,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[3,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.svm, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[4,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[4,i] <- CM$overall[[2]] # Kappa

predictions <- predict(fit.rf, validation)
CM <- confusionMatrix(predictions, as.factor(validation$C_DOP[-ind2]))
CM_Accuracy_PDO[5,i] <- CM$overall[[1]] # Accuracy
CM_Kappa_PDO[5,i] <- CM$overall[[2]] # Kappa

}
} #for j

names <- c('Jan','Feb','March','April','May','June','July','Aug','Sep','Oct','Nov','Dec','No MET')
colnames(mat_Accuracy_PDO) <- names
names <- c("LDA","CART","KNN","SVM","RF")
rownames(mat_Accuracy_PDO) <- names

names <- c('Jan','Feb','March','April','May','June','July','Aug','Sep','Oct','Nov','Dec','No MET')
colnames(mat_Kappa_PDO) <- names
names <- c("LDA","CART","KNN","SVM","RF")
rownames(mat_Kappa_PDO) <- names

####################################################################################
#################### ITERATIONS AVERAGE ############################
####################################################################################

mat_Accuracy_PDO <- matrix(0,nrow=5,ncol=13)
mat_Kappa_PDO <- matrix(0,nrow=5,ncol=13)
CM_Accuracy_PDO <- matrix(0,nrow=5,ncol=13)
CM_Kappa_PDO <- matrix(0,nrow=5,ncol=13)

in1 <- c(14,15,16)
in2 <- c(14,15,16)
for (k in 1:in1){
for (i in 1:5){
for (j in 1:13){
mat_Accuracy_PDO[i,j] <- mean(list_CAL_ACC_PDO[[k]][i,j])
mat_Kappa_PDO[i,j] <- mean(list_CAL_KAP_PDO[[k]][i,j])
}}}

for (k in 1:in2){
for (i in 1:5){
for (j in 1:13){
CM_Accuracy_PDO[i,j] <- mean(list_VAL_ACC_PDO[[k]][i,j])
CM_Kappa_PDO[i,j] <- mean(list_VAL_KAP_PDO[[k]][i,j])
}}}

########################################################################################################################################
############## CALIBRATION ####################################################################

par(mfcol=c(1,2));  
par(mar=c(5,5,4,2));  
par(oma=c(3,3,3,3));
  
#### ACCURACY FIGURE

par(new=FALSE)
pointsX <- c(0,0,0,0,0)
pointsY <- c(mat_Accuracy_PDO[1,13],mat_Accuracy_PDO[2,13],mat_Accuracy_PDO[3,13],mat_Accuracy_PDO[4,13],mat_Accuracy_PDO[5,13])
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)
#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

par(new=TRUE)
pointsX <- c(0.5,0.5,0.5,0.5,0.5)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(mat_Accuracy_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)
#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(mat_Accuracy_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)
legend("bottomright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)


#### KAPPA FIGURE

par(new=FALSE)

pointsX <- c(0,0,0,0,0)
pointsY <- c(mat_Kappa_PDO[1,13],mat_Kappa_PDO[2,13],mat_Kappa_PDO[3,13],mat_Kappa_PDO[4,13],mat_Kappa_PDO[5,13])
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

par(new=TRUE)
pointsX <- c(0.5,0.5,0.5,0.5,0.5)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(mat_Kappa_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(mat_Kappa_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)
legend("topright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)

#mtext("PDO, Calibration",side=3,line = - 2, outer = TRUE, cex=1.5)


########################################################################################################################################

############## VALIDATION ####################################################################

par(mfcol=c(1,2));  
par(mar=c(5,5,4,2));  
par(oma=c(3,3,3,3));
  
#### ACCURACY FIGURE

par(new=FALSE)
pointsX <- c(0,0,0,0,0)
pointsY <- c(CM_Accuracy_PDO[1,13],CM_Accuracy_PDO[2,13],CM_Accuracy_PDO[3,13],CM_Accuracy_PDO[4,13],CM_Accuracy_PDO[5,13])
pointsY
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Accuracy",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

par(new=TRUE)
pointsX <- c(0.9,0.9,0.9,0.9,0.9)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(CM_Accuracy_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(CM_Accuracy_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Accuracy",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

legend("topright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)


#### KAPPA FIGURE

par(new=FALSE)
pointsX <- c(0,0,0,0,0)
pointsY <- c(CM_Kappa_PDO[1,13],CM_Kappa_PDO[2,13],CM_Kappa_PDO[3,13],CM_Kappa_PDO[4,13],CM_Kappa_PDO[5,13])
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)

#Para colorear partes del gráfico:
xx1<-c(-0.4,-0.4,0.5,0.5)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col="antiquewhite",border="antiquewhite")
par(new=TRUE)
plot(pointsX,pointsY,type="p",lwd=1,xlim=range(0:12),ylim=range(0:1),las=1,xlab="Month",ylab="Kappa",col=c("chartreuse4","magenta","orange","blue","red"),xaxt = "n",cex=1.7,pch=19,cex.lab=1.7,las=1)


par(new=TRUE)
pointsX <- c(0.9,0.9,0.9,0.9,0.9)
pointsY <- c(0,0.25,0.5,0.75,1)
points <- data.frame(pointsX,pointsY)
plot(pointsX,pointsY,type="l",col="grey",lwd=3,lty=c(1),pch=25,xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),cex.lab=1.7,las=1)

par(new=TRUE)
matplot(t(m_axisX),t(CM_Kappa_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

#Para colorear partes del gráfico:
xx1<-c(0.5,0.5,12.4,12.4)
yy1<-c(0,1,1,0)
#polygon(xx1,yy1,col="gray94",border="gray92")
polygon(xx1,yy1,col=palette("Pastel 2"),border=palette("Pastel 2"))
par(new=TRUE)
matplot(t(m_axisX),t(CM_Kappa_PDO[,1:12]),type="b",lwd=3,pch=19,cex=1.5,lty=c(2:5,1),col=c("chartreuse4","magenta","orange","blue","red"),xlab="Month",ylab="Kappa",xaxt = "n",ylim=c(0:1),xlim=range(0:12),las=1,cex.lab=1.7)

legend("topright",pch=19,col=c("chartreuse4","magenta","orange","blue","red"),c("LDA","CART","KNN","SVM","RF"),cex=1.5,bg="white")
axis(1,at=0:12,labels=c("No MET","Jan","Feb","March","April","May","June","July","Aug","Sep","Oct","Nov","Dec"),cex.axis=1.1,cex.lab=1.5)

#mtext("PDO, Validation",side=3,line = - 2, outer = TRUE, cex=1.5)

####################################################################