Scraper, ou pas ?

Ce matin, je mettais en ligne un billet scraper la base d’incendies de forêts expliquant comment scraper la base Promethee, en remplaçant automatique une formulaire. Dans la soirée, sur Twitter me faisait remarquer que c’était intéressant, mais peut être inutile (sur cet exemple particulier en tous cas).

Sur Firefox, il existe un “Moniteur Réseau“, qui s’ouvre en cliquant sur le clavier CtrlMaj + E ( CommandOption + E quand on est est sur un Mac). Ça ouvre un espace en bas de la fenêtre

Si maintenant je lance a la main ma requête, on voit ce qui est fait

Une methode GET apparait. En cliquant dessus, on a toute l’information nécessaire a droite

Dans l’onglet ‘Reponse’, on voit l’information sur le fichier json créé,

on a plus d’information bien entendu, si on agrandit

et surtout si on va dans l’onglet ‘En-tetes’

on a l’URL complet de la requête ! On peut ensuite copier le lien, et l’ouvrir

On a ainsi a accès a presque toute la table. Quand on regarde l’URL, seuls 20 incendies sont renvoyés, a cause du &nbrLigne=20& mais on peut tenter de le changer. propose de mettre &nbrLigne=10000& et il affirmait que ça marchait. J’avoue avoir essayé, sans succès.

Mais je retiens le tuyau en tous cas. Parfois, on peut faire simple !

Scraper la base d’incendies de forêts

Aller, deux billets pour parler de scraping, ou harvesting serait peut être plus juste ? Disons qu’on va faire un code automatique pour récupérer des informations en ligne. Plus particulièrement, on va s’intéresser au cas ou les données sont obtenues par une requête manuelle, et le site n’a pas d’API pour générer un joli fichier csv. Pierre-Marie, de la formation Data Science pour l’Actuariat a propose un exemple dont je vais m’inspirer (avec une autre methode pour scraper car j’ai eu des soucis pour installer le logiciel).

En 1973 a été lancée en France Prométhée, une base de données sur les incendies de forêts de la région méditerranéenne (couvrant 15 départements du Sud-Est). Oui, Προμηθεύς est ce titan qui a voulu jouer avec le feu… On peut consulter la base sur http://promethee.com/default/incendies,

En remplissant les champs, on peut faire une extraction suivant tel ou tel critère. Par exemple, en 2016, en Corse, il y a eu 6 incendies de nature “accidentelle”,

On le voit, le lien de la page ne change pas… il faut donc trouver une methode pour remplir manuellement les champs, et récupérer les informations. En particulier, on voit qu’on peut avoir accès aux informations relatives aux différents sinistres en cliquant sur la loupe

Idéalement, on voudrait un code R, automatique, pour faire ça ! Une première solutions, proposée par Pierre-Marie, est de passer par RSelenium. Le package R (non maintenu sur le CRAN depuis quelques semaines) est base sur Selenium, “Selenium automates browsers” comme dit la publicité. Le code pour scraper une page comme par

require("devtools")
devtools::install_github("ropensci/RSelenium")
require(RSelenium)
remDr = remoteDriver(remoteServerAddr = "localhost" 
                      , port = 4445L
                      , browserName = "firefox")

Le soucis est que ce code ne tourne pas sur mon mac. Il y a de nombreuses pages qui parlent des spécificités de mac (en particulier le soucis du port par défaut 4445, qu’il faut changer). J’ai tout essayé, sans succès.

Une autre piste est de passer par un autre package

require("devtools")
devtools::install_github("ropensci/wdman")

wdman, pour “Webdriver Manager“, propose aussi de le faire, en s’appuyant sur diverses technologies, selenium, chromedriver, phantomJS binary ou encore geckodriver. Moyennant l’installation de chrome, la seconde option est intéressante

library(wdman)
cDrv = chrome()
eCaps = list(chromeOptions = list(
  args = c('--headless', '--disable-gpu', '--window-size=1280,800')
))
remDr = remoteDriver(browserName = "chrome", port = 4567, 
                     extraCapabilities = eCaps)
remDr$open()

qui marche ! On peut lui demander d’aller sur la base promethee, et de faire une copie d’écran

