Tag Archives: meteo

Computer

Carte météo dynamique

1 Comment

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…

Computer

Faire ses cartes météo

Leave a comment

Allez, un peu de scraping aujourd’hui (promis, il y aura d’autres billets d’ici la fin de la semaine sur le sujet). Je vais m’inspirer des codes de Romain de la formation Data Science pour l’Actuariat. Le but est de récupérer des données de températures, de précipitations, de vitesse de vent, journalières, en France. On commence par charger quelques librairies,

library(plyr)
library(stringr)
library(OpenStreetMap)
library(leaflet)
library(shiny)
library(rsconnect)
library(mapview)
library(png)
library(magick)
library(yaml)

On va récupérer les deux dernières années, mais comme toujours, le code peut facilement s’adapter.

annee = c(2016:2017)
mois = c('01','02','03','04','05','06','07','08','09','10','11','12')
aaaamm = sort(sub(pattern=' ',replacement = '', x=outer(annee,mois,paste)))
myCols = c(NA,NA,"NULL","NULL","NULL",NA,NA,"factor","NULL","NULL","NULL",NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL",NA,NA,NA,NA,NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL")
for (mois in aaaamm) {
  fichier = paste0("synop_",mois,".csv.gz")
  if (!file.exists(fichier)) {
    page = paste0('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/Archive/synop.',mois,'.csv.gz')
    download.file(page, fichier, quiet = TRUE, cacheOK = FALSE)
  }
}  
for (mois in aaaamm) {
  fichier = paste0("synop_",mois,".csv.gz")
  if (mois=='201601') {
    data = read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols)
  } else {
    data = rbind(data,read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols))
  }
}

On supprime une des colonnes qui ne sert a rien

data[,"X"] = NULL

et on bricole un peu (sinon on a des soucis avec minuit, qui est décalé de 24 heures)

data$temp = as.numeric(substr(data$date,1,14))
decalage_minuit = function(dateheure){
  r = dateheure
  if (dateheure-trunc(dateheure/1000000)*1000000==0){
    r = dateheure + 10000
  }
  return(r)
}
data$temp = sapply(data$temp,decalage_minuit)
data$jourheure_obs = strptime(data$temp, format = '%Y%m%d%H%M%S', 'UTC')
data$date = as.POSIXct((strptime(data$jourheure_obs,format = "%Y-%m-%d")))
data$heure = as.numeric(substr(data$jourheure_obs, start=12, stop=13))
data[,"temp"] = NULL
data$t = as.numeric(as.character(data$t))-273.15

oui, on doit convertir en degrés Celsius. On peut aussi travailler sur les vitesses de vent

data$dd = as.numeric(as.character(data$dd))
data$ff = as.numeric(as.character(data$ff))/1000*3600
data$ff = round(data$ff/5)*5

que l’on va convertir en km/h, et enfin les précipitations

data$rr1 = pmax(0,as.numeric(as.character(data$rr1)))
data$rr3 = pmax(0,as.numeric(as.character(data$rr3)))
data$rr6 = pmax(0,as.numeric(as.character(data$rr6)))
data$rr12 = pmax(0,as.numeric(as.character(data$rr12)))
data$rr24 = pmax(0,as.numeric(as.character(data$rr24)))

On peut d’ailleurs commencer par une simple visualisation de séries de températures

data_temperature = subset(data, !is.na(data$t))
temp_min = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=min)
names(temp_min)[2] = "temp_min"
temp_max = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=max)
names(temp_max)[2] = "temp_max"
temp_moy = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=mean)
names(temp_moy)[2] = "temp_moy"
temp_stat = cbind(temp_min,temp_max[],temp_moy)
rm(temp_min,temp_max,temp_moy)
temp_stat[,3] = NULL
temp_stat[,4] = NULL

avec le maximum observé, par jour, sur toutes les stations, le minimum, et la moyenne des stations

plot(x=temp_stat$date,y=temp_stat$temp_moy,'l',ylim=c(-15,40),xlab='Date',ylab='Température (°C)')
title('Températures françaises métropolitaines de 01/2016 à 12/2017')
points(x=temp_stat$date,y=temp_stat$temp_min,col="blue",'l')
points(x=temp_stat$date,y=temp_stat$temp_max,col="red",'l')

Pour faire des jolies cartes, il nous faut des informations sur les stations météo.

stations = read.table('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/postesSynop.csv',header=TRUE, sep=";",quote='"',fileEncoding="UTF-8")
stations = rename(stations, replace=c("ID"="numer_sta"))
stations = rename(stations, replace=c("Nom"="nom"))
stations = rename(stations, replace=c("Latitude"="latitude"))
stations = rename(stations, replace=c("Longitude"="longitude"))
stations = rename(stations, replace=c("Altitude"="altitude"))
data = merge(x=data,y=stations)

Pour la carte, on va se restreindre aux données en France métropolitaine

data = subset(data, data$latitude >= 39)
data = subset(data, data$latitude <= 53)
data = subset(data, data$longitude >= -5.3)
data = subset(data, data$longitude <= 9.8)
coordonnees_stations = unique(data[, c("nom","latitude","longitude")])

Commencons par visualiser le maximum observee le 19 juillet 2016

