Bonjour,

Pour réaliser ce graphique http://meteo.besse83.free.fr/imfix/trendtest-Moy%205%20s%E9ries-1998-30%202034.png une personne utilise un srcipt à voir ci-dessous suivant son message.

Le script qui permet de calculer l'IP à 95% de ce graphe est :

t <- read.table("C:/r/moy5.csv", sep=";", dec=",", header=T)

# note: this does not yet handle NAs in the trend well

# be careful with the bounderies - 1979 for sats, pre-1900 for ground

# need boundary tests

y1 <- 1998 # end of trend

t0 <- 30 # years of trend

y0 <- y1 - t0 # start of trend

t1 <- 3 * t0 # number of years to display

y2 <- y0 + t1 # max of x-axis

ymin <- -1.0 # min of y-axis

ymax <- 1.5 # max of y-axis

# column idx

gis <- 3

cru <- 4

noa <- 5

uah <- 6

rss <- 7

moy5 <- 8

idx <- gis

src = c("GISTEMP", "HADCRU", "NOAA", "UAH", "RSS", "Moy 5 séries")

# annualize monthly data

t$year <- as.integer(t[,1])

# get the record of interest, monthly to annual

Tann <- aggregate(t[,idx],list(t$year), FUN=mean, na.rm=T)

target <- Tann[Tann[,1]> y0 & Tann[,1]<= y1,]

# a series for x-axis

tyear <- y0:y2

# some set setup

x <- target[,1]

y <- target[,2]

new <- data.frame(x = seq(y0+1,y1))

# get trend; hi and lo for 95% ci

trend_lm <- lm(y ~ x)

trend_lo = cbind(x,predict(lm(y ~ x), interval="prediction")[,2])

trend_hi = cbind(x,predict(lm(y ~ x), interval="prediction")[,3])

lm_lo <- lm(trend_lo[,2] ~ trend_lo[,1])

lm_hi <- lm(trend_hi[,2] ~ trend_hi[,1])

# extend trend +/- ci into future years

trend_lm_lo <- tyear*lm_lo$coef[2] + lm_lo$coef[1]

trend_lm_hi <- tyear*lm_hi$coef[2] + lm_hi$coef[1]

trend_lm_mid <- tyear*trend_lm$coef[2] + trend_lm$coef[1]

# set y axis

ymax = max(trend_lm_hi) + 0.5

ymin = min(trend_lm_lo) - 0.5

# file name to write plot

png(paste("trendtest-",src[idx-2],"-",y1,"-",t0,".png",sep=""),height=480, width=540, bg="white")

# plot

dstamp=Sys.Date()

plot(tyear,rep(0,length(tyear)), ylim=c(ymin,ymax), col="white", xlab="Années", ylab="Anomalie (°C)") # hidden - it just sets the graph axis

title(main=paste("Le réchauffement s'est-il arrêté en",y1,"?\n","Test tendance/",t0,"ans",src[idx-2]),

sub=paste("http://rhinohide.wordpress.com+ modif ChristianP ", dstamp), cex.lab=0.75)

# warm, cold, no trend boxes

lines(c(y1,y1),c(ymin,ymax), type="l", col="gray") # beginning year of trend test

rect(y1,trend_lm_hi[t0+1],y2,ymax, col="lightpink")

rect(y1,trend_lm_lo[t0+1],y2,ymin, col="lightskyblue")

vtx <- min(min(tyear[trend_lm_lo > trend_lm_hi[t0+1]]), y2)

polygon(rbind(c(y1,trend_lm_lo[t0+1]),c(vtx,trend_lm_lo[vtx-y0+1]),c(y2,trend_lm_hi[t0+1]),c(y2,trend_lm_lo[t0+1])), col="light gray", lty="solid")

# label trend boxes

text(y2,(trend_lm_lo[t0+1]+trend_lm_hi[t0+1])/2, "Stabilité", pos=2, cex=0.8)

text(y1+1,(ymin+trend_lm_lo[t0+1])/2, "Refroidissement", pos=4, cex=0.8)

text(y1+1,(ymax+trend_lm_hi[t0+1])/2, "Réchauffement", pos=4, cex=0.8)

text(y0+1,(ymax+trend_lm_hi[t0+1])/3, paste("Tendance/",t0," ans \n95% intervalle de prédiction"), pos=4, cex=0.8)

# trend slope data

a1 <- format(trend_lm$coeff[2]*10,digits=2)

a2 <- format(trend_lm_hi[t0+1] - trend_lm_mid[t0+1],digits=2)

a3 <- format(-trend_lm_lo[t0+1] + trend_lm_mid[t0+1],digits=2)

text(y0+1,ymin+0.1, paste(a1,"°C par décennie\n IP +",a2,"/ -",a3), pos=4, cex=0.8)

# points in the trend

points(target[,1],target[,2], pch=16) # data points for annual global T

# trend lines and +/- prédiction interval

abline(trend_lm,col="black", lwd=2) # trendline for annual global T

lines(trend_lm_hi ~ tyear,col="red", lty=2, lwd=2) # 2 SD upper limit of monthly global T

lines(trend_lm_lo ~ tyear,col="blue", lty=2, lwd=2) # 2 SD lower limit of monthly global T

# add the trend test points

t_test <- Tann[Tann[,1]> y1,]

points(t_test[,1],t_test[,2], col="green", pch=16) # data points beyond the trend setters

dev.off()

Le fichier csv pour le script dessus est : http://meteo.besse83.free.fr/imfix/moy5.csv

Suivant ce script d'après-vous dans un fichier Excel avec les données et années de l'anomalie de la température comme ici http://cjoint.com/?0AypKsDGKy0 quels formules doit on y mettre pour pouvoir réaliser le graphique du 1er lien ci-dessus qui est réalisé avec ce script ???

Merci

Williams

Modifié par williamsss le 24/01/2015 15:34

Bonjour

Déjà ici il y a une adresse qui ne fonctionnera peut être pas

sub=paste("http://rhinohide.wordpress.com+ modif ChristianP ", dstamp), cex.lab=0.75)

http://meteo.besse83.free.fr/

w36xb2w a écrit :

Bonjour

Déjà ici il y a une adresse qui ne fonctionnera peut être pas

sub=paste("http://rhinohide.wordpress.com+ modif ChristianP ", dstamp), cex.lab=0.75)

http://meteo.besse83.free.fr/

Bonjour W36xb2w,

L'adresse http://rhinohide.wordpress.com fonctionne bien mais ci-dessus c'est le fait qu'il y a le "+" qui s'est incorporer dans l'adresse du lien comme tout est rattaché dans ce que la personne m'a donné concernant le code du script.

merci

Williams