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
   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)