Commit | Line | Data |
---|---|---|
d2ab47a7 BA |
1 | setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno") |
2 | rm(list=ls()) | |
3 | ||
4 | # Lecture des données | |
5 | pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=";",dec=".",header=T) #,row.names=1) | |
6 | n = dim(pm)[1] | |
7 | datedebut = "10/12/2008" | |
8 | ||
9 | # Chargement des données météo et indicateurs | |
10 | VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=";",dec=".",header=T) | |
11 | VarExp <- VarExp[-1,] | |
12 | ||
13 | # Lecture des dates | |
14 | dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T) | |
15 | dates = dates[,1] | |
16 | ||
17 | # Contruction des matrices de données | |
18 | pm.h <- matrix(pm[,2],ncol=24,byrow=TRUE) | |
19 | ||
20 | Nlignes = nrow(pm.h) | |
21 | Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ]) | |
22 | dates2 = dates[2:Nlignes] | |
23 | rownames(Data) = dates2 | |
24 | # df contient l'ensemble des données. | |
25 | #df <- cbind(Data,varexp[,-1]) | |
26 | df <- Data | |
27 | PMjour <- apply(df[,25:48],1,mean,na.rm=T) | |
28 | dfexp <- cbind(VarExp,PMjour) | |
29 | ||
30 | Dates = c( | |
31 | "16/03/2015", | |
32 | "19/01/2015", | |
33 | "27/04/2015") | |
34 | ||
35 | Categorie = c("Epandage", "Chauffage", "Non Polluée") | |
36 | ||
37 | RepFig = "FIGURES_Etude" | |
38 | ||
39 | ResDates = NULL | |
40 | ||
41 | nbvois=10 | |
42 | j=1 # numéro de semaine | |
43 | ij=6 # numéro du jour (0 = lundi) | |
44 | ||
45 | Err24 = NULL | |
46 | ErrPrev = NULL | |
47 | Kvois = NULL | |
48 | ||
49 | for (Hc in 5:24) { | |
50 | ||
51 | H=24+Hc | |
52 | fen=H | |
53 | #fen=12 ou 24 ou H | |
54 | L = (H-fen+1):H | |
55 | ||
56 | # Premier conditionnement : mois | |
57 | indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10 | |
58 | data = df[indcond,] | |
59 | varexp = dfexp[indcond,] | |
60 | ||
61 | nl = (1:nrow(data))[rownames(data)==Dates[j]] | |
62 | dateJPrev = rownames(data)[nl+ij] | |
63 | dataj = as.numeric(data[nl+ij, 1:48]) | |
64 | data = data[-(nl + ij), ] | |
65 | varexp = varexp[-(nl + ij), ] | |
66 | indNA = attr(na.omit(data[, 1:48]),"na.action") | |
67 | data = data[-indNA,] | |
68 | varexp = varexp[-indNA,] | |
69 | ||
70 | # Conditionnement : les jours avec PMjour +/- large | |
71 | large = 1 | |
72 | bornes = mean(dataj[25:48])+c(-large,large) | |
73 | indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2] | |
74 | data = data[indcond,] | |
75 | varexp = varexp[indcond,] | |
76 | ||
77 | D = rep(0,nrow(data)) | |
78 | for (k in 1:nrow(data)) { | |
79 | #D[k] = sqrt(sum((1:H)*(dataj[L] -data[k,L])^2)) | |
80 | D[k] = sqrt(sum((dataj[L] -data[k,L])^2)) | |
81 | } | |
82 | ind = order(D)[1:nbvois] | |
83 | w = 1/(D[ind]^2) | |
84 | w = w/sum(w) | |
85 | W = w%o%rep(1,48) | |
86 | JourMoy = apply(data[ind, 1:48], 2, mean) | |
87 | #JourMoy = apply(W*data[ind, 1:48], 2, sum) | |
88 | NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="") | |
89 | Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="") | |
90 | #erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2)) | |
91 | if(Hc==24){erreurPrev=NA}else{erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))} | |
92 | erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48])) | |
93 | png(NomFile) | |
94 | matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind), | |
95 | cex.axis=1.4, cex.main = 1.7, cex.lab=1.5, | |
96 | xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1)," / ",round(erreurPrev,1)), main=Titre) | |
97 | legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind)) | |
98 | lines(1:48, dataj, lwd=2.5) | |
99 | lines(JourMoy, lty = 2, lwd=2) | |
100 | abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2) | |
101 | xx=dev.off() | |
102 | ||
103 | ||
104 | ||
105 | Err24 = c(Err24, erreur24) | |
106 | ErrPrev = c(ErrPrev, erreurPrev) | |
107 | ResDates = cbind(ResDates, rownames(data)[ind]) | |
108 | } | |
109 | ||
110 | ||
111 | rownames(ResDates) = 1:10 | |
112 | ||
113 | Kvois = NULL | |
114 | for (Col in ncol(ResDates):1) { | |
115 | K = 0 | |
116 | for (I in 1:10){ | |
117 | for (J in 1:10){ | |
118 | if (ResDates[I,1] == ResDates[J,Col]) K = K +1 | |
119 | } | |
120 | } | |
121 | Kvois = c(Kvois, K) | |
122 | } | |
123 | ||
124 | # pdf("Erreur_Epandage_PMjour.pdf") | |
125 | ymin = min(na.omit(ErrPrev), Err24) | |
126 | ymin = min(ymin, Kvois) | |
127 | ymax = max(na.omit(ErrPrev), Err24) | |
128 | ymax = max(ymax, Kvois) | |
129 | plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5, | |
130 | ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax)) | |
131 | lines(5:24,ErrPrev, lwd=2, lty = 2) | |
132 | legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5) | |
133 | points(5:24, Kvois, cex=1.8, pch = 19) | |
134 | # xx = dev.off() | |
135 | ||
136 | length(D) |