Le Séran déborde sur la route de Talissieu à Ceyzérieu lorsque le niveau d’eau atteint environ 2,10 m au pont de la Thuillière (Talissieu). Cependant cette hauteur n’est pas disponible sur internet. La station la plus proche se trouve à Belmont-Luthézieu.
Avec quelques observations simultanées nous pouvons construire les relations entre la hauteur à Talissieu et la hauteur ou le débit à Belmont-Luthézieu pour prédire à partir des informations sur internet l’inondation à Talissieu.
Le modèle n’est pas parfait car entre la station de Belmont-Luthézieu et Talissieu le séran reçoit en affluent : le Flon, le Groin, la rivière d’Yon, le ru de l’eau morte, le Laval et un autre petit ruisseau (à Ameyzieu).
library(ggplot2)
#library(ggrepel)
seran <-read.table(textConnection("ht\thbl\tqbl
1,6 1,20 12,300
0,75 0,80 1,300
0,8 0,78 1,260
0,75 0,70 0,650
1,05 0,87 2,720
1,15 0,87 2,750
0,8 0,73 0,766
0,65 0,67 0,271
0,62 0,63 0,166
0,61 0,62 0,116
1,35 1,20 11,100
1,97 1,48 19,500
2,32 1,76 28,900
2,1 1,4 17,000
1,79 1,24 12,300
1,78 1,21 11,400
0,97 0,87 2,500
0,47 0,58 0,017
0,97 0,82 1,860
2,1 1,67 27,000
1,41 1,05 6,000
1,7 1,2 11,100
2,1 1,43 17,900
1,34 1,291 13,700
1,56 1,4 17,000
1,98 1,55 21,600
1,11 1,01 5,770
1,43 1,04 6,570
1,41 1,03 6,290
1,09 0,88 2,900
0,80 0,84 2,190
1,06 0,82 1,530
0,67 0,7 0,560
0,58 0,68 0,439
1,06 0,89 3,090
0,87 0,84 2,040
"), dec = ",", sep = "\t", header = TRUE)
#attach(seran)
debordtal <- 2.1 # hauteur inondation mesurée à Talissieu (m)
predx <- data.frame(ht = seq(from = min(seran$ht), to = max(seran$ht), by = 0.05))resh <- lm(data = seran, hbl ~ ht)
summary(resh)##
## Call:
## lm(formula = hbl ~ ht, data = seran)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12527 -0.07339 -0.02428 0.05295 0.21364
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29191 0.03855 7.572 8.54e-09 ***
## ht 0.58616 0.02863 20.474 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08868 on 34 degrees of freedom
## Multiple R-squared: 0.925, Adjusted R-squared: 0.9228
## F-statistic: 419.2 on 1 and 34 DF, p-value: < 2.2e-16
#1 - (resh$deviance / resh$null.deviance)
shapiro.test(resh$residuals) # dê > 0.05##
## Shapiro-Wilk normality test
##
## data: resh$residuals
## W = 0.94643, p-value = 0.08075
par(mfrow = c(2,2))
plot(resh)par(mfrow = c(1,1))
# valeurs débordement
(ph <- predict(resh, se.fit = TRUE, newdata = data.frame(ht = c(debordtal)), interval = "prediction"))## $fit
## fit lwr upr
## 1 1.522839 1.333465 1.712213
##
## $se.fit
## [1] 0.02862845
##
## $df
## [1] 34
##
## $residual.scale
## [1] 0.08867805
phval <- format(ph$fit[1], digits = 3, nsmall = 2)
phmin <- format(ph$fit[2], digits = 3, nsmall = 2)
phmax <- format(ph$fit[3], digits = 3, nsmall = 2)
# intervalle de confiance sur la prédiction
ph_int <- cbind(predx, predict(resh, newdata = predx, interval = "prediction", level = 0.95))Débordement à Talissieu : 2.1 m, soit à Belmont : 1.52 m [1.33, 1.71]
# reslin <- glm(seran$qbl ~ seran$ht)
# summary(reslin)
# 1 - (reslin$deviance / reslin$null.deviance)
#
# rescube <- glm(seran$qbl ~ I(seran$ht^3))
# summary(rescube)
# 1 - (rescube$deviance / rescube$null.deviance)
respol3 <- lm(data = seran, qbl ~ poly(ht, 3)) # pour avoir les puissances significatives
#respol3 <- glm(qbl ~ poly(ht, 3, raw = TRUE)) # pour avoir les coefficients "réels" de l'équation
summary(respol3)##
## Call:
## lm(formula = qbl ~ poly(ht, 3), data = seran)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7376 -1.9311 -0.0863 0.8019 6.5129
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.57097 0.44016 17.200 < 2e-16 ***
## poly(ht, 3)1 44.88388 2.64097 16.995 < 2e-16 ***
## poly(ht, 3)2 9.25583 2.64097 3.505 0.00138 **
## poly(ht, 3)3 0.07922 2.64097 0.030 0.97625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.641 on 32 degrees of freedom
## Multiple R-squared: 0.9039, Adjusted R-squared: 0.8949
## F-statistic: 100.4 on 3 and 32 DF, p-value: 2.316e-16
#1 - (respol3$deviance / respol3$null.deviance)
par(mfrow = c(2, 2))
plot(respol3)par(mfrow = c(1, 1))
# valeurs débordement
(pq <- predict(respol3, se.fit = TRUE, newdata = data.frame(ht = c(debordtal)), interval = "prediction"))## $fit
## fit lwr upr
## 1 21.73765 15.99313 27.48216
##
## $se.fit
## [1] 0.989277
##
## $df
## [1] 32
##
## $residual.scale
## [1] 2.640973
pqval <- format(pq$fit[1], digits = 3, nsmall = 2)
pqmin <- format(pq$fit[2], digits = 3, nsmall = 2)
pqmax <- format(pq$fit[3], digits = 3, nsmall = 2)
# intervalle de confiance sur la prédiction
pq_int <- cbind(predx, predict(respol3, newdata = predx, interval = "prediction", level = 0.95))Débordement à Talissieu : 2.1 m, soit à Belmont : 21.74 m3.s-1 [15.99, 27.48]
ggplot(seran, aes(ht, hbl)) +
annotate("rect", xmin = debordtal, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = 0.1, fill = "blue") +
geom_point() +
stat_smooth(method = "lm") +
geom_line(data = ph_int, aes(ht, lwr), linetype = "dashed") +
geom_line(data = ph_int, aes(ht, upr), linetype = "dashed") +
geom_ribbon(aes(ymax = as.numeric(phmax), ymin = as.numeric(phmin)), alpha = 0.1) +
annotate("segment", x = min(seran$ht), xend = max(seran$ht), y = as.numeric(phval), yend = as.numeric(phval), color = "darkgrey") +
annotate("text", x = min(seran$ht) + 0.09, y = as.numeric(phval) + 0.03, label = phval) +
xlab(expression(h[Talissieu] (m))) +
ylab(expression(h[Belmont-Luthézieu] (m))) +
ggtitle("Hauteur Séran")ggplot(seran, aes(ht, qbl)) +
geom_rect(data = seran[1, ], aes(xmin = debordtal, xmax = Inf, ymin = -Inf, ymax = Inf), alpha = 0.1, fill = "blue") +
geom_point() +
#geom_text(aes(ht,qbl), label = rownames(seran), vjust=-.25, hjust=-.25) +
#geom_text_repel(aes(ht,qbl, label = rownames(seran))) +
stat_smooth(method = "lm", formula = "y ~ poly(x, 3)", color = "red") +
geom_line(data = pq_int, aes(ht, lwr), linetype = "dashed") +
geom_line(data = pq_int, aes(ht, upr), linetype = "dashed") +
#geom_hline(yintercept = pq$fit, color = "darkgrey") +
#stat_smooth(method = "glm", formula = "y ~ poly(x, 2)") +
geom_ribbon(aes(ymax = as.numeric(pqmax), ymin = as.numeric(pqmin)), alpha = 0.1) +
annotate("segment", x = min(seran$ht), xend = max(seran$ht), y = as.numeric(pqval), yend = as.numeric(pqval), color = "darkgrey") +
annotate("text", x = min(seran$ht) + 0.09, y = as.numeric(pqval) + 0.8, label = pqval) +
xlab(expression(h[Talissieu] (m))) +
ylab(expression(paste(Q[Belmont-Luthézieu], "(", m^3%.%s^-1, ")"))) +
ggtitle("Débit Séran")(resh2 <- lm(data = seran, ht ~ hbl))##
## Call:
## lm(formula = ht ~ hbl, data = seran)
##
## Coefficients:
## (Intercept) hbl
## -0.3673 1.5780
save(resh2, phmin, phval, pqmin, pqval, file = "seran_glm.Rdata")