remDr$navigate("http://www.promethee.com/default/incendies")
remDr$screenshot(display = TRUE)

On va ensuite remplir les champs, pour notre requête, sur la page. Par exemple, pour la date de début, on peut remonter au début des années 80

webElem = remDr$findElement(using = 'name', 'dtAlerteDeb')
webElem$clearElement()
webElem$sendKeysToElement(list('01/01/1981'))

et terminer fin 2017

webElem = remDr$findElement(using = 'name', 'dtAlerteFin')
webElem$clearElement()
webElem$sendKeysToElement(list('31/12/2017'))

On se focalise seulement sur la Corse

webElem = remDr$findElement(using = 'name', 'codeg')
webElem$sendKeysToElement(list('CORSE'))

et on veut 50 observations par page

webElem = remDr$findElement(using = 'name', 'nbrLigne')
webElem$sendKeysToElement(list('50'))

Ah oui, et on va temporiser entre chaque requête

webElem = remDr$findElement(using = 'name', 'btnSubmit')
webElem$submitElement()
Sys.sleep(15)
remDr$screenshot(TRUE)

(c’est une règle assez classique quand on scrap: on évite d’aller trop vite, pour pas se faire repérer…). On va ensuite commencer a regarder ce que l’on demande

text_base = read_html(remDr$getPageSource()[[1]])
list_tr = html_nodes(text_base, "tr")
a = html_nodes(list_tr[length(list_tr)], "a")
page_max = as.numeric(html_text(a[10]))

On a environ 620 pages ! avec au moins 10 secondes pour changer de pages, et en rajoutant 15 secondes (on en parlait avant), ca fait quelques heures pour scraper.

Que veut-on récupérer en grattant la page ? On peut récupérer les coordonnées de l’incendie (latitude/longitude), les caractéristiques de l’incendie (dates, lieu, surface), et retourner le tout dans une table

scrapping = function(feux, page) {
  list_a = html_nodes(read_html(feux), "a")
  if (length(list_a) >= 2) {
    att = html_attrs(list_a[2])
    longitude = as.numeric(str_extract(att, "[0-9]\\.([0-9])*"))
    latitude = as.numeric(str_extract(att, "([0-9]){2}\\..([0-9])*"))
  } else {
    longitude = 0
    latitude = 0
  }
  list_span = html_nodes(read_html(feux), "span")
  date = html_text(list_span[2])
  commune = html_text(list_span[4])
  surface = as.numeric(html_text(list_span[5]))  
  list_div = html_nodes(read_html(feux), "div")
  if (length(list_a) >= 2) {
    infos1 = html_text(list_div[3])
    infos2 = html_text(list_div[5])
  } else {
    infos1 = html_text(list_div[2])
    infos2 = html_text(list_div[4])
  }
  alerte = gsub("First alert : |First intervention", "", str_extract(infos1, "First alert : (.)*First intervention"))
  deb_inter = gsub("First intervention : |End of intervention", "", str_extract(infos1, "First intervention : (.)*End of intervention"))
  fin_inter = gsub("End of intervention : |Origins of the alert", "", str_extract(infos1, "End of intervention : (.)*Origins of the alert"))
  code_INSEE = gsub("INSEE code : |Locality", "", str_extract(infos2, "INSEE code : (.)*Locality"))
  lieudit = str_trim(gsub("Locality : ", "", str_extract(infos2, "Locality : (.)*")))
  liste_feux = data.frame(commune, date, code_INSEE, longitude, latitude, lieudit, surface, alerte, deb_inter, fin_inter, page)
  return(liste_feux)
}

Reste a faire tourner sur les 620 pages. Pour chaque ligne, on récupère les informations !

