rm(list=ls())
# Lecture des données: pm = dataframe 2 colonnes, date-time puis PM10 horaire
-pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=";",dec=".",header=T)
+pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=",",dec=".",header=T)
#n = dim(pm)[1]
#datedebut = "10/12/2008"
# Chargement des données météo et indicateurs: VarExp matrice des données météo,
# première colonne = date, première rangée = second jour
-VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=";",dec=".",header=T)
+VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=",",dec=".",header=T)
VarExp <- VarExp[-1,]
# Lecture des dates: ??? sans doute liste des jours de 10/12/2008 à 06/2016...
-dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T)
-dates = dates[,1]
+#dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T)
+dates = VarExp[,1] #dates[,1]
# Contruction des matrices de données: pm.h = matrice des séries en ligne
pm.h <- matrix(pm[,2],ncol=24,byrow=TRUE)
Nlignes = nrow(pm.h)
# Data = matrice des couples de jours (séries de 48 PM10),
-# première ligne = jour 1/2, dernière = jour N-1/N
+# première ligne = jour 1,2, dernière = jour N-1,N
Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ])
# dates2 = dates du 2eme jour au dernier
PMjour <- apply(df[,25:48],1,mean,na.rm=T)
dfexp <- cbind(VarExp,PMjour)
-Dates = c(
-"16/03/2015",
-"19/01/2015",
-"27/04/2015")
+Dates = c( #Month/Day/Year, as in meteo file
+"3/16/2015",
+"1/19/2015",
+"4/27/2015")
Categorie = c("Epandage", "Chauffage", "Non Polluée")
ErrPrev = NULL
Kvois = NULL
-for (Hc in 5:24)
-{
+#for (Hc in 5:24)
+#{
+ Hc = 7 #predict at 7:00, from 8:00
H=24+Hc #Hc: dernier PM10 connu avant prédiction des 24-Hc restants, le jour 2
L = 1:H
# Premier conditionnement : mois
- indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4
- | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10
+ indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10
data = df[indcond,] #restriction aux couples dont le mois du 2eme jour est 2,3,4,9 ou 10
varexp = dfexp[indcond,] #de même sur les exogènes (+ PM10 moyens)
nl = (1:nrow(data))[rownames(data)==Dates[j]] #numéro de ligne où 2eme jour == 16/03/2015
- dateJPrev = rownames(data)[nl+ij] # ??? la date du 1er jour du couple de jours pile une semaine après le couple 'nl'
+ dateJPrev = rownames(data)[nl+ij] # translation : date du jour à prévoir (2eme colonne)
dataj = as.numeric(data[nl+ij, 1:48]) #extraction de cette série dans dataj
- data = data[-(nl + ij), ] #suppression de cette série dans data (????)
- varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes (date 16/03/2015 + 6 jours
+ data = data[-(nl + ij), ] #suppression de cette série dans data
+ varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes
indNA = attr(na.omit(data[, 1:48]),"na.action")
+# print(indNA)
data = data[-indNA,] #remove all pairs of days with NAs
varexp = varexp[-indNA,] #and remove corresponding exogenous variables
D[k] = sqrt(sum((dataj[L] - data[k,L])^2))
}
ind = order(D)[1:nbvois]
- w = 1/(D[ind]^2)
- w = w/sum(w)
- W = w %o% rep(1,48)
+# w = 1/(D[ind]^2)
+# w = w/sum(w)
+# W = w %o% rep(1,48)
JourMoy = apply(data[ind, 1:48], 2, mean)
#JourMoy = apply(W*data[ind, 1:48], 2, sum)
NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="")
Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="")
#erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2))
- if(Hc==24)
+ if(Hc==24) {
erreurPrev = NA
- else
+ } else {
erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))
- erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48]))
- #png(NomFile)
- matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind),
- cex.axis=1.4, cex.main = 1.7, cex.lab=1.5,
- xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1),
- " / ",round(erreurPrev,1)), main=Titre)
- legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind))
- lines(1:48, dataj, lwd=2.5)
- lines(JourMoy, lty = 2, lwd=2)
- abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2)
- #xx=dev.off()
-
- Err24 = c(Err24, erreur24)
- ErrPrev = c(ErrPrev, erreurPrev)
- ResDates = cbind(ResDates, rownames(data)[ind])
-}
-
-rownames(ResDates) = 1:10
-
-Kvois = NULL
-for (Col in ncol(ResDates):1)
-{
- K = 0
- for (I in 1:10)
- {
- for (J in 1:10)
- {
- if (ResDates[I,1] == ResDates[J,Col])
- K = K +1
- }
}
- Kvois = c(Kvois, K)
-}
-
-# pdf("Erreur_Epandage_PMjour.pdf")
-ymin = min(na.omit(ErrPrev), Err24)
-ymin = min(ymin, Kvois)
-ymax = max(na.omit(ErrPrev), Err24)
-ymax = max(ymax, Kvois)
-plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5,
- ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax))
-lines(5:24,ErrPrev, lwd=2, lty = 2)
-legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5)
-points(5:24, Kvois, cex=1.8, pch = 19)
-# xx = dev.off()
-
-#length(D)
+# erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48]))
+# #png(NomFile)
+# matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind),
+# cex.axis=1.4, cex.main = 1.7, cex.lab=1.5,
+# xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1),
+# " / ",round(erreurPrev,1)), main=Titre)
+# legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind))
+# lines(1:48, dataj, lwd=2.5)
+# lines(JourMoy, lty = 2, lwd=2)
+# abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2)
+# #xx=dev.off()
+#
+# Err24 = c(Err24, erreur24)
+# ErrPrev = c(ErrPrev, erreurPrev)
+# ResDates = cbind(ResDates, rownames(data)[ind])
+#}
+
+#rownames(ResDates) = 1:10
+#
+#Kvois = NULL
+#for (Col in ncol(ResDates):1)
+#{
+# K = 0
+# for (I in 1:10)
+# {
+# for (J in 1:10)
+# {
+# if (ResDates[I,1] == ResDates[J,Col])
+# K = K +1
+# }
+# }
+# Kvois = c(Kvois, K)
+#}
+#
+## pdf("Erreur_Epandage_PMjour.pdf")
+#ymin = min(na.omit(ErrPrev), Err24)
+#ymin = min(ymin, Kvois)
+#ymax = max(na.omit(ErrPrev), Err24)
+#ymax = max(ymax, Kvois)
+#plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5,
+# ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax))
+#lines(5:24,ErrPrev, lwd=2, lty = 2)
+#legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5)
+#points(5:24, Kvois, cex=1.8, pch = 19)
+## xx = dev.off()
+#
+##length(D)
projects / talweg.git / blobdiff