1 #setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno")
4 # Lecture des données: pm = dataframe 2 colonnes, date-time puis PM10 horaire
5 pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=",",dec=".",header=T)
7 #datedebut = "10/12/2008"
9 # Chargement des données météo et indicateurs: VarExp matrice des données météo,
10 # première colonne = date, première rangée = second jour
11 VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=",",dec=".",header=T)
14 # Lecture des dates: ??? sans doute liste des jours de 10/12/2008 à 06/2016...
15 #dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T)
16 dates = VarExp[,1] #dates[,1]
18 # Contruction des matrices de données: pm.h = matrice des séries en ligne
19 pm.h <- matrix(pm[,2],ncol=24,byrow=TRUE)
22 # Data = matrice des couples de jours (séries de 48 PM10),
23 # première ligne = jour 1,2, dernière = jour N-1,N
24 Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ])
26 # dates2 = dates du 2eme jour au dernier
27 dates2 = dates[2:Nlignes]
29 rownames(Data) = dates2
30 # df contient l'ensemble des données.
31 #df <- cbind(Data,varexp[,-1])
34 # Complétion des variables exogènes (du 2eme jour) par les PM10 moyen ce même jour
35 PMjour <- apply(df[,25:48],1,mean,na.rm=T)
36 dfexp <- cbind(VarExp,PMjour)
38 Dates = c( #Month/Day/Year, as in meteo file
43 Categorie = c("Epandage", "Chauffage", "Non Polluée")
45 RepFig = "FIGURES_Etude"
50 j=1 # numéro de semaine
51 ij=6 # numéro du jour (0 = lundi)
59 Hc = 7 #predict at 7:00, from 8:00
60 H=24+Hc #Hc: dernier PM10 connu avant prédiction des 24-Hc restants, le jour 2
63 # Premier conditionnement : mois
64 indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10
65 data = df[indcond,] #restriction aux couples dont le mois du 2eme jour est 2,3,4,9 ou 10
66 varexp = dfexp[indcond,] #de même sur les exogènes (+ PM10 moyens)
68 nl = (1:nrow(data))[rownames(data)==Dates[j]] #numéro de ligne où 2eme jour == 16/03/2015
69 dateJPrev = rownames(data)[nl+ij] # translation : date du jour à prévoir (2eme colonne)
70 dataj = as.numeric(data[nl+ij, 1:48]) #extraction de cette série dans dataj
71 data = data[-(nl + ij), ] #suppression de cette série dans data
72 varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes
73 indNA = attr(na.omit(data[, 1:48]),"na.action")
75 data = data[-indNA,] #remove all pairs of days with NAs
76 varexp = varexp[-indNA,] #and remove corresponding exogenous variables
78 # Second conditionnement : les jours avec PMjour +/- large
80 bornes = mean(dataj[25:48]) + c(-large,large)
81 indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
82 data = data[indcond,] #pollution du 2eme jour == pollution du jour courant +/- 1
83 varexp = varexp[indcond,]
86 for (k in 1:nrow(data))
88 #D[k] = sqrt(sum((1:H)*(dataj[L] - data[k,L])^2))
89 D[k] = sqrt(sum((dataj[L] - data[k,L])^2))
91 ind = order(D)[1:nbvois]
95 JourMoy = apply(data[ind, 1:48], 2, mean)
96 #JourMoy = apply(W*data[ind, 1:48], 2, sum)
97 NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="")
98 Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="")
99 #erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2))
103 erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))
105 # erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48]))
107 # matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind),
108 # cex.axis=1.4, cex.main = 1.7, cex.lab=1.5,
109 # xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1),
110 # " / ",round(erreurPrev,1)), main=Titre)
111 # legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind))
112 # lines(1:48, dataj, lwd=2.5)
113 # lines(JourMoy, lty = 2, lwd=2)
114 # abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2)
117 # Err24 = c(Err24, erreur24)
118 # ErrPrev = c(ErrPrev, erreurPrev)
119 # ResDates = cbind(ResDates, rownames(data)[ind])
122 #rownames(ResDates) = 1:10
125 #for (Col in ncol(ResDates):1)
132 # if (ResDates[I,1] == ResDates[J,Col])
136 # Kvois = c(Kvois, K)
139 ## pdf("Erreur_Epandage_PMjour.pdf")
140 #ymin = min(na.omit(ErrPrev), Err24)
141 #ymin = min(ymin, Kvois)
142 #ymax = max(na.omit(ErrPrev), Err24)
143 #ymax = max(ymax, Kvois)
144 #plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5,
145 # ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax))
146 #lines(5:24,ErrPrev, lwd=2, lty = 2)
147 #legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5)
148 #points(5:24, Kvois, cex=1.8, pch = 19)