date_sel = "2016-07-19"
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))
temp = aggregate(x=donnees_carte_0$t,by=list(numer_sta=donnees_carte_0$numer_sta),FUN="max")
donnees_carte_0 = merge(x=donnees_carte_0,y=temp)
donnees_carte_0 = rename(donnees_carte_0, replace=c("x"="TEMP"))
donnees_carte_0$TEMP = round(donnees_carte_0$TEMP)
donnees_carte = unique(donnees_carte_0[,c('numer_sta','nom','longitude','latitude','TEMP')])
get_france = get_map(c(lon=2.25,lat=46), zoom=5, col='bw')
colfunc = colorRampPalette(c("darkblue","royalblue","cyan","lightblue","orange","red","darkred"))
liste_couleurs = data.frame(TEMP=as.numeric(seq(-49,50,1)),couleur=as.character(colfunc(100)))
donnees_carte = merge(x=donnees_carte,y=liste_couleurs)
france_tmax = ggmap(get_france) + 
  ggtitle(paste0('Températures maximales du ', date_sel_OK)) +
  scale_x_continuous(limits = c(-5, 10), expand = c(0, 0)) +
  scale_y_continuous(limits = c(41.3, 51.1), expand = c(0, 0))
france_tmax = france_tmax + geom_label(data=donnees_carte,aes(x=longitude,y=latitude,label=TEMP),col=as.character(donnees_carte$couleur))
france_tmax

On peut tenter aussi une visualisation des précipitations. On se place le 13 octobre 2016 (oui, il y a eu des inondations dans la region de Montpelier ce jour la)

date_sel = "2016-10-13"
colfunc = colorRampPalette(c('white','darkblue'))
data$precip_en_mm = data$rr3
data_precipitation = subset(data, !is.na(data$precip_en_mm))
precipitations = aggregate(x=data_precipitation$precip_en_mm,by=list(date=data_precipitation$date,station=data_precipitation$nom),FUN=sum)
names(precipitations)[3] = "pp"
precipitations$pp = pmax(precipitations$pp,0)
donnees_carte = subset(precipitations, precipitations$date==date_sel)
qt_prec = quantile(round(donnees_carte$pp),seq(0, 1, 0.05))
donnees_carte$Couleur_prec = colfunc(21)[(findInterval(round(donnees_carte$pp), qt_prec, all.inside=TRUE))]
donnees_carte = rename(donnees_carte, replace=c("station"="nom"))
donnees_carte = merge(x=donnees_carte,y=coordonnees_stations)

On peut tenter un leaflet ici (ça sera plus joli)

france = leaflet() %>% addTiles() %>% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)
france = france %>% addCircles(lng=donnees_carte$longitude, lat=donnees_carte$latitude, color=donnees_carte$Couleur_prec, opacity = 1, fillColor=donnees_carte$Couleur_prec, fillOpacity = 1, radius=donnees_carte$pp*500)
france

Comme souvent, j’ai des soucis pour integrer du leaflet en wordpress, alors on se contentera dans le billet d’une copie d’ecran (mais un copier/coller du code suffit).

Et pour finir, la meme chose avec les vitesses de vent. On se place le 6 mars 2017 (la encore, la date n’est pas choisie au hasard, c’était la tempete Zeus)

date_sel = "2017-03-06"
date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))
heure_sel = 18
donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))
donnees_carte_0 = subset(donnees_carte_0, donnees_carte_0$heure==heure_sel)
donnees_carte = donnees_carte_0

Pour les precipitations, des cercles plus ou moins grands suffisaient. La, on va mettre des fleches

logo = image_trim(image_read("https://image.freepik.com/icones-gratuites/fleche-noire-vers-le-haut_318-30934.jpg"))
if (file.exists('flecheblanche90.png')==FALSE){
  for (angle in seq(10,360, by=10)) {
    temp = image_rotate(logo,angle)
    image_write(temp,paste0('flecheblanche',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='green')
    image_write(temp,paste0('flecheverte',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='yellow')
    image_write(temp,paste0('flechejaune',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='orange')
    image_write(temp,paste0('flecheorange',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='red')
    image_write(temp,paste0('flecherouge',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='black')
    image_write(temp,paste0('flechenoire',angle,'.png'))
  }
}
arrowIcons = icons(
  iconUrl = ifelse(donnees_carte$dd == 0,
            ifelse(donnees_carte$ff<39,
                    'flecheverte360.png',
            ifelse(donnees_carte$ff<79,
                    'flechejaune360.png',
            ifelse(donnees_carte$ff<119,
                    'flecheorange360.png',
            ifelse(donnees_carte$ff<159,
                    'flecherouge360.png',
            ifelse(donnees_carte$ff>159,
                    'flechenoire360.png',
                    'flecheblanche360.png'))))),
            ifelse(donnees_carte$ff<39,
                   paste0('flecheverte',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff<79,
                    paste0('flechejaune',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff<119,
                    paste0('flecheorange',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff<159,
                    paste0('flecherouge',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff>159,
                    paste0('flechenoire',donnees_carte$dd,'.png'),
           paste0('flecheblanche',donnees_carte$dd,'.png'))))))),
  iconWidth = 30, iconHeight = 30)

Cette fois, on est bon ! On peut faire la carte

france = leaflet(options = leafletOptions()) %>% addTiles() %>% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)
france = addMarkers(map = france,lng=donnees_carte$longitude,lat=donnees_carte$latitude,icon=arrowIcons,popup=as.character(donnees_carte$ff))
france