df = data.frame()
for (page in 1:page_max) {
  text_base = read_html(remDr$getPageSource()[[1]])
  list_tr = html_nodes(text_base, "tr")
  if (length(html_nodes(read_html(as.character(list_tr[14])), "a")) == 8) {
    deb = 15
  } else {
    deb = 14
  }
  for (i in deb:(deb+49)) {
    feux = as.character(list_tr[i])
    df = rbind(df, scrapping(feux, page))
  }
  if (page == 1) {
    webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))]')
  } else {
    webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))][3]')
  }
  webElem$clickElement()
  heure = (floor((page_max-page)*15/60) - floor((page_max-page)*15/60) %% 60) / 60
  minute = floor((page_max-page)*15/60) %% 60
  pb = progress_bar$new(total = 100)
  for (i in 1:100) {
    if (i==1) {cat("page", page, "sur", page_max," -  (reste", heure, "h", minute, "min)",  "\n")}
    pb$tick()
    Sys.sleep(15 / 100)
  }
  remDr$screenshot(TRUE)
}
feux = df

Ça y est, on a une base, avec presque 31 000 lignes ! on va ensuite ranger un peu, et clarifier

colnames(feux)[colnames(feux)=="longitude"] = "x"
colnames(feux)[colnames(feux)=="latitude"] = "y"
feux[, date := as.Date(date, "%d/%m/%Y")]
feux[, date_alerte := as.Date(substr(alerte, 1, 10), "%d/%m/%Y")]
feux$heure_alerte[nchar(feux$alerte) > 10] = substr(feux$alerte[nchar(feux$alerte) > 10], 12, 16)
feux$heure_alerte[!(nchar(feux$alerte) > 10)] = ""
feux[, date_deb_inter := substr(deb_inter, 1, 10)]
feux[, date_deb_inter := as.Date(date_deb_inter, "%d/%m/%Y")]
feux$heure_deb_inter[nchar(feux$deb_inter) > 10] = substr(feux$deb_inter[nchar(feux$deb_inter) > 10], 12, 16)
feux$heure_deb_inter[!(nchar(feux$deb_inter) > 10)] = ""
feux$heure_deb_inter[feux$heure_deb_inter == "00:00"] = ""
feux[, date_fin_inter := substr(fin_inter, 1, 10)]
feux[, date_fin_inter := as.Date(date_fin_inter, "%d/%m/%Y")]
feux$heure_fin_inter[nchar(feux$fin_inter) > 10] = substr(feux$fin_inter[nchar(feux$fin_inter) > 10], 12, 16)
feux$heure_fin_inter[!(nchar(feux$fin_inter) > 10)] = ""
feux$heure_fin_inter[feux$heure_fin_inter == "00:00"] = ""
feux$annee = year(feux$date_alerte)
feux$month = month(feux$date_alerte)

Pour voir un peu nos données, rien de tel que des statistiques descriptives. On peut commencer avec du barplot du nombre de feux de forêts par année

