Commit | Line | Data |
---|---|---|
9dbd32d9 | 1 | #setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno") |
d2ab47a7 BA |
2 | rm(list=ls()) |
3 | ||
9dbd32d9 | 4 | # Lecture des données: pm = dataframe 2 colonnes, date-time puis PM10 horaire |
794d99e9 | 5 | pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=",",dec=".",header=T) |
9dbd32d9 BA |
6 | #n = dim(pm)[1] |
7 | #datedebut = "10/12/2008" | |
d2ab47a7 | 8 | |
9dbd32d9 BA |
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 | |
794d99e9 | 11 | VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=",",dec=".",header=T) |
d2ab47a7 BA |
12 | VarExp <- VarExp[-1,] |
13 | ||
9dbd32d9 | 14 | # Lecture des dates: ??? sans doute liste des jours de 10/12/2008 à 06/2016... |
794d99e9 BA |
15 | #dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T) |
16 | dates = VarExp[,1] #dates[,1] | |
d2ab47a7 | 17 | |
9dbd32d9 | 18 | # Contruction des matrices de données: pm.h = matrice des séries en ligne |
d2ab47a7 | 19 | pm.h <- matrix(pm[,2],ncol=24,byrow=TRUE) |
d2ab47a7 | 20 | Nlignes = nrow(pm.h) |
9dbd32d9 BA |
21 | |
22 | # Data = matrice des couples de jours (séries de 48 PM10), | |
794d99e9 | 23 | # première ligne = jour 1,2, dernière = jour N-1,N |
d2ab47a7 | 24 | Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ]) |
9dbd32d9 BA |
25 | |
26 | # dates2 = dates du 2eme jour au dernier | |
d2ab47a7 | 27 | dates2 = dates[2:Nlignes] |
9dbd32d9 | 28 | |
d2ab47a7 BA |
29 | rownames(Data) = dates2 |
30 | # df contient l'ensemble des données. | |
31 | #df <- cbind(Data,varexp[,-1]) | |
32 | df <- Data | |
9dbd32d9 BA |
33 | |
34 | # Complétion des variables exogènes (du 2eme jour) par les PM10 moyen ce même jour | |
d2ab47a7 BA |
35 | PMjour <- apply(df[,25:48],1,mean,na.rm=T) |
36 | dfexp <- cbind(VarExp,PMjour) | |
37 | ||
794d99e9 BA |
38 | Dates = c( #Month/Day/Year, as in meteo file |
39 | "3/16/2015", | |
40 | "1/19/2015", | |
41 | "4/27/2015") | |
d2ab47a7 BA |
42 | |
43 | Categorie = c("Epandage", "Chauffage", "Non Polluée") | |
44 | ||
45 | RepFig = "FIGURES_Etude" | |
46 | ||
47 | ResDates = NULL | |
48 | ||
49 | nbvois=10 | |
50 | j=1 # numéro de semaine | |
51 | ij=6 # numéro du jour (0 = lundi) | |
52 | ||
53 | Err24 = NULL | |
54 | ErrPrev = NULL | |
55 | Kvois = NULL | |
56 | ||
794d99e9 BA |
57 | #for (Hc in 5:24) |
58 | #{ | |
59 | Hc = 7 #predict at 7:00, from 8:00 | |
9dbd32d9 BA |
60 | H=24+Hc #Hc: dernier PM10 connu avant prédiction des 24-Hc restants, le jour 2 |
61 | L = 1:H | |
62 | ||
63 | # Premier conditionnement : mois | |
794d99e9 | 64 | indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10 |
9dbd32d9 BA |
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) | |
67 | ||
68 | nl = (1:nrow(data))[rownames(data)==Dates[j]] #numéro de ligne où 2eme jour == 16/03/2015 | |
794d99e9 | 69 | dateJPrev = rownames(data)[nl+ij] # translation : date du jour à prévoir (2eme colonne) |
9dbd32d9 | 70 | dataj = as.numeric(data[nl+ij, 1:48]) #extraction de cette série dans dataj |
794d99e9 BA |
71 | data = data[-(nl + ij), ] #suppression de cette série dans data |
72 | varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes | |
9dbd32d9 | 73 | indNA = attr(na.omit(data[, 1:48]),"na.action") |
794d99e9 | 74 | # print(indNA) |
adf76c2a BA |
75 | data = data[-indNA,] #remove all pairs of days with NAs |
76 | varexp = varexp[-indNA,] #and remove corresponding exogenous variables | |
9dbd32d9 | 77 | |
adf76c2a | 78 | # Second conditionnement : les jours avec PMjour +/- large |
9dbd32d9 | 79 | large = 1 |
adf76c2a | 80 | bornes = mean(dataj[25:48]) + c(-large,large) |
9dbd32d9 | 81 | indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2] |
adf76c2a | 82 | data = data[indcond,] #pollution du 2eme jour == pollution du jour courant +/- 1 |
9dbd32d9 BA |
83 | varexp = varexp[indcond,] |
84 | ||
85 | D = rep(0,nrow(data)) | |
86 | for (k in 1:nrow(data)) | |
87 | { | |
88 | #D[k] = sqrt(sum((1:H)*(dataj[L] - data[k,L])^2)) | |
89 | D[k] = sqrt(sum((dataj[L] - data[k,L])^2)) | |
90 | } | |
91 | ind = order(D)[1:nbvois] | |
794d99e9 BA |
92 | # w = 1/(D[ind]^2) |
93 | # w = w/sum(w) | |
94 | # W = w %o% rep(1,48) | |
9dbd32d9 BA |
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)) | |
794d99e9 | 100 | if(Hc==24) { |
adf76c2a | 101 | erreurPrev = NA |
794d99e9 | 102 | } else { |
9dbd32d9 | 103 | erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48])) |
9dbd32d9 | 104 | } |
794d99e9 BA |
105 | # erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48])) |
106 | # #png(NomFile) | |
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) | |
115 | # #xx=dev.off() | |
116 | # | |
117 | # Err24 = c(Err24, erreur24) | |
118 | # ErrPrev = c(ErrPrev, erreurPrev) | |
119 | # ResDates = cbind(ResDates, rownames(data)[ind]) | |
120 | #} | |
121 | ||
122 | #rownames(ResDates) = 1:10 | |
123 | # | |
124 | #Kvois = NULL | |
125 | #for (Col in ncol(ResDates):1) | |
126 | #{ | |
127 | # K = 0 | |
128 | # for (I in 1:10) | |
129 | # { | |
130 | # for (J in 1:10) | |
131 | # { | |
132 | # if (ResDates[I,1] == ResDates[J,Col]) | |
133 | # K = K +1 | |
134 | # } | |
135 | # } | |
136 | # Kvois = c(Kvois, K) | |
137 | #} | |
138 | # | |
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) | |
149 | ## xx = dev.off() | |
150 | # | |
151 | ##length(D) |