Etude_En_Cours_R_E_bis.R

   1 #setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno")
   2 rm(list=ls())
   3
   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)
   6 #n = dim(pm)[1]
   7 #datedebut = "10/12/2008"
   8
   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)
  12 VarExp <- VarExp[-1,]
  13
  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 = dates[,1]
  17
  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)
  20 Nlignes = nrow(pm.h)
  21
  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, ])
  25
  26 # dates2 = dates du 2eme jour au dernier
  27 dates2 = dates[2:Nlignes]
  28
  29 rownames(Data) = dates2
  30 # df contient l'ensemble des données.
  31 #df <- cbind(Data,varexp[,-1])
  32 df <- Data
  33
  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)
  37
  38 Dates = c(
  39 "16/03/2015",
  40 "19/01/2015",
  41 "27/04/2015")
  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
  57 for (Hc in 5:24)
  58 {
  59     H=24+Hc #Hc: dernier PM10 connu avant prédiction des 24-Hc restants, le jour 2
  60     L = 1:H
  61
  62     # Premier conditionnement : mois
  63     indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4
  64         | 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)
  67
  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] # ??? la date du 1er jour du couple de jours pile une semaine après le couple 'nl'
  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 (date 16/03/2015 + 6 jours
  73     indNA = attr(na.omit(data[, 1:48]),"na.action")
  74     data = data[-indNA,]
  75     varexp = varexp[-indNA,]
  76
  77     # Conditionnement : les jours avec PMjour +/- large
  78     large = 1
  79     bornes = mean(dataj[25:48])+c(-large,large)
  80     indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
  81     data = data[indcond,]
  82     varexp = varexp[indcond,]
  83
  84     D = rep(0,nrow(data))
  85     for (k in 1:nrow(data))
  86     {
  87         #D[k] = sqrt(sum((1:H)*(dataj[L] - data[k,L])^2))
  88         D[k] = sqrt(sum((dataj[L] - data[k,L])^2))
  89     }
  90     ind = order(D)[1:nbvois]
  91     w = 1/(D[ind]^2)
  92     w = w/sum(w)
  93     W = w %o% rep(1,48)
  94     JourMoy = apply(data[ind, 1:48], 2, mean)
  95     #JourMoy = apply(W*data[ind, 1:48], 2, sum)
  96     NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="")
  97     Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="")
  98     #erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2))
  99     if(Hc==24)
 100         erreurPrev=NA
 101     else
 102         erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))
 103     erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48]))
 104     #png(NomFile)
 105     matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind),
 106         cex.axis=1.4, cex.main = 1.7, cex.lab=1.5,
 107         xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1),
 108             " / ",round(erreurPrev,1)), main=Titre)
 109     legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind))
 110     lines(1:48, dataj, lwd=2.5)
 111     lines(JourMoy, lty = 2, lwd=2)
 112     abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2)
 113     #xx=dev.off()
 114
 115     Err24 = c(Err24, erreur24)
 116     ErrPrev = c(ErrPrev, erreurPrev)
 117     ResDates = cbind(ResDates, rownames(data)[ind])
 118 }
 119
 120 rownames(ResDates) = 1:10
 121
 122 Kvois = NULL
 123 for (Col in ncol(ResDates):1)
 124 {
 125     K = 0
 126     for (I in 1:10)
 127     {
 128         for (J in 1:10)
 129         {
 130             if (ResDates[I,1] == ResDates[J,Col])
 131                 K = K +1
 132         }
 133     }
 134     Kvois = c(Kvois, K)
 135 }
 136
 137 # pdf("Erreur_Epandage_PMjour.pdf")
 138 ymin = min(na.omit(ErrPrev), Err24)
 139 ymin = min(ymin, Kvois)
 140 ymax = max(na.omit(ErrPrev), Err24)
 141 ymax = max(ymax, Kvois)
 142 plot(5:24,Err24, type = "l", lwd = 2, cex.axis=1.4, cex.lab = 1.5,
 143   ylab = "MAE",xlab = "Heures de prévision", ylim= c(ymin,ymax))
 144 lines(5:24,ErrPrev, lwd=2, lty = 2)
 145 legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lwd=2, cex=1.5)
 146 points(5:24, Kvois, cex=1.8, pch = 19)
 147 # xx = dev.off()
 148
 149 #length(D)