agg_annee = aggregate(alerte~annee, FUN=length, data=feux)
colnames(agg_annee)[colnames(agg_annee)=="alerte"] = "Nombre"
nb = ggplot(agg_annee) +
  geom_bar(aes(x=annee, y=Nombre, fill="Feux de forêts"), stat="identity", width = 0.8) +
  scale_colour_manual(name=NULL, values=c("orange")) +
  ggtitle("Nombre de feux par an") +
  theme(plot.title = element_text(hjust = 0.5, size=14),
        legend.position = "none",
        axis.text.x = element_text(angle=50, vjust=0.5, size=10),
        axis.text.y = element_text(size=12),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  scale_fill_discrete("") +
  scale_x_discrete(limit=seq(1981, 2017, by = 2))

ou un barplot de la surface de feux de forêts par anee

agg_surface = aggregate(surface~annee, FUN=sum, data=feux)
surf = ggplot(agg_surface) +
  geom_bar(aes(x=annee, y=surface, fill="Feux de forêts"), stat="identity", width = 0.8) +
  scale_colour_manual(name=NULL, values=c("orange")) +
  ggtitle("Surface totale incendiée par an (en ha)") +
  theme(plot.title = element_text(hjust = 0.5, size=14),
        legend.position = "none",
        axis.text.x = element_text(angle=50, vjust=0.5, size=10),
        axis.text.y = element_text(size=12),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  scale_fill_discrete("") +
  scale_x_discrete(limit=seq(1981, 2017, by = 2)) +
  scale_y_discrete(limit=seq(0, 25000, by = 5000))

Et on fait le graphique

require(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
define_region = function(row, col){
  viewport(layout.pos.row = row, layout.pos.col = col)
} 
print(nb, vp = define_region(1, 1))
print(surf, vp = define_region(1, 2))

Passons aux choses serieuses, avec des cartes (la partie la plus serieuse, on en conviendra, c’etait la phase initiale de scraping). On commence par des informations relatives a la Corse

url="https://raw.githubusercontent.com/data-PM/Promethee/master/communes_corse-shp.zip"
temp="corse.zip"
download.file(url, "corse.zip")
unzip("corse.zip")
communes_corse = sp::spTransform(shapefile("/Users/arthur/communes_corse-shp/communes_corse.shp"), CRS("+proj=longlat +datum=WGS84"))

On se limite aux feux de plus de 1000m2

feux_min1000m2 = feux[feux$surface > 0.1, ]

Maintenant, on peut faire la carte

feux_INSEE = aggregate(x = feux_min1000m2$code_INSEE, by = list(unique.values = feux_min1000m2$code_INSEE), FUN = length)
colnames(feux_INSEE)[colnames(feux_INSEE)=="unique.values"] = "insee"
communes_corse = merge(x = communes_corse, y = feux_INSEE, by = "insee", all.x = TRUE)
communes_corse$x[is.na(communes_corse$x)] = 0
max = max(communes_corse$x)
communes_corse$x2 = communes_corse$x/max * 100
corse2 = gSimplify(communes_corse, tol = 0.001, topologyPreserve = TRUE)
corse2 = SpatialPolygonsDataFrame(corse2, communes_corse@data, match.ID = FALSE)
opar = par(mar = c(0,0,1.2,0))
plot(corse2, col = "#DAE3E6", border = "#8A0641", lwd = 0.7, bg = "#B5D0D0")

Le fond de carte (avec les communes) ressemble a ca

On peut maintenant passer a la carte choroplethe

choroLayer(spdf = corse2, df = corse2@data, var = "x", border = NA, 
           col = carto.pal("wine.pal", 6), legend.pos = "topleft", 
           add = TRUE, method = "fisher-jenks", nclass = 6, legend.title.txt = NULL)
top_villes = corse2[corse2$nom %in% c("Ajaccio", "Bastia", "Corte", "Porto-Vecchio", "Propriano", "Bonifacio", "Calvi"), ]
graphics::points(coordinates(top_villes), pch = 20, cex = 0.5)
labelLayer(top_villes, top_villes@data, txt = "nom", cex = 1, pos = 2, font = 4, offset = 0.5)
layoutLayer(title = "Nombre de feux de forêts depuis 1981 (minimum 1000 m2)",
            source = "Source : 1973-2018 Prométhée ©",
            author = NULL, scale = 10, frame = FALSE, col = "#cdd2d4", coltitle = "#8A5543")

Comme souvent avec R, plusieurs formats sont envisageables. On peut passer a du ggplot2, mais on commence par aller chercher un shapefile global de la Corse

corse = shp_zip("https://raw.githubusercontent.com/data-PM/Promethee/master/corse-shp.zip", "corse")
corse = spTransform(corse, CRS("+proj=longlat +datum=WGS84"))
corse = fortify(corse)
p=data %>%
  ggplot() +
  geom_polygon(data=corse, aes(x=long, y = lat, group = group), fill="grey", alpha=0.3) +
  geom_point(data=feux_min100ha, aes(x=x, y=y, size=surface, color=surface, text=mytext, alpha=surface) ) +
  labs(title = "\nPrincipaux feux de forêts depuis 1981 \n (plus de 100 Ha)") +
  xlim(8.5, 9.6) + ylim(41.3,43.1) +
  scale_size_continuous(range=c(1,5), ) +
  scale_color_viridis(option="inferno", trans="log", name="Surface", limits=c(99, 5700), breaks=c(100, 200, 500, 1000, 2500, 5000)) +
  scale_alpha_continuous(trans="log") +
  theme_void() +
  theme(legend.title = element_text(size=10, face="bold"))
t=800
div(ggplotly(tooltip="text", width=t/1.7, height=t), align = "center")

Joli, non ?

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…

Faire ses cartes météo

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