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…