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

OpenEdition suggests that you cite this post as follows:
Arthur Charpentier (September 19, 2018). Faire ses cartes météo. Freakonometrics. Retrieved September 19, 2024 from https://doi.org/10.58079/ovbn

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.