Carte météo dynamique

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…


Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.