investigations on Bruno vs. talweg: a shift overved, and neighbors differ (a bit...

[talweg.git] / Etude_En_Cours_R_E_bis.R

diff --git a/Etude_En_Cours_R_E_bis.R b/Etude_En_Cours_R_E_bis.R
index aa0d490..1fe1387 100644 (file)
--- a/Etude_En_Cours_R_E_bis.R
+++ b/Etude_En_Cours_R_E_bis.R
@@ -1,5 +1,7 @@
+predi <- function(ij)
+{
 #setwd("/Users/bp/Desktop/CONTRATS_AirNormand/2016/RapportFinalBruno")
-rm(list=ls())
+#rm(list=ls())
 
 # 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)
@@ -48,7 +50,7 @@ ResDates = NULL
 
 nbvois=10
 j=1 # numéro de semaine
-ij=6 # numéro du jour (0 = lundi)
+#ij=6 # numéro du jour (0 = lundi)
 
 Err24 = NULL
 ErrPrev = NULL
@@ -79,6 +81,18 @@ Kvois = NULL
        large = 1
        bornes = mean(dataj[25:48]) + c(-large,large)
        indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
+       if (sum(indcond) < 10)
+       {
+               large = 2
+               bornes = mean(dataj[25:48]) + c(-large,large)
+               indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
+       }
+       while (sum(indcond) < 10)
+       {
+               large = large + 3
+               bornes = mean(dataj[25:48]) + c(-large,large)
+               indcond = varexp[,"PMjour"]>=bornes[1] & varexp[,"PMjour"]<=bornes[2]
+       }
        data = data[indcond,] #pollution du 2eme jour == pollution du jour courant +/- 1
        varexp = varexp[indcond,]
 
@@ -92,6 +106,13 @@ Kvois = NULL
 #      w = 1/(D[ind]^2)
 #      w = w/sum(w)
 #      W = w %o% rep(1,48)
+
+print("Voisins + ij")
+print(sort(D)[1:nbvois])
+print(dateJPrev)
+print(rownames(data)[ind])
+print(data[ind, 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="")
@@ -102,6 +123,8 @@ Kvois = NULL
        } else {
                erreurPrev = mean(abs(dataj[(H+1):48] - JourMoy[(H+1):48]))
        }
+
+       list(serie=dataj, prev=JourMoy, err=erreurPrev, line=(nl+ij), neighbs=ind, dates=data[ind,1])
 #      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), 
@@ -149,3 +172,4 @@ Kvois = NULL
 ## xx = dev.off()
 #
 ##length(D)
+}