Pour aller un peu plus loin par rapport au précédent billet, on peut faire une carte dynamique, pour visualiser une tempête, sur une journée. Je vais partir du fait qu’on a encore les objets du précédant billet en mémoire.
On commence par garder uniquement les données relatives a la journée qui nous intéresse,
date_sel = "2016-10-13" date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4)) donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel)) horaires = sort(unique(donnees_carte_0$heure)) |
(mais on pourrait prendre une semaine). On recupere ensuite le contour de la France,
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds") FR0=readRDS("FRA_adm0.rds") P1=FR0@polygons[[1]]@Polygons[[355]]@coords P2=FR0@polygons[[1]]@Polygons[[27]]@coords |
car seule la France métropolitaine nous intéressera (incluant la Corse, soit 2 gros polygones) dans l’esprit du précédant billet sur le zonier. On va ainsi definir un maillage sur lesquel on va lisser la pluviometrie,
grille = expand.grid(seq(min(donnees_carte_0$longitude),max(donnees_carte_0$longitude),length=101),seq(min(donnees_carte_0$latitude),max(donnees_carte_0$latitude),length=101)) paslong=(max(donnees_carte_0$longitude)-min(donnees_carte_0$longitude))/100 paslat=(max(donnees_carte_0$latitude)-min(donnees_carte_0$latitude))/100 f=function(i){ (point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P1[,1],P1[,2])>0)+(point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P2[,1],P2[,2])>0) } indic=unlist(lapply(1:nrow(grille),f)) grille=grille[which(indic==1),] |
Pour lisser, on utilise les k-plus proches voisins
knn=function(i,k=10){ d=distHaversine(grille[i,1:2],donnees_carte[, c("longitude","latitude")], r=6378.137) r=rank(d) ind=which(r<=k) weighted.mean(donnees_carte[ind,"rr3"],(1/d[ind])/sum(1/d[ind]))} |
On a ensuite la fonction suivante, pour faire une carte
carto_prec<-function(){ grille2<-grille grille2$rr=Vectorize(knn)(i=1:nrow(grille2)) bk=seq(0,50,length=21) grille2$cuty=cut(grille2$rr,breaks=bk,labels=1:20) cols = rev(carto.pal(pal1 = "blue.pal", n1=20, pal2 = "white.pal", n2=1)) plot(FR0,border=NA) polygon(P1) polygon(P2) points(grille2[,1]+paslong/2,grille2[,2]+paslat/2,col=cols[grille2$cuty],pch=19) points(donnees_carte$longitude,donnees_carte$latitude, col="black",pch=19,cex=.5) title(main = paste0(date_sel_OK," à ",heure_sel,"H"), sub = "Précipitations des 3 dernières heures", cex.main = 1.5, font.main= 4, cex.sub = 1, font.sub = 3) legend(8.2, 50, legend=seq(0,50,length=5), title='en mm', fill=cols[seq(1,20,length=5)], cex=0.8) } |
On va ensuite utiliser 8 fichiers – car on regarde toutes les 3 heures
for (hh in horaires) { heure_sel<-hh donnees_carte <- subset(donnees_carte_0, donnees_carte_0$heure==heure_sel) png(paste0("CartePrec",date_sel,"_",hh,".png")) carto_prec() dev.off() } |
(mais ça peut se raffiner, avec un lissage spatio-temporel). Et on va ensuite créer notre animation en concaténant les images
frames <- image_morph( c(image_scale(image_read("CartePrec2016-10-13_1.png")), image_scale(image_read("CartePrec2016-10-13_3.png")), image_scale(image_read("CartePrec2016-10-13_6.png")), image_scale(image_read("CartePrec2016-10-13_9.png")), image_scale(image_read("CartePrec2016-10-13_12.png")), image_scale(image_read("CartePrec2016-10-13_15.png")), image_scale(image_read("CartePrec2016-10-13_18.png")), image_scale(image_read("CartePrec2016-10-13_21.png"))), frames = 16) frames.anim <- image_animate(frames) image_write(frames.anim, paste0("CartePrecDay",date_sel,".gif")) |

Il faut avouer que ça a de gueule… On voit clairement la tempête arriver…
Hier matin, au moment de quitter la maison, j’ai jeté un œil à la météo (car à Montréal le temps est très changeant, pas seulement d’un jour sur l’autre, comme je l’avis mentionné en arrivant en septembre, 
la variable indicatrice indiquant s’il pleut, ou pas (1 s’il pleut, et 0 sinon) à l’instant
. Supposons qu’il existe un processus latent sous-jacent,
tel que

qui peut encore s’écrire
ou, ce qui sera pleut être plus utile ensuite
Avec notre modèle latent, cette probabilité peut s’écrire
c’est à dire en notant
la fonction de répartition (jointe) du vecteur aléatoire
cette probabilité s’écrit
Or Wassily Hoeffding a montré qu’il était possible d’encadrer une fonction de répartition, avec les bornes suivantes: pour tout 


où 

(qui est associée à la probabilité d’avoir de la pluie pendant une heure en moyenne) et de la corrélation
(qui est associée à la dynamique du processus). Pour simplifier, si on prend 


















