From: Benjamin Auder Date: Tue, 25 Apr 2017 21:15:41 +0000 (+0200) Subject: commented a bit Bruno's file X-Git-Url: https://git.auder.net/variants/current/doc/css/assets/%3C?a=commitdiff_plain;h=9dbd32d9f73bbc91c65c2c6a05095f084c6fdb16;p=talweg.git commented a bit Bruno's file --- diff --git a/Etude_En_Cours_R_E_bis.R b/Etude_En_Cours_R_E_bis.R index 69b1da0..40c7954 100644 --- a/Etude_En_Cours_R_E_bis.R +++ b/Etude_En_Cours_R_E_bis.R @@ -1,29 +1,37 @@ -setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno") +#setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno") rm(list=ls()) -# Lecture des données -pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=";",dec=".",header=T) #,row.names=1) -n = dim(pm)[1] -datedebut = "10/12/2008" +# Lecture des données: pm = dataframe 2 colonnes, date-time puis PM10 horaire +pm = read.table("DATA/mesures_horaires_hloc_pm10_a_filer.csv",sep=";",dec=".",header=T) +#n = dim(pm)[1] +#datedebut = "10/12/2008" -# Chargement des données météo et indicateurs +# Chargement des données météo et indicateurs: VarExp matrice des données météo, +# première colonne = date, première rangée = second jour VarExp <- read.table("DATA/meteo_extra_jourMois.csv",sep=";",dec=".",header=T) VarExp <- VarExp[-1,] -# Lecture des dates -dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T) +# Lecture des dates: ??? sans doute liste des jours de 10/12/2008 à 06/2016... +dates = read.table("DATA/Dates_jours.csv",header=F,as.is=T) dates = dates[,1] -# Contruction des matrices de données +# Contruction des matrices de données: pm.h = matrice des séries en ligne pm.h <- matrix(pm[,2],ncol=24,byrow=TRUE) - Nlignes = nrow(pm.h) + +# Data = matrice des couples de jours (séries de 48 PM10), +# première ligne = jour 1/2, dernière = jour N-1/N Data = cbind(pm.h[1:(Nlignes -1), ], pm.h[2:Nlignes, ]) + +# dates2 = dates du 2eme jour au dernier dates2 = dates[2:Nlignes] + rownames(Data) = dates2 # df contient l'ensemble des données. #df <- cbind(Data,varexp[,-1]) df <- Data + +# Complétion des variables exogènes (du 2eme jour) par les PM10 moyen ce même jour PMjour <- apply(df[,25:48],1,mean,na.rm=T) dfexp <- cbind(VarExp,PMjour) @@ -46,79 +54,84 @@ Err24 = NULL ErrPrev = NULL Kvois = NULL -for (Hc in 5:24) { - -H=24+Hc -fen=H -#fen=12 ou 24 ou H -L = (H-fen+1):H - -# Premier conditionnement : mois -indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10 -data = df[indcond,] -varexp = dfexp[indcond,] - -nl = (1:nrow(data))[rownames(data)==Dates[j]] -dateJPrev = rownames(data)[nl+ij] -dataj = as.numeric(data[nl+ij, 1:48]) -data = data[-(nl + ij), ] -varexp = varexp[-(nl + ij), ] -indNA = attr(na.omit(data[, 1:48]),"na.action") -data = data[-indNA,] -varexp = varexp[-indNA,] - -# Conditionnement : les jours avec PMjour +/- large -large = 1 -bornes = mean(dataj[25:48])+c(-large,large) -indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2] -data = data[indcond,] -varexp = varexp[indcond,] - -D = rep(0,nrow(data)) -for (k in 1:nrow(data)) { - #D[k] = sqrt(sum((1:H)*(dataj[L] -data[k,L])^2)) - D[k] = sqrt(sum((dataj[L] -data[k,L])^2)) -} -ind = order(D)[1:nbvois] -w = 1/(D[ind]^2) -w = w/sum(w) -W = w%o%rep(1,48) -JourMoy = apply(data[ind, 1:48], 2, mean) -#JourMoy = apply(W*data[ind, 1:48], 2, sum) -NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="") -Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="") -#erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2)) -if(Hc==24){erreurPrev=NA}else{erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))} -erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48])) -png(NomFile) -matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind), - cex.axis=1.4, cex.main = 1.7, cex.lab=1.5, - xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1)," / ",round(erreurPrev,1)), main=Titre) -legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind)) -lines(1:48, dataj, lwd=2.5) -lines(JourMoy, lty = 2, lwd=2) -abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2) -xx=dev.off() - - - - Err24 = c(Err24, erreur24) - ErrPrev = c(ErrPrev, erreurPrev) -ResDates = cbind(ResDates, rownames(data)[ind]) +for (Hc in 5:24) +{ + H=24+Hc #Hc: dernier PM10 connu avant prédiction des 24-Hc restants, le jour 2 + L = 1:H + + # Premier conditionnement : mois + indcond <- dfexp[,"Mois"] == 2 | dfexp[,"Mois"] == 3 | dfexp[,"Mois"] == 4 + | dfexp[,"Mois"] == 9 | dfexp[,"Mois"] == 10 + data = df[indcond,] #restriction aux couples dont le mois du 2eme jour est 2,3,4,9 ou 10 + varexp = dfexp[indcond,] #de même sur les exogènes (+ PM10 moyens) + + nl = (1:nrow(data))[rownames(data)==Dates[j]] #numéro de ligne où 2eme jour == 16/03/2015 + dateJPrev = rownames(data)[nl+ij] # ??? la date du 1er jour du couple de jours pile une semaine après le couple 'nl' + dataj = as.numeric(data[nl+ij, 1:48]) #extraction de cette série dans dataj + data = data[-(nl + ij), ] #suppression de cette série dans data (????) + varexp = varexp[-(nl + ij), ] #idem dans les variables exogènes (date 16/03/2015 + 6 jours + indNA = attr(na.omit(data[, 1:48]),"na.action") + data = data[-indNA,] + varexp = varexp[-indNA,] + + # Conditionnement : les jours avec PMjour +/- large + large = 1 + bornes = mean(dataj[25:48])+c(-large,large) + indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2] + data = data[indcond,] + varexp = varexp[indcond,] + + D = rep(0,nrow(data)) + for (k in 1:nrow(data)) + { + #D[k] = sqrt(sum((1:H)*(dataj[L] - data[k,L])^2)) + D[k] = sqrt(sum((dataj[L] - data[k,L])^2)) + } + ind = order(D)[1:nbvois] + w = 1/(D[ind]^2) + w = w/sum(w) + W = w %o% rep(1,48) + JourMoy = apply(data[ind, 1:48], 2, mean) + #JourMoy = apply(W*data[ind, 1:48], 2, sum) + NomFile = paste("Voisins_Epandage_PMjour_Hc_",Hc,".png",sep="") + Titre = paste("Jour à prévoir : ",dateJPrev," - ", length(ind)," voisins",sep="") + #erreur = sqrt(sum((dataj[25:48] - JourMoy[25:48])^2)) + if(Hc==24) + erreurPrev=NA + else + erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48])) + erreur24 = mean(abs(dataj[25:48] - JourMoy[25:48])) + #png(NomFile) + matplot(t(data[ind, 1:48]), type = "l", lwd=1.4, lty=1, col=1:length(ind), + cex.axis=1.4, cex.main = 1.7, cex.lab=1.5, + xlab="Heures locales", ylab=paste0("PM10 - Erreurs = ",round(erreur24,1), + " / ",round(erreurPrev,1)), main=Titre) + legend("top",rownames(data)[ind],ncol=2,lwd=1.4, lty=1, col=1:length(ind)) + lines(1:48, dataj, lwd=2.5) + lines(JourMoy, lty = 2, lwd=2) + abline(v=c(24.5,H+0.5), lty = 2, lwd=1.2) + #xx=dev.off() + + Err24 = c(Err24, erreur24) + ErrPrev = c(ErrPrev, erreurPrev) + ResDates = cbind(ResDates, rownames(data)[ind]) } - rownames(ResDates) = 1:10 Kvois = NULL -for (Col in ncol(ResDates):1) { -K = 0 -for (I in 1:10){ - for (J in 1:10){ - if (ResDates[I,1] == ResDates[J,Col]) K = K +1 - } -} -Kvois = c(Kvois, K) +for (Col in ncol(ResDates):1) +{ + K = 0 + for (I in 1:10) + { + for (J in 1:10) + { + if (ResDates[I,1] == ResDates[J,Col]) + K = K +1 + } + } + Kvois = c(Kvois, K) } # pdf("Erreur_Epandage_PMjour.pdf") @@ -133,4 +146,4 @@ legend("topright", legend=c("Erreur 24h", "Erreur prévision"), lty = c(1,2), lw points(5:24, Kvois, cex=1.8, pch = 19) # xx = dev.off() -length(D) +#length(D)