#ponizej znajduje sie przykladowa analiza EOF wraz z regresją wielokrotną
#PC SLP do temeraptury powietrza

#1
#zaladowanie danych z Reanalizy NCEP NCAR

slp = slp.NCEP(lon = c(-40,40), lat = c(35,75))

#przykladowe wykreslenie dla Stycznia

map(slp, it='Jul')

#2
#obliczenie EOF dla stycznie

eof_slp_7 = EOF(slp, it='Jul', verbose = TRUE)

#wykresy
plot(eof_slp_7, new = FALSE, pattern = 1)

#wykreslenie mapy dla wybranego EOF

map(eof_slp_1, pattern = 1, new = FALSE)

# 3 
#przygotowanie danych (serii czasowych)
#zlozenie danych PC z okreslonej liczby EOFow

#ponizej nalezy wpisac liczbe EOF ktora nas interesuje

n = 8

PCs = eof_slp_7[ ,1:8]
PCs = as.data.frame(PCs)
PCs = cbind(ROK = c(1948:2012), PCs)

# 4
#wczytanie danych odnosnie temperatury powietrza w Europie Srodkowej
temp = read.table("pc_reg_temp.txt", header = TRUE)

#wybieramy tylko styczen
temp_7 = temp[which(temp$MC == 7), ]

#5
#teraz trzeba przyciac dane do wspolnego okresu i polacyc w jedna data.frame
#wspolny okres to 1951-2012 tak wiec najpier wybierzemy te dane z PCs


PCs_5112 = PCs[which(PCs$ROK >= 1951 & PCs$ROK <= 2012), ]

#6
#teraz temperatura

temp_7_5112 = temp_7[which(temp_7$ROK >= 1951 & temp_7$ROK <= 2012), ]

#7
#gotowe, a teraz polaczymy nasze dane w jedna dataframe

data = cbind(PCs_5112[ , ], temp_7_5112[ , 3])

#zmienmy sobie r?wniez nazwy kolumn
colnames(data) = c("ROK", "PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "TEMP")
write.table(data, "dane_reg.txt", sep = "\t", row.names = FALSE)
#teraz podlaczzmy dane w celu latwiejszego z nich korzystania


attach(data)

#8
#model z PC1

model1 = lm(TEMP~PC1)

summary(model1)


plot(TEMP, model1$fitted.values, ylab = "TEMP (model)")
abline(0,1, col = 2, lwd = 2)

qqplot(TEMP, model1$fitted.values)
abline(0,1, col = 2, lwd = 2)

plot(1951:2012, TEMP, type = "l", xlab = "ROK")
lines(1951:2012, model1$fitted.values, col = 2, lwd = 2)

#8
#teraz model ze wszystkimmi wybranymi PCs


model2 = lm(TEMP~PC1+PC2+PC3+PC4+PC5+PC6)
summary(model2)


model3 = lm(TEMP~PC1+PC3+PC5+PC6)
summary(model3)





plot(1951:2012, TEMP, type = "l", xlab = "ROK")
lines(1951:2012, model1$fitted.values, col = 4, lwd = 2)
lines(1951:2012, model2$fitted.values, col = 2, lwd = 2)
lines(1951:2012, model3$fitted.values, col = 3, lwd = 2)



qqplot(TEMP, model3$fitted.values)
abline(0,1, col = 2, lwd = 2)

AIC(model1, model2, model3)

#9
# pokazanie wyników modeli na diagramie Taylora
#niezbędna jest paczka plotrix

taylor.diagram(TEMP, model1$fitted.values, add = FALSE, col = 2, sd.arcs = TRUE, ref.sd = TRUE)
taylor.diagram(TEMP, model2$fitted.values, add = TRUE, col = 3, sd.arcs = TRUE, ref.sd = TRUE)
taylor.diagram(TEMP, model3$fitted.values, add = TRUE, col = 4, sd.arcs = TRUE, ref.sd = TRUE)


detach(data)